Performing Arm-Based Network Meta-Analysis in R with the pcnetmeta Package

Network meta-analysis is a powerful approach for synthesizing direct and indirect evidence about multiple treatment comparisons from a collection of independent studies. At present, the most widely used method in network meta-analysis is contrast-based, in which a baseline treatment needs to be specified in each study, and the analysis focuses on modeling relative treatment effects (typically log odds ratios). However, population-averaged treatment-specific parameters, such as absolute risks, cannot be estimated by this method without an external data source or a separate model for a reference treatment. Recently, an arm-based network meta-analysis method has been proposed


Introduction
In diverse scientific fields, such as social and medical research, summaries of cumulative knowledge are increasingly based on the results of meta-analyses (Hunter and Schmidt 1996;Lindholm et al. 2005;Cooper et al. 2009).Meta-analysis is a statistical method for combining and contrasting a collection of estimated effect sizes, such as odds ratios, from multiple independent studies (DerSimonian and Laird 1986).Various approaches for We use an illustrative example to show the difference between the contrast-and arm-based approaches.Suppose that the outcome in a network meta-analysis is binary, and y ik and n ik are the numbers of events and participants, respectively, in treatment group k in the ith study.For such data, both the contrast-and arm-based models use the binomial likelihood y ik ~ Binomial(n ik , p ik ); they differ in the way they model the underlying absolute risks p ik in each study's treatment group.Specifically, the contrast-based method needs to specify a baseline treatment b(i) in the ith study.For convenience, we simply denote b(i) as b.Then, the Bayesian hierarchical model for this approach is (Lu andAdes 2004, 2009): where g(•) is a link function and X ik is a dummy variable taking the value 0 if k = b or 1 if k ≠ b.Also, μ i is the baseline effect for treatment b in the ith study, and δ ibk is the relative effect of treatment k compared with the baseline b on the g-transformed scale.Note that this model treats the μ i 's as nuisances and uses non-informative priors for them.This model is described as contrast-based because it focuses on the overall relative effects d hk between treatment pairs (h, k), which are estimated using the evidence consistency equation d hk = d bk − d bh .This model does not permit a back-transformation from the relative effects to absolute effects, unless the absolute effect of a given "reference" treatment group can be accurately estimated from external data, or can be estimated using a separate model to analyze the existing data in the "reference" treatment group (Welton et al. 2012;Dias et al. 2013b).
Population-averaged absolute effects are preferred in some situations such as costeffectiveness analysis and patient decisions (Dias et al. 2013b).For example, consider two scenarios comparing treatments A and B according to one-year survival rates: (i) p A = 0.8 vs. p B = 0.5; (ii) p A = 0.004 vs. p B = 0.001.Both scenarios yield an odds ratio 4.0, but patients would prefer treatment A in scenario (i) more strongly than in scenario (ii).Therefore, an absolute effect or absolute difference is preferred in this case.Compared with the contrast-based method, the arm-based model provides a straightforward way to estimate absolute effects and various types of relative effects.The model is specified as (Zhang et al. 2014(Zhang et al. , 2015a)): where Σ K is the variance-covariance matrix of the vector of random effects specific to study i.The μ k 's are treatment-specific parameters reflecting absolute effects.Based on the absolute effects, various types of relative effects can be obtained.In some cases (e.g., the missingness of treatment arms is not at random), the effect size estimates produced by the arm-based method may be less biased than those given by the contrast-based method (Zhang et al. 2015a).Moreover, the arm-based method can use information contained in single-arm studies, in which only one treatment group is available or of interest.Single-arm studies cannot be included in the contrast-based model but they may provide valuable information for treatment comparisons and enhance the robustness of a network meta-analysis (Lin et al. 2016a).Although the arm-based approach has many advantages, the convergence of Markov chain Monte Carlo (MCMC) algorithms for parameter estimation may be slower compared with the contrast-based approach.The estimates using the arm-based models may not converge well if some treatments are only included in a few (say, less than three) studies.Also, some researchers have concerns about the arm-based model, for example, that absolute effects tend to be highly variable compared to relative effects, and that pooling arm-level data may not fully respect the randomization process in randomized controlled trials (Dias and Ades 2016).However, these concerns are mainly raised because the arm-and contrastbased models use different assumptions about treatment effects.Specifically, the contrastbased model assumes that relative effects are exchangeable across studies, while the armbased model assumes absolute effects are exchangeable (Hong et al. 2016b).Although the assumption of exchangeable relative effects is popular in meta-analysis, the assumption of exchangeable absolute effects is also accepted in the literature (see, e.g., van Houwelingen et al. 1993;Shuster et al. 2007;Senn 2010;Chu et al. 2012).More details of the arm-based model are discussed in Dias and Ades (2016) and Hong et al. (2016b).
Plenty of software packages are dedicated to conducting traditional meta-analysis (e.g., Rosenberg et al. 2000;Borenstein et al. 2005;Viechtbauer 2014;Schwarzer 2015), but very limited software is available specifically for network meta-analysis.The R package netmeta (Rücker et al. 2015) provides models in a frequentist framework described in Rücker (2012).
The R package gemtc (van Valkenhoef and Kuiper 2015) and the Stata (StataCorp 2015) command network perform contrast-based analyses.Neither package provides estimates for population-averaged treatment-specific parameters.This article introduces the R package pcnetmeta (Lin et al. 2016b), which performs network meta-analysis using the arm-based model and provides estimates for various effect sizes.This package is available from Comprehensive R Archive Network (R Core Team 2015) at http://CRAN.R-project.org/package=pcnetmeta.It uses MCMC techniques on the R platform through JAGS (Plummer 2015b(Plummer , 2016)).JAGS is a program for analyzing Bayesian hierarchical models using MCMC simulation, which is available for diverse computer platforms including Windows and Mac OS X.
The package pcnetmeta provides user-friendly functions to perform network meta-analysis for various types of data.Convergence of the MCMC routine can be assessed by the function outputs.The package also provides functions to draw network plots which illustrate the comparisons between multiple treatments.In addition, plots for 95% credible intervals of treatment-specific and relative effect sizes are provided to visually display treatment effects and their comparisons.This article is organized as follows.Section 2 presents an overview of arm-based Bayesian hierarchical models as implemented in the pcnetmeta package.Section 3 illustrates the use of the package with several examples, and discusses the output structures.Finally, Section 4 closes with suggested future improvements.

