nma.ab: Arm-Based Network Meta-Analysis for Binary Outcomes

Description

nma.ab performs arm-based network meta-analysis proposed by Zhang et al (2014). It estimates absolute risk (AR), risk difference (RD), odds ratio (OR), relative risk (RR), log odds ratio (LOR), and log relative risk (LRR).

Usage

nma.ab(s.id, t.id, event.n, total.n, data, trtname,
       param = c("AR", "LOR", "LRR", "RD", "rank.prob"),
       model = "het_cor", 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 = "", dic = FALSE, postdens = FALSE)

Arguments

s.id

a numeric or character vector indicating study ID, or the corresponding column name in the argument data.

t.id

a numeric or character vector indicating treatment ID, or the corresponding column name in the argument data.

event.n

a numeric vector of non-negative integers, indicating event number for a certain treatment in the corresponding study, or the corresponding column name in the argument data.

total.n

a numeric vector of positive integers, indicating total number of participants for a certain treatment in the corresponding study, or the corresponding column name in the argument data.

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 correspond

trtname

a vector of character strings indicating the treatment names for the corresponding treatment IDs according their order in t.id. If not specified, t.id is used as treatment names.

param

a vector of character strings indicating the effect sizes to be estimated. The default includes "AR" (absolute risk), "LOR" (log odds ratio), "LRR" (log relative risk), "RD" (risk difference), "ran

model

a character string indicating which Bayesian hierarchical model to be applied in the arm-based network meta-analysis. This argument can be set as "hom_ind", "het_ind", "hom_eqcor", "het_eqcor", or

prior.type

prior distribution of variances/covariances of random effects. If model is "hom_eqcor", "het_eqcor", "hom_ind", or "het_ind", it can be set as "unif" (uniform prior for standard

a, b

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", "het_eqcor", "hom_ind"

positive number, specifying the upper bound of uniform prior for standard deviations of random effects if using prior.type as "unif" for model "hom_eqcor", "het_eqcor", "hom_ind", or "

higher.better

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 higher event rate implying better treatment

digits

a positive integer specifying the digits after the decimal point for the effect size estimates. The default is 4.

n.adapt

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

n.iter

the total number of iterations in each MCMC chain. The default is 100,000.

n.burnin

the number of iterations for burn-in period. The default is n.iter/2.

n.chains

the number of MCMC chains. The default is 3.

n.thin

a positive integer indicating thinning rate. The default is the thinning rate which yields no more than 100,000 iterations remaining in each chain.

conv.diag

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 re

trace

a vector of character strings of effect sizes. The character strings should be selected from those specified in param (except "rank.prob"), and trace plots of the specified effect sizes will be saved in users' current working dir

dic

a logical value indicating whether the deviance information criterion (DIC) to be calculated. The default is FALSE.

postdens

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<

Value

nma.ab returns a list with estimates of effect sizes 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.

Details

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. $\mu_{k}$'s are treatment-specific fixed effects, and random effects $\nu_{ik}$'s are correlated within each study with the covariance matrix $\mathbf{\Sigma}_{K}$. Using an inverse-Wishart prior for $\mathbf{\Sigma}_{K}$, the above model corresponds to model = "het_cor". Denote $\sigma_{k}$ as the standard deviation of $\nu_{ik}$ and $\mathbf{D} = diag(\sigma_{1}, \ldots, \sigma_{K})$, then the correlation matrix $\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". Moreover, if $\mathbf{R}_{K}$ is set to be the $K \times K$ identity matrix, then the model becomes "het_ind" for heterogeneous variances of random effects, and "hom_ind" for homogeneous variance. Note that the models "het_ind" and "hom_ind" lead to independence between different treatments. Therefore, they may break study randomization, and we do not recommend these two models. In addition, for the models "hom_eqcor", "het_eqcor", "hom_ind" and "het_ind", setting prior.type as "invgamma" implies using inverse-gamma priors with shape and scale parameters $a, b$ for $\sigma_{k}^2$ or $\sigma^2$, and "unif" implies uniform priors $U(0, c)$ for $\sigma_{k}$ or $\sigma$.

References

Gelman A and Rubin, DB (1992). "Inference from iterative simulation using multiple sequences." Stat Sci 7(4), 457--472. Lu G and Ades AE (2004). "Combination of direct and indirect evidence in mixed treatment comparisons." Stat Med 23(20), 3105--24. Spiegelhalter DJ, Best NG, Carlin BP, and Van Der Linde A (2002). "Bayesian measures of model complexity and fit." J R Stat Soc Series B Stat Methodol 64(4), 583--639. 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." Clin Trials 11(2), 246--62.

Examples

Run this code

data(smoke)

# For the smoke cessation data,
# higher event rate indicates better treatment

# use the model = "het_cor"
#het.cor.out <- nma.ab(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
hom.eqcor.out <- nma.ab(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)

Run the code above in your browser using DataLab