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)(n\times P) storing PP biomarker measurements on a set of nn observations. Missing values (NA) are allowed.

quantities

A vector of length nn storing the food quantities allocated to each of the nn 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\sigma_{\alpha^2}), see details. The default value is sigmaAlpha = 1.

nuZ1, nuZ2

Are two vectors of length DD 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\mathbf{Y} of dimension (n×P)(n\times P), storing PP different biomarker measurements on nn independent observations. The number of food quantities considered in the intervention study is denoted by DD, with the corresponding set being X=(X1,,Xd,,XD)\mathbf{X}=(X_1, \dots, X_d, \dots, X_D) and Xd<Xd+1X_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:

yip=αp+βpzi+ϵip,i=1,,n,p=1,,P, 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 ziz_i denotes the latent intake of observation ii, with z=(z1,,zi,,zn)\mathbf{z}=(z_1, \dots, z_i, \dots, z_n). The αp\alpha_p and βp\beta_p parameters characterize, respectively, the intercept and the scaling effect for biomarker pp. We assume that these parameters are distributed a priori according to 0-truncated Gaussian distributions, with parameters (μα,σα2)(\mu_{\alpha}, \sigma_{\alpha}^2) and (μβ,σβ2)(\mu_{\beta}, \sigma_{\beta}^2 ) respectively. The error term ϵp\epsilon_p is the variability associated with biomarker pp. We assume that these errors are normally distributed with 0 mean and variance σp2\sigma_p^2, which serves as a proxy for the precision of the biomarker.

A mixture of DD 0-truncated Gaussian distributions is assumed as prior distribution for the latent intakes. Components are centered around food quantity values XdX_d, and component-specific variances θd2\theta_d^2 model food quantity-specific intake variability, with lower values suggesting higher consumption-compliance. Mixture weights are observation-specific and denoted with πi=(πi1,,πiD)\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 95%95\% credible interval lower and upper bounds for the model parameters:

  • ALPHA_E is a matrix of dimension (4×P)(4\times P) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95%95\% credible interval lower (3rd row) and upper bounds (4th row) for the PP intercept parameters, (α1,,αP)(\alpha_1, \dots, \alpha_P).

  • BETA_E is a matrix of dimension (4×P)(4\times P) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95%95\% credible interval lower (3rd row) and upper bounds (4th row) for the PP scaling coefficient parameters, (β1,,βP)(\beta_1, \dots, \beta_P).

  • SigmaErr_E is a matrix of dimension (4×P)(4\times P) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95%95\% credible interval lower (3rd row) and upper bounds (4th row) for the PP error variance parameters, (σ12,,σP2)(\sigma_1^2, \dots, \sigma_P^2).

  • SigmaD_E is a matrix of dimension (4×D)(4\times D) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95%95\% credible interval lower (3rd row) and upper bounds (4th row) for the DD components' variance parameter, (σ12,,σD2)(\sigma_1^2, \dots, \sigma_D^2).

  • Z_E is a matrix of dimension (4×n)(4\times n) storing the posterior estimates of medians (1st row), standard deviations (2nd row) and 95%95\% credible interval lower (3rd row) and upper bounds (4th row) for the nn latent intakes, (z1,,zn)(z_1, \dots, z_n).

  • THETA_Est is an array of ((P+1)×(D1)×4)((P+1)\times (D-1) \times 4) dimensions composed of 44 ((P+1)×(D1))((P+1)\times (D-1)) matrices, storing the posterior estimates of medians (1st matrix), standard deviations (2nd matrix) and 95%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 PP rows store those of the weights scaling coefficient parameters, (γ,θ1,,θD1)(\gamma, \theta_1, \dots,\theta_{D-1}).

  • sigmaBeta_E is a vector containing the posterior estimates of medians, standard deviations and 95%95\% credible interval lower and upper bounds for the scaling coefficients' variance parameter (σβ2\sigma_{\beta}^2).

  • muAlpha_E is a vector containing the posterior estimates of medians, standard deviations and 95%95\% credible interval lower and upper bounds for the intercepts' mean parameter (μα\mu_{\alpha}).

  • muBeta_E is a vector containing the posterior estimates of medians, standard deviations and 95%95\% credible interval lower and upper bounds for the scaling coefficients' mean parameter (μβ\mu_{\beta}).

  • varPHp Estimated error variance parameter values, (νP1,νP2)(\nu_{P1}^{*}, \nu_{P2}^{*}), see References.

constants

A list with 11 components, storing constant model quantities:

  • nuZ1, nuZ2 are two vectors of length DD storing hyperparameters for the components' variance parameters, see References.

  • sigmaAlpha is a scalar and it corresponds to the variance of the intercept parameters (σα2\sigma_{\alpha^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 DD 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 PP biomarker median values.

  • y_Var is a vector storing the observed PP biomarker variance values.

chains

If posteriors = TRUE, a list with posterior distributions of model parameters is returned:

  • ALPHA_c is a matrix of dimension (niterburnIn)×P(niter-burnIn)\times P containing the estimated posterior distributions for the intercept parameters, (α1,,αP)(\alpha_1, \dots, \alpha_P).

  • BETA_c is a matrix of dimension (niterburnIn)×P(niter-burnIn)\times P containing the estimated posterior distributions for the scaling coefficient parameters, (β1,,βP)(\beta_1, \dots, \beta_P).

  • SigmaErr_c is a matrix of dimension (niterburnIn)×P(niter-burnIn)\times P containing the estimated posterior distributions for the error variance parameters, (σ12,,σP2)(\sigma_1^2, \dots, \sigma_P^2).

  • SigmaD_c is a matrix of dimension (niterburnIn)×D(niter-burnIn)\times D containing are the estimated posterior distributions for the components' variance parameters, (σ12,,σD2)(\sigma_1^2, \dots, \sigma_D^2).

  • Z_c is a matrix of dimension (niterburnIn)×n(niter-burnIn)\times n containing the estimated posterior distributions for the latent intakes, (z1,,zn)(z_1, \dots, z_n).

  • THETA_c is an array of (P+1)×(D1)×(niterburnIn) (P+1)\times(D-1)\times(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 PP posterior distributions are those of the weights scaling coefficient parameters. In each one of the (niterburnIn)(niter-burnIn) matrices, the first row reports the values for the components' weights intercept parameter, while the other PP rows store those of the weights scaling coefficient parameters, (γ,θ1,,θD1)(\gamma, \theta_1, \dots,\theta_{D-1}).

  • sigmaBeta_c is a vector containing the estimated posterior distribution for the scaling coefficients' variance parameter (σβ2\sigma_{\beta}^2).

  • muAlpha_c is a vector containing the estimated posterior distribution for the intercepts' mean parameter (μα\mu_{\alpha}).

  • muBeta_c is a vector containing the estimated posterior distribution for the scaling coefficients' mean parameter (μβ\mu_{\beta}).

  • weights_info is a list containing the acceptance probability values for the weights' parameters, (γ,θ1,,θD1)(\gamma, \theta_1, \dots,\theta_{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]