Arm-based model for binary outcomes
Suppose a network meta-analysis reviews I studies on K treatments, where each study only investigated a subset of the K treatments.Let T i (1 ≤ i ≤ I) be the set of treatments compared in the ith study.Also, in the ith study, the total number of events/participants in treatment group k (k ∈ T i ) is denoted by y ik /n ik .Zhang et al. (2014) specified the following arm-based model using the probit link function: (1) where Φ(•) denotes the standard normal cumulative distribution function.In this model, μ k is a fixed effect for Treatment k, and the random effects (ν i1 , ν i2 , …, ν iK ) ⊤ are correlated within each study with variance-covariance matrix Σ K .Based on this model, the absolute risk of Treatment k can be estimated as (Zeger et al. 1988), where is the kth diagonal element in Σ K .Since this estimate is a marginal expectation of p ik given μ k and σ k , we can interpret p k as the population-averaged absolute risk of Treatment k.With the estimates p k for absolute risk, we can further estimate the risk difference (RD), odds ratio (OR), and risk ratio (RR), which are defined as RD kl = p k − p l , , and RR kl = p k /p l , respectively.Other link functions may be also considered in the arm-based model, but they do not yield simple expressions for population-averaged absolute risks, though some approximations exist.For example, using the logit link we can approximate , where (Zeger et al. 1988).This article and the package pcnetmeta use the probit link function for simplicity.

