multiMarker: A latent variable model to infer food intake from multiple...
In multiMarker: Latent Variable Model to Infer Food Intake from Multiple Biomarkers
Description
Usage
Arguments
Details
Value
References
Examples
Description
Implements the multiMarker model via an MCMC algorithm.
Usage
1
2
3
4
5
6
Arguments
y
A matrix of dimension (n x 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 x 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=(X_1, ..., X_d, ..., 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} = α_p + β_p z_i +ε_{ip}, for i=1,...,n, and p = 1, ...,P,
where the latent variable z_i denotes the latent intake of observation i, with (z_1, ..., z_i, .., z_n).
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 σ_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 θ_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}). 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 x 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 x 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 x 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, (σ_1^2, ..., σ_P^2).
SigmaD_E is a matrix of dimension 4 x 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, (σ_1^2, ..., σ_D^2).
Z_E is a matrix of dimension 4 x 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, (z_1, ..., z_n).
THETA_Est is an array of ((P+1) x (D-1) x 4) dimensions composed of 4 ((P+1) x (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) x P containing the estimated posterior distributions for the intercept parameters, (α_1, ..., α_P).
BETA_c is a matrix of dimension (niter-burnIn) x P containing the estimated posterior distributions for the scaling coefficient parameters, (β_1, ..., β_P).
SigmaErr_c is a matrix of dimension (niter-burnIn) x P containing the estimated posterior distributions for the error variance parameters, (σ_1^2, ..., σ_P^2).
SigmaD_c is a matrix of dimension (niter-burnIn) x D containing are the estimated posterior distributions for the components' variance parameters, (σ_1^2, ..., σ_D^2).
Z_c is a matrix of dimension (niter-burnIn) x n containing the estimated posterior distributions for the latent intakes, (z_1, ..., z_n).
THETA_c is an array of (P+1) x (D-1) x (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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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)
multiMarker documentation built on Dec. 15, 2020, 5:36 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value References Examples
Description
Implements the multiMarker model via an MCMC algorithm.
Usage
1 2 3 4 5 6 |
Arguments
y |
A matrix of dimension (n x P) storing P biomarker measurements on a set of n observations. Missing values ( |
quantities |
A vector of length n storing the food quantities allocated to each of the n observations in the intervention study data. Missing values ( |
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 (σ_{α^2}), see details. The default value is |
nuZ1, nuZ2 |
Are two vectors of length D storing hyperparameters for the components' variance parameters. The default values are |
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 Y of dimension (n x 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=(X_1, ..., X_d, ..., 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} = α_p + β_p z_i +ε_{ip}, for i=1,...,n, and p = 1, ...,P,
where the latent variable z_i denotes the latent intake of observation i, with (z_1, ..., z_i, .., z_n). 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 σ_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 θ_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}). 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:
|
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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | 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)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.