Estimate a Bayesian bivariate hierarchical model fitted within INLA.
meta4diag(data=NULL, model.type = 1,
var.prior = "Invgamma", var2.prior = "Invgamma", cor.prior = "Normal",
var.par = c(0.25, 0.025), var2.par, cor.par = c(0,5),
wishart.par = c(4, 1, 1, 0),
init = c(0.01,0.01,0), link="logit", quantiles=c(0.025,0.5,0.975),
modality = NULL, covariates = NULL,
verbose = FALSE, nsample=FALSE,num.threads = 1, seed=0L)A data frame contains at least 4 columns specifying the number of True Positive(TP), False Negative(FN), True Negative(TN) and False Positive(FP). The additional columns other than studynames will be considered as potential covariates and the name or the column number of the potential covariates can be set in the arguments modality and covariates to use them in the model.
A numerical value specifying the model type, options are 1(default), 2, 3 and 4. model.type=1 indicates that the Sensitivity(se) and Specificity(sp) will be modelled in the bivariate model, i.e. \(g(se)\) and \(g(sp)\) are bivariate normal distributed. model.type=2,3,4 indicate that the Sensitivity(se) and False Negative Rate(1-sp), False Positive Rate(1-se) and Specificity(sp), False Positive Rate(1-se) and False Negative Rate(1-sp) are modelled in the bivariate model, respectively.
A string specifying the prior density for the first variance component, options are "PC" for penalised complexity prior, "Invgamma" for inverse gamma prior, "Tnormal" for truncated normal prior, "Unif" for uniform prior which allow the standard deviation uniformaly distributed on [0,1000], "Hcauchy" for half-cauchy prior and "table" for user specific prior. var.prior can also be set to "Invwishart" for inverse wishart prior for covariance matrix. When var.prior="Invwishart", no matter what var2.prior and cor.prior are given, the inverse Wishart prior covariance matrix is used for covariance matrix and the wishart.par must be given. The definition of the priors is as following,
var.prior="Invgamma": This is a prior for a variance \(\sigma^2\). The inverse gamma prior has density,
$$\pi(\sigma^2)=\frac{1}{\Gamma(a)b^a}(\sigma^2)^{-a-1}exp(-\frac{1}{b\sigma^2}),$$
for \(\sigma^2>0\) where: \(a>0\) is the shape parameter, \(b>0\) is the rate (1/scale) parameter.
The parameters are here c(a, b), see arguments var.par.
var.prior="Tnormal": This is a prior for a variance \(\sigma^2\) and defined as follows.
The standard deviation \(\sigma=\sqrt{\sigma^2}\) is Gaussian distributed with mean \(m\) and variance \(v\) but truncated to be positive.
The parameters are here c(m, v), see arguments var.par.
var.prior="PC": This is a prior for a variance \(\sigma^2\) and defined as follows.
The left tail of the distribution of standard deviation \(\sigma\) has behavior
$$P(\sigma>u)=\alpha,$$
which means it is unlikely that the standard deviation \(\sigma\) to be larger than a value \(u\) with probability \(\alpha\).
The parameters are here c(u, alpha), see arguments var.par.
var.prior="HCauchy": This is a prior for a variance \(\sigma^2\) and defined as follows.
The standard deviation \(\sigma=\sqrt{\sigma^2}\) is half-Cauchy distributed with density,
$$\pi(\sigma)=\frac{2\gamma}{\pi(\sigma^2+\gamma^2)},$$
where \(\gamma>0\) is the rate parameter.
The parameters are here c(gamma), see arguments var.par.
var.prior="Unif": This is a prior for a variance \(\sigma^2\) and defined as follows.
The standard deviation \(\sigma=\sqrt{\sigma^2}\) is uniform distributed on \((0,\infty)\).
No parameters need to be given for this prior, see arguments var.par.
var.prior="Table": This is a prior for a variance \(\sigma^2\) and defined as follows.
Users have to specify a data.frame with 2 columns,
one indicates the values of \(\sigma^2\) and the other one indicates the values of \(\pi(\sigma^2)\).
The parameters are this data frame, see arguments var.par.
var.prior="Invwishart": Instead of specifying separate prior distributions for the variance and correlation parameters we could also assume that the covariance matrix \(\Sigma\)
$$\Sigma \sim Wishart^{-1}_{p}(\nu,R),p=2,$$
where the Wishart distribution has density
$$\pi(\Sigma)=\frac{|R|^{\frac{\nu}{2}}}{2^{\frac{p\nu}{2}}\Gamma_{p}(\frac{\nu}{2})}|\Sigma|^{-\frac{\nu+p+1}{2}}exp(-\frac{1}{2}Trace(\frac{R}{\Sigma})), \nu>p+1,$$
Then,
$$E(\Sigma)=\frac{R}{\nu-p-1}.$$
The parameters are \(\nu, R_{11},R_{22} and R_{12}\), where
$$R=\left(\begin{array}{cc}R_{11} & R_{12} \\ R_{21} & R_{22}\end{array}\right)$$
The parameters are here c(nu, R11, R22, R12), see arguments var.par.
See var.prior.
A string specifying the prior density for the correlation, options are "PC" for penalised complexity prior, "Invgamma" for inverse gamma prior, "beta" for beta prior and "table" for user specific prior. cor.prior can also be set to "Invwishart" for inverse wishart prior for covariance matrix. When cor.prior="Invwishart", no matter what var.prior and var2.prior are given, the inverse Wishart prior for covariance matrix is used and the wishart.par must be given. The definition of the priors is as following,
cor.prior="Normal": This is a prior for a correlation \(\rho\) and defined as follows.
The correlation parameter is constrained to \([-1, 1]\). We reparameterise the correlation parameter \(\rho\) using Fisher's z-transformation as
$$\rho^{\star}=logit(\frac{\rho+1}{2}),$$
which assumes values over the whole real line and assign the following prior distribution to \(\rho\),
$$\rho \sim Gaussian(\mu,\sigma^2).$$
The prior variance of \(2.5\) and prior mean of \(0\) corresponds, roughly, to a uniform prior on \([-1,1]\) for \(\rho\) .
The parameters are here c(mean, variance), see arguments cor.par.
cor.prior="PC": This is a prior for a correlation \(\rho\) and defined as follows.
The prior is defined around at a reference point with value \(\rho_{0}\).
To define the density behavior, three strategy can be applied.
The first strategy is to define the left tail behavior and the density weight on the left-hand side of the reference point \(\rho_{0}\),
$$P(\rho<u_{1})=a_{1},$$ and $$P(\rho<\rho_{0})=\omega,$$
which means it is unlikely that the value of \(\rho\) is smaller than a small value \(u_{1}\) with probability \(a_{1}\) and
the probability that \(\rho\) is smaller than \(\rho_0\) is \(\omega\).
The parameters for the first strategy are here c(1, rho0, omega, u1, a1, NA, NA), see arguments cor.par.
The second strategy is to define the right tail behavior and the density weight on the left-hand side of the reference point \(\rho_{0}\),
$$P(\rho>u_{2})=a_{2},$$ and $$P(\rho<\rho_{0})=\omega,$$
which means it is unlikely that the value of \(\rho\) is larger than a big value \(u_{2}\) with probability \(a_{2}\) and
the probability that \(\rho\) is smaller than \(\rho_0\) is \(\omega\).
The parameters for the second strategy are here c(2, rho0, omega, NA, NA, u2, a2), see arguments cor.par.
The third strategy is to define both tail behaviors,
$$P(\rho<u_{1})=a_{1},$$ and $$P(\rho>u_{2})=a_{2}.$$
The parameters for the third strategy are here c(3, rho0, NA, u1, a1, u2, a2), see arguments cor.par.
The parameters of the PC prior for the correlation here is c(strategy, rho0, omega, u1, a1, u2, a2), see arguments cor.par.
cor.prior="beta": This is a prior for a correlation \(\rho\) and defined as follows.
The correlation parameter \(\rho\) has a \(Beta(a,b)\) distribution scaled to have domain in \((-1, 1)\):
$$\pi(\rho)=0.5\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\rho^{a-1}(1-\rho)^{b-1}$$,
where \(a,b>0\) are the shape parameter.
The parameters are here c(a, b), see arguments cor.par.
cor.prior="Table": This is a prior for a correlation \(\rho\) and defined as follows.
Users have to specify the data.frame with 2 columns, one indicates the values of \(\rho\) and the other one indicates the values of \(\pi(\rho)\).
The parameters are this data frame, see arguments cor.par.
cor.prior="Invwishart": Instead of specifying separate prior distributions for the hyperparameters we could also assume that the covariance matrix \(\Sigma\)
$$\Sigma \sim Wishart^{-1}_{p}(\nu,R),p=2,$$
where the inverse Wishart distribution has density
$$\pi(\Sigma)=\frac{|R|^{\frac{\nu}{2}}}{2^{\frac{p\nu}{2}}\Gamma_{p}(\frac{\nu}{2})}|\Sigma|^{-\frac{\nu+p+1}{2}}exp(-\frac{1}{2}Trace(\frac{R}{\Sigma})), \nu>p+1,$$
Then,
$$E(\Sigma)=\frac{R}{\nu-p-1}.$$
The parameters are \(\nu, R_{11},R_{22} and R_{12}\), where
$$R=\left(\begin{array}{cc}R_{11} & R_{12} \\ R_{21} & R_{22}\end{array}\right)$$
The parameters are here c(nu, R11, R22, R12), see arguments cor.par.
A numerical vector specifying the parameter of the prior density for the first variance component.
var.par=c(rate, shape) when var.prior="Invgamma".
var.par=c(u, alpha) when var.prior="PC".
var.par=c(m, v) when var.prior="Tnormal".
var.par=c(gamma) when var.prior="Hcauchy".
var.par=c() when var.prior="Unif".
var.par is a data frame with 2 columns, one indicates the values of \(\sigma^2\) and the other one indicates the values of \(\pi(\sigma^2)\) when var.prior="Table".
var.par doesn't need to be given when var.prior="Invwishart".
See also argument var.prior.
A numerical vector specifying the parameter of the prior density for the second variance component. If not given, function will copy the setting for the first variance component. The definition of the priors is the same as for var.par.
A numerical vector specifying the parameter of the prior density for the correlation. See also examples.
cor.par=c(mean, variance) when cor.prior="normal".
cor.par=c(strategy, rho0, omega, u1, a1, u2, a2) when cor.prior="PC".
cor.par=c(a, b) when var.prior="beta".
cor.par is a data frame with 2 columns, one indicates the values of \(\rho\) and the other one indicates the values of \(\pi(\rho)\) when cor.prior="Table".
cor.par doesn't need to be given when cor.prior="Invwishart".
See also argument cor.prior.
A numerical vector specifying the parameter of the prior density for the covariance matrix. wishart.par must be given and wishart.par=c(nu, R11, R22, R12) when any of var.prior, var2.prior or cor.prior is "Invwishart". See also examples.
A numerical vector specifying the initial value of the first variance, the second variance and correlation.
A string specifying the link function used in the model. Options are "logit", "probit" and "cloglog".
A vector of quantiles, p(0), p(1),... to compute for each posterior marginal. The function returns, for each posterior marginal, the values x(0), x(1),... such that $$Prob(X<x)=p.$$ The default value are c(0.025, 0.5, 0.975). Not matter what other values are going to be given, the estimates for these 3 quantiles are always returned.
Boolean (default:FALSE) indicating whether the program should run in a verbose model.
A string specifying the modality variable, which is a categorical variable, such as test threshold. Default value is NULL. See also examples. At the moment, only one categorical covariate variable can be used in the model.
A vector, which could be either the name of columns or the number of columns, specifying the continuous covariates variables, such as disease prevalence or average individual patients status of each study. Default value is NULL. See also examples.
A numerical value specifying the number of posterior samples, default is 5000. The posterior samples are used to compute the marginals and estimates values of non-linear functions, such as log radios and diagnostic odds radios. If nsample is given, summary.summarized.statistics, summary.fitted.LRpos, summary.fitted.LRneg, summary.fitted.DOR and samples of \(E(se)\), \(E(sp)\), \(E(1-se)\) and \(E(1-sp)\) will be returned.
Maximum number of threads the inla-program will use. xFor Windows this defaults to 1, otherwise its defaults to NULL (for which the system takes over control).
A numerical value specifying the random seed to control the RNG for generating posterior samples if nsample > 0. If you want reproducible results, you ALSO need to control the seed for the RNG in R by controlling the variable .Random.seed or using the function set.seed.
meta4diag returns a meta4diag object with components:
The provided input data.
The internal data that could be used in INLA from function makeData().
Prior distributions for the variance components and correlation from function makePriors().
Names of the jointly modelled accuracies in the model. For example, se and sp or (1-se) and sp.
The cpu time used for running the model.
The matched call.
Matrix containing the mean and standard deviation (plus, possibly quantiles) of the fixed effects of the model.
A list containing the posterior marginal densities of the fixed effects of the model.
Matrix containing the mean and standard deviation (plus, possibly quantiles) of the mean of accuracies transformed with the link function, i.e. E(g(Se)), E(g(Sp)), E(g(1-Se)) and E(g(1-Sp)). Dynamic name for this output. (...) indicates the name of link function used in runModel(), i.e. if link function is "logit", the full name of this output is "summary.expected.logit.accuracy".
A list containing the posterior marginal densities of the mean of accuracies transformed with the link function, i.e. E(g(Se)), E(g(Sp)), E(g(1-Se)) and E(g(1-Sp)). Dynamic name for this output. (...) indicates the name of link function used in runModel(), i.e. if link function is "logit", the full name of this output is "marginals.expected.logit.accuracy".
Matrix containing the mean and standard deviation (plus, possibly quantiles) of the mean of the accuracies, i.e. E(Se), E(Sp), E(1-Se) and E(1-Sp).
A list containing the posterior marginal densities of of the mean of the accuracies, i.e. E(Se), E(Sp), E(1-Se) and E(1-Sp).
A matrix containing the mean and sd (plus, possibly quantiles) of the hyperparameters of the model.
A list containing the posterior marginal densities of the hyperparameters of the model.
A correlation matrix between the mean of the accuracies transformed with the link function. Dynamic name for this output. (...) indicates the name of link function used in runModel().
A covariance matrix between the mean of the accuracies transformed with the link function. Dynamic name for this output. (...) indicates the name of link function used in runModel().
A matrix containing the mean and sd (plus, possibly quantiles) of the linear predictors one transformed accuracy in the model. The accuracy type depends on the model type. See argument model.type. For example, the possible accuracy type could be \(g(se)\), \(g(sp)\), \((se)\) or \((sp)\), where \(g()\) is the link function.
A list containing the posterior marginals of the linear predictors of one transformed accuracy in the model. The accuracy type depends on the model type. See argument model.type. For example, the possible accuracy type could be \(g(se)\), \(g(sp)\), \((se)\) or \((sp)\), where \(g()\) is the link function.
Some other settings that maybe useful retruned by meta4diag.
The deviance information criteria and effective number of parameters.
A list of three elements: cpo$cpo are the values of the conditional predictive ordinate (CPO), cpo$pit are the values of the
probability integral transform (PIT) and cpo$failure
indicates whether some assumptions are violated. In short, if
cpo$failure[i] > 0 then some assumption is violated, the higher the
value (maximum 1) the more seriously.
A list of two elements: waic$waic is the Watanabe-Akaike information criteria, and waic$p.eff is the estimated effective number of parameters.
The log marginal likelihood of the model
A INLA object that from function runModel() which implements INLA.
A matrix of the fixed effects samples if nsample is given.
A matrix of the hyperparameter samples if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of overall sensitivity samples if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of overall specificity samples if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of mean positive and negative likelihood ratios and mean diagnostic odds ratios if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of study specific sensitivity samples if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of study specific specificity samples if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of positive likelihood ratios for each study if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of negative likelihood ratios for each study if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of diagnostic odds ratios for each study if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of risk difference for each study if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of log diagnostic odds ratios for each study if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of log positive likelihood ratios for each study if nsample is given.
A matrix containing the mean and sd (plus, possibly quantiles) of log negative likelihood ratios for each study if nsample is given.
The bivariate model has two levels, in the first level, the observed number of individuals in a specific group in a 2 by 2 table is binomial distributed, for example, the numbers of individuals in the group of true positive and true negative of a study \(i\) are modelled jointly,
$$TP_{i} | Se_{i} \sim Binomial(TP_{i}+ FN_{i},Se_{i})$$
$$TN_{i} | Sp_{i} \sim Binomial(TN_{i}+ FP_{i},Sp_{i})$$
In the second level, two transformed accuracies with some link function (see argument link) are bivariate Gaussian distributed, continuous with the previous example,
$$g(Se_{i}) = \mu + V_{i}\alpha + \phi_{i}$$
$$g(Sp_{i}) = \nu + U_{i}\beta + \psi_{i}$$
where \(\phi_{i}\) and \(\psi_{i}\) are bivariate Gaussian distributed with mean 0 and covariance matrix \(\Sigma\). The \(se\) and \(sp\) in the example could be changed to \(se\) and \((1-sp)\), \((1-se)\) and \(sp\) or \((1-se)\) and \((1-sp)\), see argument model.type.
The function meta4diag() depends on four internal functions which are also given in the meta4diag package in order to flexibly implement a list of dataset with the same prior setting. The four internal functions are makeData(), makePriors(), runModel() and makeObject(). Details can be seen the corresponding documentations and examples.
After running the function meta4diag(), a meta4diag object will be returned which contains various estimated results for later analysis, such as the posterior marginals, estimated value, standard deviations and the coresponding quaniles of the accuracies. See Values.
Rue H., Martino S, and Chopin N. (2009). Approximate Bayesian Inference for Latent Gaussian Models Using Integrated Nested Laplace Approximations. Journal of the Royal Statistical Society B 71: 319--392. (www.r-inla.org)
Simpson DP, Martins TG, Riebler A, Fuglstad GA, Rue H, Sorbye SH (2014) Penalised Model Component Complexity: A principled, Practical Approach to Constructing Priors. Arxiv e-prints. 1403.4630
Guo, J., Riebler, A. and Rue H. (2017) Bayesian bivariate meta-analysis of diagnostic test studies with interpretable priors. Statistics in Medicine 36(19): 3039--3058.
makeData, makePrior, runModel, forest, SROC, crosshair.
# NOT RUN {
if(requireNamespace("INLA", quietly = TRUE)){
require("INLA", quietly = TRUE)
data(Catheter)
meta4diag(data = Catheter, model.type = 1, var.prior = "invgamma", cor.prior = "normal",
var.par = c(0.25, 0.025), cor.par = c(0, 5), init = c(0.01, 0.01, 0),
link = "logit", quantiles = c(0.025, 0.5, 0.975), verbose = FALSE, covariates = NULL,
nsample = FALSE)
}
# }
Run the code above in your browser using DataLab