Arm-based model for continuous outcomes
Researchers can also perform network meta-analysis on studies with continuous outcomes (e.g., Kasiske et al. 1998;Philbrook et al. 2007;Zhang et al. 2015b).We continue to use i and k to index the studies and treatments and n ik as the total number of participants receiving Treatment k in the ith study.The summary data include sample mean ȳ ik and within-study sample standard deviation s ik .We can specify the arm-based model for continuous outcomes as (2) In this model, μ k is of interest because it can be interpreted as the overall effect of Treatment k by noticing that E[θ ik |μ k , σ k ] = μ k , and we may further estimate the effect difference between Treatments k and l, which is defined as d kl = μ k − μ l .Dias et al. (2013a) discussed methods for count datasets in network meta-analysis.Their models are contrast-based, but we can consider the corresponding arm-based versions.

Arm-based model for count datasets
In some network meta-analyses, the available data are in the form of counts over a certain time period.The total number of person-years at risk is supplied rather than the total number of participants.Let y ik be the number of events in treatment group k in the ith study, and E ik is the corresponding exposure time in person-years.Suppose λ ik is the treatment-specific rate for treatment group k in the ith study, and we are interested in the population-averaged treatment-specific rate.We consider the following arm-based model with a Poisson likelihood and the log link function: (3) A key assumption for this model is that the rates λ ik are constant over the follow-up period.
Based on this model, the population-averaged treatment specific rate can be estimated as .
A similar model is available for studies that report the proportion of patients developing an event within a given follow-up period, where the follow-up time may differ for each study (e.g., Psaty et al. 2003;Elliott and Meyer 2007).In this case, the event probability depends on the length of follow-up.Specifically, we have y ik and n ik as the number of events and participants, respectively, in treatment group k in the ith study, and the study-specific followup times are denoted f i .Following Dias et al. (2013a), we assume a latent Poisson process with rate λ ik for each treatment group in each study; therefore the time T ik until an event occurs is distributed as exponential with rate λ ik and survivor function Pr(T ik > f i ) = exp(−λ ik f i ).Thus, the event probability is p ik = 1 − exp(−λ ik f i ).Using the complementary log-log link cloglog(t) = log(−log(1 − t)) for p ik , we have cloglog(p ik ) = log(f i ) + log(λ ik ).
Again, we can model log(λ ik ) using model (3) above, and an arm-based model is constructed as follows: (4) The parameter of interest is still the population-averaged treatment-specific rate .

Parameter estimation for arm-based models
Let ν i = (ν i1 , ν i2 , …, ν iK ) ⊤ .The full likelihood function for the arm-based model is where L c (y ik ; μ k , ν ik ) is the conditional likelihood given the random effects ν ik .For example, for binary outcomes, the conditional likelihood is where g −1 (•) is the inverse of the link function.To obtain the maximum likelihood estimates, the maximization problem is subject to the condition that Σ K is positive definite.When the outcome is binary or count, the full likelihood function cannot be expressed in a closed form; for continuous outcomes, it may have a closed form.However, maximizing the likelihood may be unstable and converge slowly if the number of treatments K is large and the number of collected studies I is small.Alternatively, we can, and pcnetmeta does, apply MCMC to obtain Bayesian estimates for parameters of interest.In arm-based models, vague N(0, 1000) priors are used for the treatment-specific fixed effects μ k .As suggested in Gelman and Hill (2007), we may assign an inverse-Wishart prior to the unstructured variance-covariance matrix Σ K with the scale matrix being the K × K identity matrix I K and degrees of freedom K + 1, i.e., .This has the effect of setting a uniform prior on the individual correlation parameters.Alternatively, the separation strategy by Cholesky decomposition can be used to specify a vague prior to Σ K (Barnard et al. 2000; Lu and Ades 2009; Wei and  Higgins 2013).We denote this model as HET-COR because the variances of the random effects are heterogeneous.
To reduce model complexity, we may assume an exchangeable correlation structure for Σ K (Zhang et al. 2015a), that is, To guarantee that R ex is positive definite, ρ must be greater than .A vague uniform prior on can be used for ρ.We denote models with this exchangeable correlation matrix as HET-EQCOR for heterogeneous variances σ k .If we further assume homogeneity of variances, that is, σ k = σ for all k = 1, 2, …, K, the model is denoted as HOM-EQCOR.
To avoid overfitting, we can use the deviance information criterion (DIC) proposed by Spiegelhalter et al. (2002) for model selection.A smaller penalized deviance implies a better model.
Finally, to implement the hierarchical models in JAGS, the package pcnetmeta automatically generates initial values for the parameters and specifies different random number generators (RNGs) for different chains.Zero is used to initialize the random effects ν ik ; the initial values for the variance parameters Σ K are the means of their prior distributions.For example, since in the HET-COR model, the initial value for is the mean of this Wishart distribution, (K + 1)I K .For the fixed effects μ k , the initial values are the naïve estimates of the corresponding absolute effects computed by simply pooling the data in each treatment group.For example, consider the continuous outcome in model ( 2), the initial values for μ k 's are Σ {i:k∈T i } ȳ ik n ik /Σ {i:k∈T i } n ik .

