version 1.0

cran · Jul 2, 2020 · 93ba593 · 93ba593
commit 93ba593
Show file tree

Hide file tree

Showing 10 changed files with 1,106 additions and 0 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -0,0 +1,25 @@
+Package: multiMarker
+Version: 1.0
+Date: 2020-05-29
+Title: Latent Variable Model to Infer Food Intake from Multiple
+        Biomarkers
+Description: A latent variable model based on factor analytic and mixture of experts models, designed to infer food intake from multiple biomarkers data. The model is framed within a Bayesian hierarchical framework, which provides flexibility to adapt to different biomarker distributions and facilitates prediction of the intake along with its associated uncertainty. Details are in D'Angelo, et al. (2020) <arXiv:2006.02995>.
+Authors@R: c( person("Silvia", "D\'Angelo", role = c("aut", "cre"), 
+           email = "silvia.dangelo@ucd.ie"), 
+	         person("Claire", "Gormley", role = "ctb"),
+	         person("Lorraine", "Brennan", role = "ctb")
+	         )
+Maintainer: Silvia D'Angelo <silvia.dangelo@ucd.ie>
+Imports: truncnorm, ordinalNet
+Depends: R (>= 3.0)
+License: GPL (>= 2)
+Encoding: UTF-8
+ByteCompile: true
+LazyData: true
+NeedsCompilation: no
+Packaged: 2020-06-09 15:09:16 UTC; Silvia
+Author: Silvia D'Angelo [aut, cre],
+  Claire Gormley [ctb],
+  Lorraine Brennan [ctb]
+Repository: CRAN
+Date/Publication: 2020-07-02 11:10:07 UTC
diff --git a/MD5 b/MD5
@@ -0,0 +1,9 @@
+b48e2883799779d93d73fbf26f171bb0 *DESCRIPTION
+23463272c0d14b1e504b892d3b6e8f43 *NAMESPACE
+6ec04e6729dbadd0d772b82471e0c2e3 *R/internal_functions.R
+121d2c22bce282eb42ff0690ed47435f *R/multiMarker.R
+90666372ab6c650e685ea33db50c2678 *R/predict.multiMarker.R
+70df53897e7538451f65245e5a2647e4 *R/zzz_figlet.R
+1946fec930ca1761f396bd7e85e33200 *man/internal_functions.Rd
+83de2c5550f37eb58a0d7aeaff91af8e *man/multiMarker.Rd
+ff9fd482375c0980959aa55cba711c73 *man/predict.multiMarker.Rd
diff --git a/NAMESPACE b/NAMESPACE
@@ -0,0 +1,14 @@
+# Export everything
+# exportPattern("^[^\\.]")
+
+export(multiMarker)
+
+importFrom("stats", "rgamma", "runif", "sd", "var",
+           "coef", "lm", "median", "quantile", "rnorm")
+importFrom("truncnorm", "rtruncnorm")
+importFrom("ordinalNet", "ordinalNet")
+importFrom("utils", "flush.console")
+importFrom("utils", "setTxtProgressBar", "txtProgressBar")
+
+S3method("print", "multiMarker")
+S3method(predict, multiMarker)
diff --git a/R/internal_functions.R b/R/internal_functions.R
@@ -0,0 +1,217 @@
+##--- utilities ---##
+msdci <- function(p_chain){ # after burn in
+
+  out <- c( median(p_chain), sd(p_chain),
+            as.numeric( quantile(p_chain, c(0.025, 0.975), na.rm = TRUE) ) )
+  return(out)
+}
+
+
+##############################
+## initialization functions ##
+##############################
+
+alpha_beta_iniz <- function( x_D, y, D, n_D, P){
+
+  tmpN <- c(1,cumsum(n_D))
+  meanY_PD <-  sapply(1:P, function(p) # P cols and D rows
+    sapply(1:D, function(d) 
+      mean(y[tmpN[d]:tmpN[d+1], p], na.rm = TRUE)
+    ))
+
+  ncolM <- (D-1)*D - sum( seq(1,D-1) )
+  A_Mat <- B_Mat <- matrix( nrow = P, ncol = ncolM )
+
+  # computing alpha and beta values 
+  for( p in 1:P){
+    a <- 1
+    for( d in D:2){
+      for( j in (d-1):1){
+        alpha <- max(0, ( meanY_PD[d,p] - meanY_PD[j,p]*(x_D[d]/x_D[j]) )/( 1 - (x_D[d]/x_D[j])) )
+        beta <- (meanY_PD[j,p] - alpha)/x_D[j]
+        if( (beta <= 0) | is.nan(beta)){beta <- 0.001}
+        A_Mat[p, a] <- alpha
+        B_Mat[p, a] <- beta
+        a <- a+1
+      }
+    }
+  } # close P loop (rows are the Ps)
+
+  alpha_iniz <- rowMeans(A_Mat, na.rm = TRUE)
+  beta_iniz <- rowMeans(B_Mat, na.rm = TRUE)
+
+  out <- list( alpha_iniz =  alpha_iniz,  beta_iniz =  beta_iniz)
+  return(out)
+}
+
+
+sigma2_err_iniz <- function(beta, y, D, P){
+
+  temp1 <- diag( beta%*%t(beta) )/D
+  temp2 <- diag( var(y) )
+
+  sigma2_err_iniz <- rep( sum( temp2 - temp1, na.rm = TRUE ), P )
+  sigma2_err_iniz[which(sigma2_err_iniz <= 0)] <- 0.3
+
+  out <- list( sigma2_err_iniz  = sigma2_err_iniz )
+  return(out)
+}
+
+z_par_iniz <- function(y, sigma2_err, sigmaAlpha, 
+                       sigmaBeta, P, D, n_D){
+  tmpN <- c(1,cumsum(n_D))
+
+  temp <- matrix(NA, ncol = P, nrow = D)
+  varY_PD <-  sapply(1:P, function(p) # P cols and D rows
+    sapply(1:D, function(d) 
+      var(y[tmpN[d]:tmpN[d+1], p], na.rm = TRUE)
+    ))
+
+  for ( p in 1:P){
+    for( d in 1:D){
+      temp[d,p] <- (varY_PD[d,p] -sigmaAlpha - sigma2_err[p])/(sigmaBeta)
+    }
+  }
+
+  out <- rowMeans(temp)
+  out[which(out <= 0)] <- 0.2
+  return(out)
+}
+
+#################################
+## full conditionals functions ##
+#################################
+
+#-- all nuisance parameters (those of alpha, beta and mix components)--#
+
+variance_fc <- function(param, l_param, nu1, nu2,
+                        tau, mu_param, mu_mu_param){
+  nu1_star <- (l_param )/2 + nu1
+  if( nu1_star <= 1){ nu1_star <- 1}
+  if( is.infinite(nu1_star)){ nu1_star <- 100}
+  nu2_star <- nu2 + (tau * sum( (param -mu_param)^2, na.rm = T ) +
+                       (mu_param - mu_mu_param)^2)/(2 * tau)
+  if( nu2_star <= 1){ nu2_star <- 1}
+  if( is.infinite(nu2_star)){ nu2_star <- 100}
+  out <- 1/rgamma( 1, shape = nu1_star, rate = nu2_star)
+  return(out)
+}
+
+variance_fc_d <- function(param, l_param, nu1, nu2,
+                          tau, mu_param, scale_Fact_d){
+  nu1_star <- (l_param )/2 + nu1
+  if( nu1_star <= 1){ nu1_star <- 1}
+  if( is.infinite(nu1_star)){ nu1_star <- 100}
+  nu2_star <- nu2 + sum((param - mu_param)^2, na.rm = TRUE)/(2*scale_Fact_d)
+  if( nu2_star <= 1){ nu2_star <- 1}
+  if( is.infinite(nu2_star)){ nu2_star <- 100}
+  out <- 1/rgamma( 1, shape = nu1_star, rate = nu2_star)
+  OUT <- list( out = out,  nu1_star =  nu1_star, nu2_star = nu2_star )
+  return(OUT)
+}
+
+variance_fc_di <- function(param, nu1, nu2,
+                           mu_param, p_est_d, n){
+  nu1_star <- nu1 + sum( p_est_d * n ,na.rm = T)/2
+  nu2_star <- nu2 + sum(p_est_d * (param - mu_param)^2, na.rm = TRUE)/(2)
+  out <- 1/rgamma( 1, shape = nu1_star, rate = nu2_star)
+  return(out)
+}
+
+
+
+mean_fc <- function( tau, sigma2, param, l_param,
+                     mu_mu_param){
+  sigma2_star <- (tau * sigma2)/(tau * l_param + 1)
+  mu_star <- (tau * sum(param, na.rm = T) + mu_mu_param)/(tau * sigma2)
+  mu_star <- sigma2_star * mu_star
+  out <- rtruncnorm(1, 0, Inf, mu_star, sqrt(sigma2_star))
+  return(out)
+}
+
+#--main parameters (alpha, beta)--#
+
+alpha_fc <- function( sigma2_a, beta_p, n, mu_a,
+                      z, sigma2_err_p, y_p){
+  sigma2_star <- (sigma2_err_p * sigma2_a)/(n * sigma2_a + sigma2_err_p)
+  mu_star <- (sum(y_p -  beta_p * z, na.rm = T))/(sigma2_err_p) + mu_a/sigma2_a
+  mu_star <- sigma2_star * mu_star
+  out <- rtruncnorm(1, 0, Inf, mu_star, sqrt(sigma2_star)) 
+  OUT <- list( out = out,  mu_star =  mu_star, sigma2_star = sigma2_star )
+  return(OUT)
+}
+
+beta_fc <- function( sigma2_b, alpha_p, n, mu_b,
+                     z, sigma2_err_p, y_p){
+
+  sigma2_star <- (sigma2_err_p * sigma2_b)/( sum(z^2, na.rm = T) *
+                                               sigma2_b + sigma2_err_p)
+  mu_star <- (sum( (y_p -  alpha_p) * z, na.rm = T))/
+    (sigma2_err_p) + mu_b/sigma2_b
+  mu_star <- sigma2_star * mu_star
+  out <- rtruncnorm(1, 0, Inf, mu_star, sqrt(sigma2_star))
+  OUT <- list( out = out, mu_star =  mu_star, sigma2_star = sigma2_star )
+  return(OUT)
+}
+
+#--latent intakes (z) for a given component (d)--#
+
+z_fc <- function( sigma2_d, mu_d, sigma2_err, 
+                  beta, y_i, alpha, P, scale_Fact_d){
+
+  sigma2_d <- sigma2_d * scale_Fact_d
+  sigma2_star <- sum( ((beta^2) * sigma2_d + 
+                         sigma2_err / P)/(sigma2_d * sigma2_err), na.rm = T )
+  sigma2_star <- (sigma2_star)^(-1)  
+  mu_star <- sum((beta * (y_i - alpha)) / sigma2_err, na.rm = T) +
+    (mu_d / sigma2_d)
+  mu_star <- sigma2_star * mu_star
+
+  sigma2_star <-  sigma2_star
+  out <- rtruncnorm( 1,0, Inf, mu_star, sqrt(sigma2_star))
+  return(out)  
+}
+
+#--error variances (sigma_p)--#
+
+sigma2_err_fc <- function( n, theta1, theta2, y_p,
+                           alpha_p, beta_p, z){
+  par1 <- (n / 2) + theta1
+  par2 <- 0.5 * sum( (y_p - alpha_p -(beta_p * z))^2, na.rm = T ) +theta2
+  out <- 1/rgamma(1, shape = par1 , rate = par2)
+  OUT <- list( out = out,  nu1_star = par1, nu2_star = par2)
+  return(OUT)
+}
+
+#--log post MCMC weigths--#
+logPost_MCMC_wcum <- function(prob_nc, D, labelsMat, z, doses){ # non cum probs
+  out <- sum( labelsMat * log(prob_nc), na.rm = T ) 
+  out <- out - sum(abs(z-doses))
+
+  return(out)
+}
+
+cauchit_probs <- function(yNew, theta, D){ # Ynew is for a single obs i
+  scaling_y <- yNew * theta[-1,1]  # vector length P
+  tmp <- sapply(1:(D-1), function(d)
+    (1/pi)*( (atan(theta[1,d] + sum(scaling_y, na.rm = T)) + pi/2) )
+  ) # length D-1
+  tmp2 <- c( tmp[1], diff(tmp) )
+  TMP <- c( tmp2, 1 - sum(tmp2,na.rm = T))
+  if (sum(c(TMP) < 1) == T ){ TMP <- -Inf}
+  return(TMP) # output is a vector of length D
+}
+
+
+dim_range_prob <- function(prob_c, D){
+  opz1 <- which(prob_c >=0.5)
+  if( length(opz1) != 0 ){ 
+    out <- opz1}else{
+      tmp2 <- cbind(seq(1,D), prob_c)
+      tmp2 <- tmp2[order(tmp2[,2], decreasing = T),]
+      tmp3 <- which(cumsum(tmp2[,2]) >=0.5)[1]
+      out <- tmp2[1:tmp3,1]
+    }
+  return(out)
+}
+