Performs the arm-based network meta-analysis proposed by Zhang et al. (2014), including estimating the absolute risk (AR), risk difference (RD), odds ratio (OR), risk ratio (RR), log odds ratio (LOR), and log risk ratio (LRR).
nma.ab.bin(s.id, t.id, event.n, total.n, data, trtname,
param = c("AR", "LOR", "LRR", "RD", "rank.prob"),
model = "het_cor", link = "probit",
prior.type, a = 0.001, b = 0.001, c = 10,
higher.better = FALSE, digits = 4, n.adapt = 5000,
n.iter = 100000, n.burnin = floor(n.iter/2), n.chains = 3,
n.thin = max(1, floor((n.iter - n.burnin)/100000)),
conv.diag = FALSE, trace = NULL, dic = FALSE, postdens = FALSE,
mcmc.samples = FALSE)nma.ab.bin returns a list with estimates of effect measures specified in param. If the argument dic = TRUE, the deviance information criterion (DIC) statistic will be returned in the output list. In addition, if conv.diag = TRUE, a .txt file containing the point estimates of the potential scale reduction factor and their upper confidence limits by Gelman and Rubin (1992) will be saved in users' current working directory. If postdens = TRUE, the posterior densities of treatment-specific absolute risks will be saved as a .pdf file. If trace is specified, the trace plots are saved as .png files.
a numeric or character vector specifying study ID, or the corresponding column name in the argument data.
a numeric or character vector specifying treatment ID, or the corresponding column name in the argument data.
a numeric vector of non-negative integers, specifying event number for a certain treatment in the corresponding study, or the corresponding column name in the argument data.
a numeric vector of positive integers, specifying total number of participants for a certain treatment in the corresponding study, or the corresponding column name in the argument data.
an optional data frame containing the dataset for network meta-analysis. If data is specified, the previous arguments, s.id, t.id, event.n, and total.n, should be specified as the corresponding column names in data; otherwise, the previous arguments use environment variables.
a vector of character strings specifying the treatment names for the corresponding treatment IDs according to their order in t.id. If not specified, t.id is used as treatment names.
a vector of character strings specifying the effect measures to be estimated. The default includes "AR" (absolute risk), "LOR" (log odds ratio), "LRR" (log risk ratio), "RD" (risk difference), "rank.prob" (treatment rank probability). "AR" is automatically added into param even if it is not specified. In addition to the default measures, "OR" (odds ratio) and "RR" (risk ratio) can be added into param. If the logit link is used (link = "logit"), "OR.med" (overall median odds ratio), "LOR.med" (overall median log odds ratio), and "rank.prob.med" (treatment rank probability based on overall median odds ratio) can be also added; see "Details". If model is "hom_eqcor" or "het_eqcor", param can also include "\rho" (the common correlation coefficient between treatments).
a character string specifying which Bayesian hierarchical model to be applied in the arm-based network meta-analysis. This argument can be set as "hom_eqcor", "het_eqcor", or "het_cor" (the default). See "Details" for the models.
a character string specifying the link function in the Bayesian hierarchical model for binary outcomes. It can be either "probit" (the default) or "logit".
prior distribution of variances and/or covariances of random effects. If model is "hom_eqcor" or "het_eqcor", it can be set as "unif" (uniform prior for standard deviations, the default) or "invgamma" (inverse gamma prior for variances). If model is "het_cor", prior.type can be "invwishart" (the default) or "chol". Specifying "invwishart" yields an inverse-Wishart prior for the variance-covariance matrix of random effects; by specifying "chol", non-informative priors are assigned to variance and correlation components using the separation strategy by Cholesky decomposition. See "Details".
positive numbers, specifying the shape and scale parameters of inverse gamma priors for variance(s) of random effects if using prior.type as "invgamma" for model "hom_eqcor" or "het_eqcor". The default values for both parameters are 0.001.
positive number, specifying the upper bound of uniform prior for standard deviation(s) of random effects if using prior.type as "unif" for model "hom_eqcor" or "het_eqcor". The default is 10.
an optional logical value which needs to be specified when estimating the treatment rank probabilities (i.e., "rank.prob" is included in the argument param). TRUE indicates that a higher event rate implies a better treatment, and vice versa. The default is FALSE.
a positive integer specifying the digits after the decimal point for the effect measure estimates. The default is 4.
the number of iterations for adaptation in Markov chain Monte Carlo (MCMC) algorithm. The default is 5,000. If a warning "adaptation incomplete" appears, users may increase n.adapt. This argument and the following n.iter, n.burnin, n.chains, and n.thin are passed to the functions in R package rjags.
the total number of iterations in each MCMC chain. The default is 100,000.
the number of iterations for burn-in period. The default is n.iter/2.
the number of MCMC chains. The default is 3.
a positive integer specifying the thinning rate. The default is the thinning rate which yields no more than 100,000 iterations remaining in each chain.
a logical value indicating whether to perform MCMC convergence diagnostic. The default is FALSE. If TRUE, n.chains must be greater than 1, and a .txt file, which contains the point estimates of the potential scale reduction factor and their upper confidence limits (see Gelman and Rubin 1992), will be written in users' current working directory.
a vector of character strings of effect measures. The character strings should be selected from those specified in param (except "rank.prob"), and trace plots of the specified effect measures will be saved in users' current working directory. The default is not drawing trace plots (NULL).
a logical value indicating whether to calculate the deviance information criterion (DIC) value. The default is FALSE.
a logical value indicating whether to draw the posterior density plots for treatment-specific absolue risks (ARs). If TRUE, a .pdf file containing the density plot will be saved in users' current working directory. The default is FALSE.
a logical value indicating whether to save MCMC posterior samples in the output object. The default is FALSE.
Suppose that a network meta-analysis collects \(I\) studies on \(K\) treatments, where each study investigates a subset of the \(K\) treatments. Label the studies from \(i = 1\) to \(I\) and the treatments from \(k = 1\) to \(K\). Let \(T_{i}\) be the subset of the \(K\) treatments that is compared in the \(i\)th study. Also, in the \(i\)th study, let \(n_{ik}\) be the number of participants allocated to treatment group \(k\) (\(k \in T_{i}\)), and \(y_{ik}\) be the number of events. The arm-based model is constructed as (Zhang et al. 2014): $$y_{ik} \sim Bin (n_{ik}, p_{ik}) \qquad k \in T_{i};$$ $$\Phi^{-1} (p_{ik}) = \mu_{k} + \nu_{ik};$$ $$(\nu_{i1}, \nu_{i2}, \ldots, \nu_{iK})^{T} \sim N (\boldsymbol{0}, \mathbf{\Sigma}_{K}),$$ where \(\Phi (\cdot)\) is the standard normal cumulative distribution function, and \(\mathbf{\Sigma}_{K}\) is a \(K \times K\) positive definite covariance matrix. The \(\mu_{k}\)'s are treatment-specific fixed effects, and the random effects \(\nu_{ik}\) are correlated within each study with the covariance matrix \(\mathbf{\Sigma}_{K}\). The marginal absolute risk of treatment \(k\) is \(p_{k} = E[p_{ik}]\); other effect measures are calculated based on these absolute risks.
An unstructured covariance matrix \(\mathbf{\Sigma}_{K}\) in the model above corresponds to model = "het_cor". The inverse-Wishart prior can be assigned to \(\mathbf{\Sigma}_{K}\). Alternatively, using the separation strategy by Cholesky decomposition (prior.type = "chol"), uniform priors \(U(0, c)\) are assigned to the standard deviations in \(\mathbf{\Sigma}_{K}\) and non-informative priors are assigned to the correlation components (Barnard et al., 2000; Lu and Ades, 2009; Wei and Higgins, 2013; Lin and Chu, 2018). Denote \(\sigma_{k}\) as the standard deviation of \(\nu_{ik}\) and \(\mathbf{D} = diag(\sigma_{1}, \ldots, \sigma_{K})\), then the correlation matrix is \(\mathbf{R}_{K} = \mathbf{D}^{-1} \mathbf{\Sigma}_{K} \mathbf{D}^{-1}\). If we assume that all of the off-diagonal elements in \(\mathbf{R}_{K}\) are equal, say to \(\rho\), then this model corresponds to model = "het_eqcor". If we further assume the homogeneity of variances of the random effects, that is, \(\sigma_{k} = \sigma\) for \(k = 1, 2, \ldots, K\), then the model is "hom_eqcor". In addition, for the models "hom_eqcor" and "het_eqcor", setting prior.type as "invgamma" implies using inverse-gamma priors with shape and scale parameters, \(a\) and \(b\), for \(\sigma_{k}^2\) or \(\sigma^2\), and "unif" implies uniform priors \(U(0, c)\) for \(\sigma_{k}\) or \(\sigma\).
In addition to the probit link as used in the model above, one may also use the logit link (Chu et al., 2012), which is adopted more commonly in contrast-based models. Using the logit link, \(\mu_k\) represents the overall median log odds across studies; thus, \(\mu_k - \mu_h\) is the overall median log odds ratio between treatments \(h\) and \(k\), \(\exp(\mu_k - \mu_h)\) is the overall median odds ratio, and they futher yield the rank probabilities. The contrast-based models usually report such overall median (log) odds ratios, rather than the marginal effect measures obtained from the arm-based models. When the argument link is "logit", users can additionally specify "LOR.med", "OR.med", and "rank.prob.med" to obtain these results based on the overall median odds ratios.
Barnard J, McCulloch R, Meng XL (2000). "Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage." Statistica Sinica, 10(4), 1281--1311.
Chu H, Nie L, Chen Y, Huang Y, Sun W (2012). "Bivariate random effects models for meta-analysis of comparative studies with binary outcomes: methods for the absolute risk difference and relative risk." Statistical Methods in Medical Research, 21(6), 621--633. <tools:::Rd_expr_doi("10.1177/0962280210393712")>
Gelman A, Rubin DB (1992). "Inference from iterative simulation using multiple sequences." Statistical Science, 7(4), 457--472. <tools:::Rd_expr_doi("10.1214/ss/1177011136")>
Lin L, Chu H (2018). "Bayesian multivariate meta-analysis of multiple factors." Research Synthesis Methods, 9(2), 261--272. <tools:::Rd_expr_doi("10.1002/jrsm.1293")>
Lin L, Zhang J, Hodges JS, Chu H (2017). "Performing arm-based network meta-analysis in R with the pcnetmeta package." Journal of Statistical Software, 80(5), 1--25. <tools:::Rd_expr_doi("10.18637/jss.v080.i05")>
Lu G, Ades AE (2004). "Combination of direct and indirect evidence in mixed treatment comparisons." Statistics in Medicine, 23(20), 3105--3124. <tools:::Rd_expr_doi("10.1002/sim.1875")>
Lu G, Ades AE (2009). "Modeling between-trial variance structure in mixed treatment comparisons." Biostatistics, 10(4), 792--805. <tools:::Rd_expr_doi("10.1093/biostatistics/kxp032")>
Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A (2002). "Bayesian measures of model complexity and fit." Journal of the Royal Statistical Society, Series B (Statistical Methodology), 64(4), 583--639. <tools:::Rd_expr_doi("10.1111/1467-9868.00353")>
Wei Y, Higgins JPT (2013). "Bayesian multivariate meta-analysis with multiple outcomes." Statistics in Medicine, 32(17), 2911--2934. <tools:::Rd_expr_doi("10.1002/sim.5745")>
Zhang J, Carlin BP, Neaton JD, Soon GG, Nie L, Kane R, Virnig BA, Chu H (2014). "Network meta-analysis of randomized clinical trials: Reporting the proper summaries." Clinical Trials, 11(2), 246--262. <tools:::Rd_expr_doi("10.1177/1740774513498322")>
nma.ab.cont, nma.ab.py, nma.ab.followup
data("smoke")
# For the smoke cessation data,
# higher event rate indicates better treatment
# use the model = "het_cor"
#set.seed(1234)
#het.cor.out <- nma.ab.bin(s.id, t.id, r, n, data = smoke,
# trtname = c("NC", "SH", "IC", "GC"), param = c("AR", "OR", "RR", "LOR",
# "LRR", "RD", "rank.prob"), model = "het_cor", higher.better = TRUE,
# n.iter = 200000, n.thin = 1, conv.diag = TRUE, dic = TRUE,
# trace = c("AR", "LOR"), postdens = TRUE)
# use the model = "hom_eqcor"
# increase n.iter to reach convergence
set.seed(123)
hom.eqcor.out <- nma.ab.bin(s.id, t.id, r, n, data = smoke,
param = c("AR", "LRR"), model = "hom_eqcor", prior.type = "unif", c = 10,
higher.better = TRUE, n.adapt = 1000, n.iter = 100, n.chains = 1)
# use the logit link
set.seed(1234)
hom.eqcor.out2 <- nma.ab.bin(s.id, t.id, r, n, data = smoke,
param = c("AR", "OR", "OR.med", "rank.prob", "rank.prob.med"),
model = "hom_eqcor", link = "logit", prior.type = "unif", c = 10,
higher.better = TRUE, n.adapt = 1000, n.iter = 100, n.chains = 1)
Run the code above in your browser using DataLab