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 (n×P) storing P biomarker measurements on a set of n observations. Missing values (NA ) are allowed.
|
quantities |
A vector of length n storing the food quantities allocated to each of the n observations in the intervention study data. Missing values (NA ) are not allowed.
|
niter |
The number of MCMC iterations. The default value is niter = 10000 .
|
burnIn |
A numerical value, the number of iterations of the chain to be discarded when computing the posterior estimates. The default value is burnIn = 3000 .
|
posteriors |
A logical value indicating if the full parameter chains should also be returned in output. The default value is posteriors = FALSE .
|
sigmaAlpha |
Intercepts' hyperparameter (σα2 ), see details. The default value is sigmaAlpha = 1 .
|
nuZ1 , nuZ2 |
Are two vectors of length D storing hyperparameters for the components' variance parameters. The default values are nuZ1 = nuZ2 = NULL , corresponding to nuZ1 = (D, D-1, ...,1) and nuZ2 = (D, D,...,D) .
|
nuSigmaP1 , nuSigmaP2 |
Scalar hyperparameters for the error's variance parameters. The default values are nuSigmaP1 = nuSigmaP2 = NULL , corresponding to nuSigmaP1 = 1 and nuSigmaP2 = n .
|
sigmaWprior |
A scalar corresponding to the components' weights hyperparameter. The default value is sigmaWprior = 0.1
|
nuBeta1 , nuBeta2 |
Scalar hyperparameters for the scaling coefficient's variance parameters. The default values are nuBeta1 = 2 and nuBeta2 = 3 .
|
tauBeta |
A scalar factor for the scaling coefficient's variance parameters. The default value is tauBeta = 0.1 .
|
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 Y
of dimension (n×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 X=(X1,…,Xd,…,XD)
and Xd<Xd+1
.
We assume that the biomarker measurements are related to an unobserved, continuous intake value, leading to the following factor analytic model:
yip=αp+βpzi+ϵip,∀i=1,…,n,p=1,…,P,
where the latent variable zi
denotes the latent intake of observation i
, with z=(z1,…,zi,…,zn)
.
The αp
and β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 (μα,σα2)
and (μβ,σβ2)
respectively. The error term ϵp
is the variability associated with biomarker p
. We assume that these errors are normally distributed with 0 mean and variance σp2
, 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 Xd
, and component-specific variances θd2
model food quantity-specific intake variability, with lower values suggesting higher consumption-compliance. Mixture weights are observation-specific and denoted with πi=(πi1,…,π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 95% credible interval lower and upper bounds for the model parameters:
ALPHA_E is a matrix of dimension (4×P) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95% credible interval lower (3rd row) and upper bounds (4th row) for the P intercept parameters, (α1,…,αP) .
BETA_E is a matrix of dimension (4×P) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95% credible interval lower (3rd row) and upper bounds (4th row) for the P scaling coefficient parameters, (β1,…,βP) .
SigmaErr_E is a matrix of dimension (4×P) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95% credible interval lower (3rd row) and upper bounds (4th row) for the P error variance parameters, (σ12,…,σP2) .
SigmaD_E is a matrix of dimension (4×D) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95% credible interval lower (3rd row) and upper bounds (4th row) for the D components' variance parameter, (σ12,…,σD2) .
Z_E is a matrix of dimension (4×n) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95% credible interval lower (3rd row) and upper bounds (4th row) for the n latent intakes, (z1,…,zn) .
THETA_Est is an array of ((P+1)×(D−1)×4) dimensions composed of 4 ((P+1)×(D−1)) matrices, storing the posterior estimates of medians (1st matrix), standard deviations (2nd matrix) and 95% credible interval lower (3rd matrix) and upper bounds (4th matrix) for the components' weights parameters. In each matrix, the first row reports the values for the components' weights intercept parameter, while the other P rows store those of the weights scaling coefficient parameters, (γ,θ1,…,θD−1) .
sigmaBeta_E is a vector containing the posterior estimates of medians, standard deviations and 95% credible interval lower and upper bounds for the scaling coefficients' variance parameter (σβ2 ).
muAlpha_E is a vector containing the posterior estimates of medians, standard deviations and 95% credible interval lower and upper bounds for the intercepts' mean parameter (μα ).
muBeta_E is a vector containing the posterior estimates of medians, standard deviations and 95% credible interval lower and upper bounds for the scaling coefficients' mean parameter (μβ ).
varPHp Estimated error variance parameter values, (νP1∗,νP2∗) , see References.
|
constants |
A list with 11 components, storing constant model quantities:
nuZ1, nuZ2 are two vectors of length D storing hyperparameters for the components' variance parameters, see References.
sigmaAlpha is a scalar and it corresponds to the variance of the intercept parameters (σα2 ).
nuSigmaP1, nuSigmaP2 are scalar hyperparameters for the error's variance parameters, see References.
nuBeta1, nuBeta2 are scalar hyperparameters for the scaling coefficient's variance parameters, see References.
tauBeta is a scalar factor for the scaling coefficient's variance parameters, see References.
x_D is a vector storing the values for the D food quantities.
P is a scalar indicating the number of biomarkers in the data.
D is a scalar indicating the number of food quantities in the data.
n is a scalar indicating the number of observations the data.
sigmaWprior is a scalar and it corresponds to the components' weights hyperparameter, see References.
y_Median is a vector storing the observed P biomarker median values.
y_Var is a vector storing the observed P biomarker variance values.
|
chains |
If posteriors = TRUE , a list with posterior distributions of model parameters is returned:
ALPHA_c is a matrix of dimension (niter−burnIn)×P containing the estimated posterior distributions for the intercept parameters, (α1,…,αP) .
BETA_c is a matrix of dimension (niter−burnIn)×P containing the estimated posterior distributions for the scaling coefficient parameters, (β1,…,βP) .
SigmaErr_c is a matrix of dimension (niter−burnIn)×P containing the estimated posterior distributions for the error variance parameters, (σ12,…,σP2) .
SigmaD_c is a matrix of dimension (niter−burnIn)×D containing are the estimated posterior distributions for the components' variance parameters, (σ12,…,σD2) .
Z_c is a matrix of dimension (niter−burnIn)×n containing the estimated posterior distributions for the latent intakes, (z1,…,zn) .
THETA_c is an array of (P+1)×(D−1)×(niter−burnIn) dimensions containing the estimated posterior distributions for the components' weights parameters. The first one corresponds to that of the weights intercept parameter, while the other P posterior distributions are those of the weights scaling coefficient parameters. In each one of the (niter−burnIn) matrices, the first row reports the values for the components' weights intercept parameter, while the other P rows store those of the weights scaling coefficient parameters, (γ,θ1,…,θD−1) .
sigmaBeta_c is a vector containing the estimated posterior distribution for the scaling coefficients' variance parameter (σβ2 ).
muAlpha_c is a vector containing the estimated posterior distribution for the intercepts' mean parameter (μα ).
muBeta_c is a vector containing the estimated posterior distribution for the scaling coefficients' mean parameter (μβ ).
weights_info is a list containing the acceptance probability values for the weights' parameters, (γ,θ1,…,θD−1) .
|
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)
[Package
multiMarker version 1.0.1
Index]