Using the R package pcnetmeta
The R package pcnetmeta provides user-friendly functions to perform arm-based network meta-analysis using the models described above.Users can download its source file at http:// CRAN.R-project.org/package=pcnetmeta, or directly install it within R by typing install.packages("pcnetmeta").Note that the pcnetmeta depends on the R packages rjags (Plummer 2016) and coda (Plummer 2015a).The pcnetmeta package does not include a copy of the JAGS library, so users must install JAGS separately.JAGS is freely available at its homepage http://mcmc-jags.sourceforge.net/.Also, the package pcnetmeta requires JAGS version ≥ 4.0.0; the earlier versions of JAGS may not guarantee exact reproducibility of the results.In this section, we introduce the basic usage of this package.

Data structure for network meta-analysis
To begin, we briefly introduce the necessary dataset structures.In the package pcnetmeta, four datasets, smoke, parkinson, dietaryfat, and diabetes, are provided as illustrative examples.
The dataset smoke contains 24 studies on smoking cessation with binary outcomes, reported in Hasselblad (1998) and Lu and Ades (2006).This network meta-analysis compares four treatments, labeled as: 1) no contact (NC); 2) self-help (SH); 3) individual counseling (IC); and 4) group counseling (GC).We display the dataset's first few rows below.The column s.id contains IDs for the 24 studies, and t.id labels the treatments included in each study.For example, Study 1 compares Treatments 1, 3, and 4. The columns r and n are the number of events (successful cessation) and participants, respectively.The datasets dietaryfat and diabetes serve as examples for models (3) and ( 4).This article uses diabetes to illustrate estimation of treatment-specific rates and rate ratios.This dataset was analyzed by Elliott and Meyer (2007) to assess the effects of antihypertensive agents on incident diabetes, and includes the follow-up times (in years) for each study.Twenty-two clinical studies are included, covering six different treatments: 1) diuretic; 2) placebo; 3) beta blocker (BB); 4) calcium-channel blocker (CCB); 5) angiotensinconverting-enzyme inhibitor (ACEI); and 6) angiotensin-receptor blocker (ARB).Users can apply data() in R to load the datasets dietaryfat and diabetes; we do not display the detailed dataset.
Note that NA is not allowed in the dataset for the package pcnetmeta, because the published articles collected in a network meta-analysis typically report all summary results (such as both mean and variance for continuous outcomes).Also, each row in the dataset represents one treatment group in a study, so a single-arm study is straightforwardly input as a single row in the dataset for analysis using the arm-based models in the pcnetmeta package.

