meta4diag {meta4diag}

Function to analyse diagnostic meta-analysis with Bayesian methods using INLA.

Description

Estimate a Bayesian bivariate hierarchical model fitted within INLA.

Usage

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)

Arguments

data

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.

model.type

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.

var.prior

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.

var2.prior

See var.prior.

cor.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.

var.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.

var2.par

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.

cor.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.

wishart.par

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.

init

A numerical vector specifying the initial value of the first variance, the second variance and correlation.

link

A string specifying the link function used in the model. Options are "logit", "probit" and "cloglog".

quantiles

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.

verbose

Boolean (default:FALSE) indicating whether the program should run in a verbose model.

modality

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.

covariates

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.

nsample

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.

num.threads

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).

seed

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.

Details

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.

Value

meta4diag returns a meta4diag object with components:

data

The provided input data.

outdata

The internal data that could be used in INLA from function makeData().

priors.density

Prior distributions for the variance components and correlation from function makePriors().

names.fitted

Names of the jointly modelled accuracies in the model. For example, se and sp or (1-se) and sp.

cpu.used

The cpu time used for running the model.

call

The matched call.

summary.fixed

Matrix containing the mean and standard deviation (plus, possibly quantiles) of the fixed effects of the model.

marginals.fixed

A list containing the posterior marginal densities of the fixed effects of the model.

summary.expected.(...).accuracy

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".

marginals.expected.(...).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".

summary.expected.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).

marginals.expected.accuracy

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).

summary.hyperpar

A matrix containing the mean and sd (plus, possibly quantiles) of the hyperparameters of the model.

marginals.hyperpar

A list containing the posterior marginal densities of the hyperparameters of the model.

correlation.expected.(...).accuracy

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().

covariance.expected.(...).accuracy

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().

summary.predictor.(...)

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.

marginals.predictor.(...)

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.

misc

Some other settings that maybe useful retruned by meta4diag.

dic

The deviance information criteria and effective number of parameters.

cpo

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.

waic

A list of two elements: waic$waic is the Watanabe-Akaike information criteria, and waic$p.eff is the estimated effective number of parameters.

mlik

The log marginal likelihood of the model

inla.result

A INLA object that from function runModel() which implements INLA.

samples.fixed

A matrix of the fixed effects samples if nsample is given.

samples.hyperpar

A matrix of the hyperparameter samples if nsample is given.

samples.overall.Se

A matrix containing the mean and sd (plus, possibly quantiles) of overall sensitivity samples if nsample is given.

samples.overall.Sp

A matrix containing the mean and sd (plus, possibly quantiles) of overall specificity samples if nsample is given.

summary.overall.statistics

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.

samples.study.specific.Se

A matrix containing the mean and sd (plus, possibly quantiles) of study specific sensitivity samples if nsample is given.

samples.study.specific.Sp

A matrix containing the mean and sd (plus, possibly quantiles) of study specific specificity samples if nsample is given.

summary.study.specific.LRpos

A matrix containing the mean and sd (plus, possibly quantiles) of positive likelihood ratios for each study if nsample is given.

summary.study.specific.LRneg

A matrix containing the mean and sd (plus, possibly quantiles) of negative likelihood ratios for each study if nsample is given.

summary.study.specific.DOR

A matrix containing the mean and sd (plus, possibly quantiles) of diagnostic odds ratios for each study if nsample is given.

summary.study.specific.RD

A matrix containing the mean and sd (plus, possibly quantiles) of risk difference for each study if nsample is given.

summary.study.specific.LDOR

A matrix containing the mean and sd (plus, possibly quantiles) of log diagnostic odds ratios for each study if nsample is given.

summary.study.specific.LLRpos

A matrix containing the mean and sd (plus, possibly quantiles) of log positive likelihood ratios for each study if nsample is given.

summary.study.specific.LLRneg

A matrix containing the mean and sd (plus, possibly quantiles) of log negative likelihood ratios for each study if nsample is given.

Author(s)

Jingyi Guo and Andrea Riebler

References

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.

See Also

makeData, makePrior, runModel, forest, SROC, crosshair.

Examples

## 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)
}

## End(Not run)