multiMarker {multiMarker} | R Documentation |
A latent variable model to infer food intake from multiple biomarkers.
Description
Implements the multiMarker model via an MCMC algorithm.
Usage
multiMarker(y, quantities,
niter = 10000, burnIn = 3000,
posteriors = FALSE, sigmaAlpha = 1,
nuZ1 = NULL, nuZ2 = NULL,
nuSigmaP1 = NULL, nuSigmaP2 = NULL, sigmaWprior = 0.000001,
nuBeta1 = 2, nuBeta2 = 3, tauBeta = 0.1)
Arguments
y |
A matrix of dimension |
quantities |
A vector of length |
niter |
The number of MCMC iterations. The default value is |
burnIn |
A numerical value, the number of iterations of the chain to be discarded when computing the posterior estimates. The default value is |
posteriors |
A logical value indicating if the full parameter chains should also be returned in output. The default value is |
sigmaAlpha |
Intercepts' hyperparameter ( |
nuZ1 , nuZ2 |
Are two vectors of length |
nuSigmaP1 , nuSigmaP2 |
Scalar hyperparameters for the error's variance parameters. The default values are |
sigmaWprior |
A scalar corresponding to the components' weights hyperparameter. The default value is |
nuBeta1 , nuBeta2 |
Scalar hyperparameters for the scaling coefficient's variance parameters. The default values are |
tauBeta |
A scalar factor for the scaling coefficient's variance parameters. The default value is |
Details
The function facilitates inference of food intake from multiple biomarkers via MCMC, according to the multiMarker model (D'Angelo et al., 2020). The multiMarker model first learns the relationship between the multiple biomarkers and food quantity data from an intervention study and subsequently allows inference on the latent intake when only biomarker data are available.
Consider a biomarker matrix \mathbf{Y}
of dimension (n\times P)
, storing P
different biomarker measurements on n
independent observations.
The number of food quantities considered in the intervention study is denoted by D
, with the corresponding set being \mathbf{X}=(X_1, \dots, X_d, \dots, X_D)
and X_d < X_{d+1}
.
We assume that the biomarker measurements are related to an unobserved, continuous intake value, leading to the following factor analytic model:
y_{ip} = \alpha_p + \beta_p z_i +\epsilon_{ip}, \quad \forall \quad i=1,\dots,n, \quad p = 1, \dots,P,
where the latent variable z_i
denotes the latent intake of observation i
, with \mathbf{z}=(z_1, \dots, z_i, \dots, z_n)
.
The \alpha_p
and \beta_p
parameters characterize, respectively, the intercept and the scaling effect for biomarker p
. We assume that these parameters are distributed a priori according to 0-truncated Gaussian distributions, with parameters (\mu_{\alpha}, \sigma_{\alpha}^2)
and (\mu_{\beta}, \sigma_{\beta}^2 )
respectively. The error term \epsilon_p
is the variability associated with biomarker p
. We assume that these errors are normally distributed with 0 mean and variance \sigma_p^2
, which serves as a proxy for the precision of the biomarker.
A mixture of D
0-truncated Gaussian distributions is assumed as prior distribution for the latent intakes. Components are centered around food quantity values X_d
, and component-specific variances \theta_d^2
model food quantity-specific intake variability, with lower values suggesting higher consumption-compliance. Mixture weights are observation-specific and denoted with \pi_i = (\pi_{i1}, \dots,\pi_{iD} )
. Given the inherent ordering of the food quantities in the intervention study, an ordinal regression model with Cauchit link function is employed to model the observation-specific weights.
A Bayesian hierarchical framework is employed for the modelling process, allowing quantification of the uncertainty in intake estimation, and flexibility in adapting to different biomarker data distributions. The framework is implemented through a Metropolis within Gibbs Markov chain Monte Carlo (MCMC) algorithm. Hyperprior distributions are assumed on the prior parameters with the corresponding hyperparameter values fixed based on the data at hand, following an empirical Bayes approach.
For more details on the estimation of the multiMarker model, see D'Angelo et al. (2020).
Value
An object of class 'multiMarker'
containing the following components:
estimates |
A list with 9 components, storing posterior estimates of medians, standard deviations and
|
constants |
A list with 11 components, storing constant model quantities:
|
chains |
If
|
References
D'Angelo, S. and Brennan, L. and Gormley, I.C. (2020). Inferring food intake from multiple biomarkers using a latent variable model. arxiv
Examples
library(truncnorm)
oldpar <- par(no.readonly =TRUE)
#-- Simulate intervention study biomarker and food quantity data --#
P <- D <- 3; n <- 50
alpha <- rtruncnorm(P, 0, Inf, 4, 1)
beta <- rtruncnorm(P, 0, Inf, 0.001, 0.1)
x <- c(50, 100, 150)
labels_z <- sample(c(1,2,3), n, replace = TRUE)
quantities <- x[labels_z]
sigma_d <- 8
z <- rtruncnorm(n, 0, Inf, x[labels_z], sigma_d)
Y <- sapply( 1:P, function(p) sapply( 1:n, function(i)
max(0, alpha[p] + beta[p]*z[i] + rnorm( 1, 0, 5) ) ) )
#-- Visualize the data --#
par(mfrow= c(2,2))
boxplot(Y[,1] ~ quantities, xlab = "Food quantity", ylab = "Biomarker 1")
boxplot(Y[,2] ~ quantities, xlab = "Food quantity", ylab = "Biomarker 2")
boxplot(Y[,3] ~ quantities, xlab = "Food quantity", ylab = "Biomarker 3")
#-- Fit the multiMarker model --#
# Number of iterations (and burnIn) set small for example.
modM <- multiMarker(y = Y, quantities = quantities,
niter = 100, burnIn = 30,
posteriors = TRUE)
# niter and burnIn values are low only for example purposes
#-- Extract summary statistics for model parameters --#
modM$estimates$ALPHA_E[,3] #estimated median, standard deviation,
# 0.025 and 0.975 quantiles for the third intercept parameter (alpha_3)
modM$estimates$BETA_E[,2] #estimated median, standard deviation,
# 0.025 and 0.975 quantiles for the second scaling parameter (beta_2)
#-- Examine behaviour of MCMC chains --#
par(mfrow= c(2,1))
plot(modM$chains$ALPHA_c[,3], type = "l",
xlab = "Iteration (after burnin)", ylab = expression(alpha[3]) )
abline( h = mean(modM$chains$ALPHA_c[,3]), lwd = 2, col = "darkred")
plot(modM$chains$BETA_c[,2], type = "l",
xlab = "Iteration (after burnin)", ylab = expression(beta[2]) )
abline( h = mean(modM$chains$BETA_c[,2]), lwd = 2, col = "darkred")
# compute Effective Sample Size
# library(LaplacesDemon)
# ESS(modM$chains$ALPHA_c[,3]) # effective sample size for alpha_3 MCMC chain
# ESS(modM$chains$BETA_c[,2]) # effective sample size for beta_2 MCMC chain
par(oldpar)