Plotting the network
The function nma.networkplot() in package pcnetmeta provides a visual overview of treatment comparisons in network datasets.Calling this function produces a network graph in an R plot window.Each vertex in the network plot represents a treatment and each edge between two nodes stands for a direct comparison between the corresponding two treatments.The usage of the function is as follows: R> args(nma.networkplot)R> nma.networkplot(s.id,t.id, data = diabetes, trtname = c("Diuretic", "Placebo", "BB", "CCB", "ACEI", "ARB")) Figure 1 shows the resulting graphs.In the left panel for the smoking cessation data, every pair of treatment nodes is connected by an edge, so all pairs of treatments are directly compared.For the dataset parkinson in Figure 1b, we did not specify treatment names, so the function used the treatment IDs from 1 to 5 as the names.Its network plot shows no edge between some pairs of treatments, e.g., Treatments 2 and 3.This means that no study directly compares Treatments 2 and 3.As a result, if the aim is to compare the effects of Treatments 2 and 3, only indirect evidence is available, e.g., from the comparisons Treatments 2 vs. 1 and 3 vs. 1. Figure 1c is the network plot for the dataset diabetes; all pairs of treatments are directly compared except ACEI and ARB.

Performing arm-based network meta-analysis
The major functions in package pcnetmeta are nma.ab.bin(), nma.ab.cont(), nma.ab.py(), and nma.ab.followup(), which perform arm-based network metaanalysis for different types of data using the models introduced in Section 2. In particular, nma.ab.bin() analyzes binary outcomes, while nma.ab.cont() is used for continuous outcomes.These two functions are based on models (1) and (2), respectively.Functions nma.ab.py() and nma.ab.followup() can be used when exposure times or follow-up times are available, and are based on the models (3) and ( 4), respectively.The commands in each of the following subsections may take around 5{20 minutes on an Intel 2.60 GHz processor.The actual runtime depends on the complexity of the treatment network and the user's processor.As in nma.networkplot(), the arguments s.id and t.id are numeric or character vectors indicating study and treatment IDs.Users also specify each study's number of events and participants using event.nand total.nrespectively.The argument model can be specified as "het_cor", "het_eqcor", or "hom_eqcor", which corresponds to the models described in Section 2.4.When model = "het_cor" (the default), users can specify prior.type= "invwishart" (the default) to assign an inverse-Wishart prior to the variance-covariance matrix of random effects.Alternatively, by assigning prior.type= "chol", the separation strategy by Cholesky decomposition is used for the variancecovariance matrix, and uniform priors U(0, c) are assigned to the standard deviations and vague priors are assigned to the correlation components (Barnard et al. 2000;Lu and Ades 2009;Wei and Higgins 2013).When "het_eqcor" and "hom_eqcor" are used, the correlation matrix of the random effects has an exchangeable correlation structure.For the "hom_eqcor" and "het_eqcor" models, two types of priors for the random-effect variances can be used.A popular prior for variances is inverse-Gamma(ε, ε) where ε can be set to a low value (e.g., 0.001) (Spiegelhalter et al. 2003).However, the posterior may be sensitive to the choice of ε, and a uniform prior for standard deviation is preferred in some cases (Gelman 2006).Users choose the prior by specifying prior.type= "unif" (the default) or "invgamma", and the prior parameters The argument param is a character string vector which indicates the effect sizes to be estimated.As in Section 2.1, param can include absolute risk ( "AR"), odds ratio ( "OR"), log odds ratio ( "LOR"), risk ratio ( "RR"), log risk ratio ( "LRR"), risk difference ( "RD").In addition, researchers may be interested in treatment ranks (Salanti et al. 2011).Users can estimate the rank probabilities of different treatments (i.e., probabilities of the treatment having ranks 1, 2, …, K) by adding "rank.prob" to the argument param.When "rank.prob" is added, users need to specify the logical argument higher.better;TRUE indicates that a higher event probability implies a better treatment.For example, the event in the dataset smoke is smoking cessation, and a higher smoking cessation probability implies a better treatment.Many outcomes in medical studies are the events of developing a disease (e.g., Thijs et al. 2008), in which case a better treatment should lead to a lower event probability.
The arguments n.adapt, n.iter, n.burnin, n.chains, and n.thin control the MCMC algorithm run by rjags (Plummer 2016).The argument n.adapt is the number of iterations for adaptation (the default is 5,000); this is used to maximize the sampling efficiency.The argument n.iter determines the number of iterations in each MCMC chain, and n.burnin is the number of burn-in iterations at the beginning of each chain which are discarded.The argument n.chains is the number of MCMC chains; the default is 3. Additionally, n.thin is the thinning rate for MCMC chains, which is used to save memory and computation time if n.iter is large.For example, if n.iter is 10 6 and n.thin is 10, then only one sample would be kept in every 10 samples in each chains, and the remaining number of iterations is 10 5 .
The argument conv.diagspecifies whether to compute potential scale reduction factors (PSRFs) proposed by Gelman and Rubin (1992) for convergence diagnostics.The argument dic = TRUE, the function will provide the deviance information criterion (DIC) statistic proposed by Spiegelhalter et al. (2002); conventionally the model with smallest DIC is considered the best among the candidate models.The posterior density plot of treatmentspecific effect sizes can be obtained as a .pdffile by setting postdens = TRUE.
The function nma.ab.bin() returns a list with effect size estimates, which lists the posterior mean, standard deviation, median, and a 95% credible interval (CI) with 2.5% and 97.5% quantiles as the lower and upper bounds.
Here is an example to demonstrate the function's usage.We called the function nma.ab.bin() on the dataset smoke as follows: R> data("smoke") R> set.seed( 12345) R> smoke.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, digits = 3, When a JAGS model is compiled, it may require an initial sampling phase during which the samplers adapt their behavior to maximize their efficiency (e.g., a Metropolis-Hastings random walk algorithm may change its step size) (Plummer 2016).The warning of "adaptation incomplete" may occasionally occur if the number of iterations for the adaptation process (i.e., the argument n.adapt) is not sufficient, so the MCMC algorithm may not achieve the maximum efficiency.This warning generally has little impact on the posterior estimates of the treatment effects.To avoid this warning, users may increase n.adapt.
The results are saved in the object smoke.out, a list containing AbsoluteRisk, OddsRatio, LogOddsRatio, RelativeRisk, LogRelativeRisk, RiskDifference, TrtRankProb, and DIC.We can use these effect size names to display the corresponding estimates.For example, the estimates of absolute risks (posterior mean and standard deviation, and posterior median and 95% credible interval) can be displayed as The argument digits in the function nma.ab.bin() can be used to change the number of digits to the right of the decimal point.Here, we used digits = 3.Each list element in the object smoke.outconsists of two sublists: Mean_SD contains posterior sample means with sample standard deviations; Median_CI contains posterior medians with 95% CIs.For example, users can output the medians and 95% CIs for the log odds ratio as follows: Since the log odds ratio is a relative effect size comparing a pair of treatments, the output estimates are displayed in a K × K matrix, where K is the number of treatments.In this example, K = 4.The element in the ith row and jth column is the estimated log odds ratio of treatment i compared to treatment j.To statistically test the difference between treatments, one can examine whether the effect size under the null hypothesis (i.e., the two treatments do not differ) is within the corresponding 95% CI; under the null hypothesis, ORs and RRs are 1, while LORs, LRRs, and RDs are 0. From the output, the 95% CI of the log odds ratio for SH vs. NC is (−0.032,1.610) which contains 0; therefore, there is not sufficient evidence to reject the null hypothesis.However, the 95% CI of the log odds ratio comparing GC to NC is (0.399, 1.970), which does not contain 0; therefore, these two treatments differ significantly.Also, in this example, we included "rank.prob" in the argument param to estimate treatment rank probabilities, and dic was specified as TRUE to calculate the DIC statistic.
Recall that for the smoking cessation dataset, a higher event probability implies a better treatment, so we specified the argument higher.betteras TRUE.Users can access the treatment rank probabilities and DIC statistics using From the output, treatment GC has the highest probability of being the best treatment (66.2%).As for the DIC statistic, D.bar is the posterior expectation of the deviance, which reflects the model fit; it is usually lower when more parameters are used in the model.However, complex models may lead to overfitting.To balance between the number of parameters and fitting effects, pD is used to penalize D.bar; it reflects the number of effective parameters used in the model.DIC is the penalized deviance, calculated as the sum of D.bar and pD; a model with smaller DIC is preferred.
Trace plots of LOR are generated because trace = "LOR" was specified.Figure 2 shows the trace plots of the LOR comparing IC and GC.Since we used the default n.chains = 3, three trace plots are drawn.Each trace plot shows evidence that the posterior samples of LOR are drawn from the stationary distribution.
A posterior density plot (Figure 3a) for treatment-specific absolute risks is also generated by specifying postdens = TRUE.This density plot is smoothed by the R function density().This plot shows visualized treatment effects, and we may also evaluate treatment differences.For example, NC clearly has lower event probability than IC and GC, and its posterior density only overlaps with the densities of IC and GC in tiny regions.
Function nma.ab.cont() for continuous outcomes For continuous outcomes, the arguments of nma.ab.cont() are mostly similar to those of nma.ab.bin().The major difference is that users need to specify the summaries of continuous outcomes (sample means and standard deviations) for the arguments mean and sd in nma.ab.cont().Also, the effect sizes to be estimated include continuous treatmentspecific effects and their differences.Users can specify the argument param as "mu" to estimate treatment-specific effects and "diff" to estimate effect differences.Also, "rank.prob"can be included in param to estimate treatment rank probabilities.The network dataset parkinson is used as an example: R> data("parkinson") R> set.seed(12345)R> parkinson.out <-nma.ab.cont(s.id, t.id, mean, sd, n, data = parkinson, + model = "hom_eqcor", prior.type = "unif", digits = 3, n.adapt = 10000, + n.iter = 100000, n.thin = 1, conv.diag = TRUE, trace = "mu", + postdens = TRUE) In this example, we used the model "hom_eqcor", which assumes an exchangeable correlation structure for the correlations between treatment effects.We display the medians and the corresponding 95% CIs for treatment-specific effects and effect differences as follows.
R> The effect difference in the ith row and jth column is calculated as the effect of treatment i minus that of treatment j.To statistically test the difference between two treatments, users may check whether 0 is within the corresponding 95% CI for the effect difference.For instance, from the output, the 95% CI for the effect difference between Treatments 1 and 2 is (−2.860,0.762), which does not contain 0. Therefore, Treatments 1 and 2 differ significantly.The posterior density plot (Figure 3b) for treatment-specific effects is obtained by specifying postdens = TRUE.From the density plot, the overlap region of densities for Treatments 1 and 2 is fairly small, and this supports the above conclusion.
Functions nma.ab.py() and nma.ab.followup() for count datasets For models (3) and ( 4), treatment effects can be related to the follow-up times of participants.Some studies report the total exposure times in person-years (e.g., Hooperet al. 2000) for each treatment group, and some report the mean follow-up time for each study (e.g., Psaty et al. 2003;Elliott and Meyer 2007).The functions nma.ab.py() and nma.ab.followup() can be used for these two types of datasets, corresponding to models (3) and (4), respectively.In these two functions, the argument param can include "rate" (treatment-specific rate), "lograte" (log rate), "ratio" (rate ratio), "logratio" (log rate ratio), and "rank.prob"(treatment rank probabilities).Since the two functions are similar, this article focuses on using the dataset diabetes to illustrate the usage and output of nma.ab.followup().From the output, all of the 95% CIs contain 0, so the difference between any pair of treatments is not significant.The posterior density plot for log rates is shown in Figure 3c, which also indicates that treatment-specific density curves do not differ much.

Plotting 95% credible intervals
When presenting network meta-analysis results, it is helpful to report the 95% CIs for effect sizes of interest.The package pcnetmeta provides functions absolute.plot()and contrast.plot()to draw the 95% CIs for treatment-specific and relative effect sizes, respectively.Users can simply call these functions on objects obtained from  The argument reference specifies the reference treatment to be compared against.The right panels of Figure 4 display the contrast plots.
Figure 4a shows the treatment-specific 95% CI plot for the smoking cessation data.Clearly, the 95% CIs of IC and GC do not overlap with the 95% CI of NC.This indicates significant differences between IC vs. NC and GC vs. NC.This can be confirmed by Figure 4b: The 95% CIs of ORs comparing IC vs. NC and GC vs. NC do not intersect with the vertical line at OR = 1. Figure 4c is the treatment-specific 95% CI plot for the Parkinson's disease data.Note that the 95% CIs for Treatments 1 and 2 overlap only in a tiny region.Correspondingly, the 95% CI for the effect difference between the two treatments in Figure 4d does not intersect with the vertical line at 0. Finally, Figures 4e and 4f show the 95% CI plots for treatment-specific and relative effects for the diabetes dataset.All of the treatment-specific 95% CIs have large overlap with other CIs.The 95% CIs of the rate ratios compared with placebo intersect with the vertical line at 1; therefore, the active treatments do not differ significantly from placebo.

Plotting treatment rank probabilities
Function rank.prob() is used to graph the probabilities of each treatment having each of the different possible ranks among the treatments.Users can call this function for the objects Figure 5 shows the plots of treatment rank probabilities.In the plots, each vertical bar represents probabilities that a specific treatment has different possible ranks.A darker area indicates the probability of being a higher rank, thus the black areas show the probabilities of being the best treatment.Therefore, from Figures 5a and 5b, treatments GC and Trt2 have much higher probabilities of being the best treatment, compared with other treatments in their respective studies.Figure 5c shows the rank probabilities plot for the dataset diabetes.Treatment ARB has highest probability of being the best treatment, although the probability for Treatment ACEI is close to the highest probability, so ARB and ACEI do not differ much.

Discussion
This article presents an overview of the R package pcnetmeta.Arm-based models are introduced to demonstrate the underlying methods of the functions.Practical usage of various functions is illustrated with examples of real network meta-analyses.Also, the package provides several plots for interpretation of network meta-analysis outputs.
MCMC convergence diagnostics have been extensively discussed in the literature (Cowles and Carlin 1996;Kass et al. 1998).The PSRFs and trace plots provided by the package pcnetmeta are used to examine whether the MCMC chains are drawn from stationary distributions; however, additional techniques are required to determine the effective sample size for adequate convergence (Robert and Casella 2004, p. 500).By specifying the argument mcmc.samples= TRUE in nma.ab.bin(), nma.ab.cont(), nma.ab.followup(), and nma.ab.py(), the MCMC posterior samples are saved in the output objects.Functions in other packages developed for MCMC convergence and samplesize adequacy, such as the R package mcmcse (Flegal and Hughes 2012), can be called for these posterior samples.
The current version of pcnetmeta does not detect inconsistency in arm-based network metaanalysis.Future work would add functions for network consistency assessment (Zhao et al. 2016).Moreover, both the contrast-based and arm-based methods can be extended to handle   Posterior density plots generated by the functions in package pcnetmeta.Plots of treatment rank probabilities generated by the function rank.prob().

Figure 2 .
Figure 2.Trace plots generated by R function nma.ab.bin() for the log odds ratio comparing IC and GC in the smoking cessation data.

Figure 4 .
Figure 4.The left three panels show the plots for treatment-specific (absolute) effects generated by the function absolute.plot();the right three panels show the plots for relative effects generated by the function contrast.plot().
Users need to input s.id and t.id for study and treatment IDs respectively.The argument title gives the graph title, and trtname specifies the treatment names.If trtname is not specified, the treatment IDs given in t.id are used.The argument alphabetic is a logical value indicating whether to sort the treatment nodes alphabetically in the clockwise direction according to the treatment names; if alphabetic = FALSE, the nodes are sorted according to treatment IDs ( t.id).The logical argument weight.edge=TRUEcauses the edge thickness to be drawn proportional to the number of direct treatment comparisons;weight.node=TRUE causes the node size to be proportional to the number of direct comparisons which contain that treatment node.The following code produces network plots for the datasets smoke, parkinson, and diabetes respectively.