 Software
 Open Access
 Published:
BUGSnet: an R package to facilitate the conduct and reporting of Bayesian network Metaanalyses
BMC Medical Research Methodology volume 19, Article number: 196 (2019)
Abstract
Background
Several reviews have noted shortcomings regarding the quality and reporting of network metaanalyses (NMAs). We suspect that this issue may be partially attributable to limitations in current NMA software which do not readily produce all of the output needed to satisfy current guidelines.
Results
To better facilitate the conduct and reporting of NMAs, we have created an R package called “BUGSnet” (Bayesian inference Using Gibbs Sampling to conduct a Network metaanalysis). This R package relies upon Just Another Gibbs Sampler (JAGS) to conduct Bayesian NMA using a generalized linear model. BUGSnet contains a suite of functions that can be used to describe the evidence network, estimate a model and assess the model fit and convergence, assess the presence of heterogeneity and inconsistency, and output the results in a variety of formats including league tables and surface under the cumulative rank curve (SUCRA) plots. We provide a demonstration of the functions contained within BUGSnet by recreating a Bayesian NMA found in the second technical support document composed by the National Institute for Health and Care Excellence Decision Support Unit (NICEDSU). We have also mapped these functions to checklist items within current reporting and best practice guidelines.
Conclusion
BUGSnet is a new R package that can be used to conduct a Bayesian NMA and produce all of the necessary output needed to satisfy current scientific and regulatory standards. We hope that this software will help to improve the conduct and reporting of NMAs.
Background
Indirect treatment comparisons (ITC) and network metaanalysis (NMA) are approaches for quantitatively summarizing an evidence base in which there are more than two treatments of interest. Unlike traditional pairwise metaanalysis, ITC/NMA can incorporate indirect evidence that arises when a group of studies evaluating different treatments share a common comparator. The incorporation of such evidence within an NMA has several advantages over pairwise metaanalysis [1, 2]. Unlike pairwise metaanalysis, an NMA allows for the comparison of two or more treatments that have never been directly compared provided that the studies examining such treatments are linked via a common comparator (i.e. an indirect comparison) [1, 2]. Another important advantage of NMA over pairwise metaanalysis is that it may provide greater statistical precision through its incorporation of indirect evidence which is not taken into account within pairwise metaanalysis [1, 2]. Lastly, an NMA can be used to rank a set of treatments for a given disease indication with respect to their clinically efficacy or harm and can be used to quantify the uncertainty surrounding such which is useful when determining policies, guidelines, and costs surrounding the choice of treatment [2].
The number of publications using NMA has increased dramatically within the past decade [3]. Despite this increase, several reviews have noted shortcomings with respect to the quality of the conduct and reporting of NMAs [4,5,6,7,8,9]. In particular, several authors have noted that a considerable proportion of NMAs do not provide a descriptive overview of the network or its structure, fail to adequately describe the statistical methods employed and whether or not their underlying assumptions were assessed and met, and lack a comprehensive summary of the results including effect estimates and measures of uncertainty regarding treatment ranks [4,5,6,7,8,9]. To improve the conduct, reporting, and appraisal of NMAs, a number of guidelines have been published which include the International Society of Pharmacoeconomics and Outcomes – Academy of Managed Care Pharmacy – National Pharmaceutical Council (ISPORAMCPNPC) questionnaire for assessing the relevance and credibility of an NMA [10], the Preferred Reporting Items for Systematic Reviews and MetaAnalyses (PRISMA) extension for reporting systematic reviews incorporating NMAs of health care interventions [11], and the National Institute for Health and Care Excellence Decision Support Unit (NICEDSU) reviewer’s checklist for appraising the synthesis of evidence within a submission to a health technology assessment agency (technical support document 7) [12].
Although the dissemination and uptake of such guidelines will hopefully help to address some of the foregoing issues, we suspect that the such issues may, in part, be related to limitations in current userfriendly software and tools used to conduct NMA. As previously noted, current software packages do not readily produce all of the output necessary to satisfy current reporting guidelines in a format that is suitable for submission to a journal or health technology assessment agency [13, 14]. Individuals must therefore rely upon multiple software packages, modify existing software, or generate code de novo in order to adhere to scientific and regulatory standards [14]. The resulting increase in time, effort, and expertise has likely impacted the quality and reporting of NMAs done to date. Furthermore, we have found that the documentation and help files of current software packages sometimes suffer from a lack of clarity regarding their implementation and use. In addition, the current lack of approachable tutorials that demonstrate how to use current NMA software could be a hindrance to users with limited programming expertise. To address these limitations, we have developed an R package called “BUGSnet” (Bayesian inference Using Gibbs Sampling to conduct a Network metaanalysis) aimed at improving the reporting and conduct of NMA/ITC. BUGSnet improves over its two main competing software packages for conducting a contrastbased Bayesian NMA: GeMTC [15] and NetMetaXL [16]. While NetMetaXL does produce much of the output necessary to satisfy reporting guidelines, it is limited in the types of analyses it can carry out. Specifically, one cannot use NetMetaXL to analyze outcomes that are not dichotomous, to conduct metaregression, or to analyzing evidence bases with more than 15 treatments [16]. While GeMTC provides an enhanced suite of functions for conducting NMA relative to NetMetaXL, its reporting capabilities are limited. For example, GeMTC does not readily produce key reporting items for an NMA such as tabular overview of the evidence base or a SUCRA plot and league table of the NMA results on the original scale.
Implementation
BUGSnet is a suite of functions that will carry out a Bayesian NMA while generating all items needed to satisfy the statistical components of the PRISMA, ISPORAMCPNPC, and NICEDSU checklists in a format that is suitable for publication or submission to a decisionmaking organization. These statistical components can be broadly categorized into: description of network (graphical and tabular), detection of heterogeneity, network metaanalysis (including metaregression), model assessment, detection of inconsistency and reporting of the results. An overview of BUGSnet’s functions and the corresponding checklist items that they address is presented in Table 1.
BUGSnet is implemented within R software. BUGSnet requires that the user have installed Just Another Gibbs Sampler (JAGS) on their computer [18, 19]. Information as to how to install JAGS can be found at the program’s sourceforge homepage: http://mcmcjags.sourceforge.net/. BUGSnet is hosted and can be accessed at the following URL: https://bugsnetsoftware.github.io/. We encourage users to submit feedback on existing code and to provide suggestions for additional functions that should be added to BUGSnet at the aforementioned homepage. Detailed vignettes describing the stepbystep use of BUGSnet to conduct an NMA on various types of outcomes are currently available in the R package documentation and on the BUGSnet homepage and additional applied examples are forthcoming.
Data preparation
The first step to using BUGSnet is to process the data using the data.prep() function where the user specifies the name of the columns variables that correspond to the study IDs and treatment arms. This way, the user does not have to enter this information over and over in subsequent functions.
Description of network
Current guidelines recommend that authors report plot of the evidence network [10,11,12]. The net.plot() and the net.tab() functions allow the user to describe the network of studies in a graphical and tabular format respectively.
With respect to the network graph, the size of the nodes and edges within the network plot are scaled such that they reflect the number of studies examining a specific treatment and the number of comparisons between any two given treatments respectively as per current recommendations. In addition, we have introduced an option that allows the user to highlight specific interventions of interest within the network graph and to label the edges with the names of the studies that have investigated these particular treatments. The colour, size, and layout of the network graph is highly customizable to ensure that the resulting figure will meet industry and journal standards.
The net.tab() function produces descriptive tables that are based on the tables produced by NetMetaXL – an excelbased software for conducting Bayesian NMAs [16]. While the tables produced by NetMetaXl are excellent descriptors of the network geometry, this software is currently only capable of handling dichotomous outcomes and is limited to 15 treatments [16]. We have expanded upon the tabular reporting of NetMetaXL by allowing such tables to summarize other types of outcomes including continuous, dichotomous, and count outcomes. An additional feature of our function is a report on whether the network is connected or not.
Homogeneity
Current guidelines recommend a careful exploration of heterogeneity within the network, typically prior to conducting the NMA [10,11,12]. Researchers should identify which characteristics are likely to be important modifiers of the treatment effects a priori using content expertise or a literature review [20]. Once identified, one can use the data.plot() function within BUGSnet to assess the heterogeneity of these modifiers within an evidence network. Specifically, this function generates a graph that allows the user to display a characteristic of interest within each treatment arm, grouped by study ID or treatment.
In addition, BUGSnet also provides an option within the pma() function to produce a table summarizing a Cochrane chisquare test, the tausquared statistic, and the Isquared statistic for assessing betweenstudy heterogeneity within each possible pairwise comparison within the network in which there is direct evidence [21].
Network metaanalysis
BUGSnet implements a Bayesian contrastbased NMA using a generalised linear model as described in the NICEDSU technical support document 2 [17]. The BUGS code used to generate these models within the BUGSnet package borrows heavily from this source [17]. Within BUGSnet, the nma.model() function is used to generate the BUGS model that one wishes to fit which includes aspects such as the link function and the likelihood distribution appropriate for the outcome of interest, the choice of using a fixed effects or a random effects model, and the inclusion of covariates if one wishes to conduct a metaregression. After the NMA model has been generated, one can run a Bayesian network metaanalysis with the function nma.run(). In the nma.run() function, the user can specify the number of burnins, iterations, and adaptations for the Markov Chain Monte Carlo (MCMC) algorithm and which variables they wish to monitor.
Bayesian inference
BUGSnet conducts NMA using Bayesian inference. There were several practical and theoretical reasons for choosing to implement the package within a Bayesian as opposed to a frequentist framework as noted by others: 1) Bayesian methods are more popular among researchers who conduct network metaanalyses; 2) Bayesian methods for network metaanalysis have been developed to a further degree; 3) Bayesian methods allow one to better handle data from trials with multiple arms and trials in which there are arms with zero events; 4) Bayesian methods are currently better suited for modeling uncertainty surrounding the heterogeneity between studies; 5) Bayesian methods present results as probabilities and are thus more suitable for ranking treatment efficacy and for incorporation into healtheconomic decision modeling [1, 22].
NMA models
BUGSnet can handle continuous, dichotomous, and count data (with or without varying followup times) as well as data from studies with more than two treatment arms. In what follows, we describe the NMA models that are implemented within BUGSnet. Suppose that we have data from studies i = 1, …, M. In arm k of study i, treatment t_{ik} ∈ {1, …, T} was used. The set {1, …, T} represents the set of treatments that were assessed across the M studies, where treatment 1 is a reference treatment. Let a_{1}, …, a_{M} represent the number of arms in studies 1, …, M. Let R_{ik} be the measured aggregate response in arm k of study i (e.g. proportion of individuals who were alive at oneyear, average blood pressure, etc.). Those responses are modeled as conditionally independent using an appropriate distribution F which is chosen based on the type of outcome at hand. For continuous outcomes, where the aggregate responses take the from of the sample mean and standard error in each arm, the distribution F is the normal distribution; \( {R}_{ik}\sim Normal\left({\varphi}_{ik},{se}_{ik}^2\ \right) \), where φ_{ik} is the mean and \( {se}_{ik}^2 \) is the observed standard error of the responses in arm k of study i. When outcome is dichotomous, the distribution F is the binomial distribution; R_{ik}~Binomial(n_{ik}, φ_{ik} ), where φ_{ik} is the probability of experiencing the event and n_{ik} is the sample size in arm k of study i. When outcomes take the form of counts and the event rates can be assumed to be constant over the duration of followup, one can use the Poisson distribution; R_{ik}~Poisson(e_{ik}φ_{ik} ), where e_{ik} is the observed persontime at risk and φ_{ik} is the event rate in arm k of study i. The latent parameters φ_{ik} ’s are transformed using an appropriate link function g(·) so g(φ_{ik}) ≡ θ_{ik} can be modeled with a linear model. Table 2 summarizes the link functions g(·) and family distributions F implemented within BUGSnet based on the type of outcome data. Following the NICEDSU technical support document 2 [17], the linear model used is generally of the contrastbased form:
where μ_{i} represents the fixed effect of the treatment from arm 1 in study i (a control treatment) and δ_{ik} represents the (fixed or random) effect of the treatment from arm k of study i relative to the treatment in arm 1 and δ_{i1} = 0 for i = 1, …,M. In BUGSnet, two exceptions to this model occur. First, when exploring a dichotomous outcome from studies with differing lengths of followup time, one can use a binomial family distribution with the complementary loglog link and the linear model includes the observed followup time f_{i} in trial i: θ_{ik} = log(f_{i}) + μ_{i} + δ_{ik} [17]. Second, when exploring a dichotomous outcome with a binomial family distribution and a log link, the linear model takes the form θ_{ik} = min(μ_{i} + δ_{ik}, −10^{−16}) to ensure that θ_{ik} is negative and the probabilities φ_{ik} are between 0 and 1.
In a random effect model, the \( {\boldsymbol{\delta}}_i'\mathrm{s}={\left({\delta}_{i2},\dots, {\delta}_{i{a}_i}\right)}^{\top } \) are modeled as conditionally independent with distributions
where \( {\mathbf{d}}_i={\left({d}_{\left({t}_{i1},{t}_{i2}\right)},\dots, {d}_{\left({t}_{i1},{t}_{i{a}_i}\right)}\right)}^{\top } \) and \( {d}_{\left({t}_{i1},{t}_{ik}\right)}={d}_{\left(1,{t}_{ik}\right)}{d}_{\left(1,{t}_{i1}\right)} \) is the difference in the treatment effect of treatments t_{i1} and t_{ik} on the g(·) scale and d_{(1, 1)} = 0. For Σ we adopt the usual compound symmetry structure described in (16), with variances σ^{2} and covariances 0.5σ^{2}, where σ^{2} represents the betweentrial variability in treatment effects (heterogeneity). Independent priors are used on σ, d_{(1, 2)}, …. ,d_{(1, T)} and μ_{1}, …, μ_{M}. For ease of implementation, in BUGSnet, the distribution (1) is decomposed into a series of conditional distributions [17].
In a fixed effect model, the δ_{ik} ’ s are treated as “fixed” (to use frequentist jargon) and are defined as \( {\delta}_{ik}={d}_{\left({t}_{i1},{t}_{ik}\right)}={d}_{\left(1,{t}_{ik}\right)}{d}_{\left(1,{t}_{i1}\right)} \) with d_{(1, 1)} = 0. Independent priors are used on d_{(1, 2)}, …. ,d_{(1, T)} and μ_{1}, …, μ_{M}. In both the fixed and randomeffects model, the posterior quantities of interest are all the mean treatment contrasts \( {d}_{\left({t}_{i1},{t}_{ik}\right)} \) which can be determined from d_{(1, 2)}, …. ,d_{(1, T)} through the transitivity relation \( {d}_{\left({t}_{i1},{t}_{ik}\right)}={d}_{\left(1,{t}_{ik}\right)}{d}_{\left(1,{t}_{i1}\right)}. \)
Metaregression
Let x_{ik} be a continuous covariate available in arms k = 1, …, a_{i} of studies i = 1, …, M. Network metaregression is implemented in BUGSnet via the linear model
where \( \overline{x} \) is the average of the x_{ik} ’s across studies and the \( {\beta}_{\left({t}_{i1},{t}_{ik}\right)}={\beta}_{\left(1,{t}_{ik}\right)}{\beta}_{\left(1,{t}_{i1}\right)} \) are regression coefficients for the effect of the covariate on the relative effect of treatments t_{i1} and t_{ik}, with β_{(1, 1)} = … = β_{(T, T)} = 0. A prior is used on β_{(1, 2)}, …, β_{(1, K)}. When conducting a metaregression analysis, the output plots and tables described in the Output section (league heat plot, league table, etc.) can also be produced but the user will need to specify a value for the covariate at which to produce treatment comparisons. Those treatment comparisons are calculated internally within BUGSnet by computing posterior quantities of interest at a specific covariate value x^{0} as \( {d}_{\left({t}_{i1},{t}_{ik}\right)}+{\beta}_{\left({t}_{i1},{t}_{ik}\right)}\left({x}^0\overline{x}\right), \) and using the transitivity relations \( {d}_{\left({t}_{i1},{t}_{ik}\right)}={d}_{\left(1,{t}_{ik}\right)}{d}_{\left(1,{t}_{i1}\right)} \) and \( {\beta}_{\left({t}_{i1},{t}_{ik}\right)}={\beta}_{\left(1,{t}_{ik}\right)}{\beta}_{\left(1,{t}_{i1}\right)}. \)
Choice of priors
By default, BUGSnet implements the vague priors described in Table 3. Our choice of priors was based on the justifications made by van Valkenhoef et al. (2012) [15] which allow a prior variance to be easily calculated from the data without any user input. These priors are the same as the ones implemented in the GeMTC R package [15]. The user also has the option within the nma.model() function to specify their own prior which is useful for conducting sensitivity analyses, namely for the comparison of prior distributions on the random effects standard deviation, σ, to insure that they do not have a significant effect on the posterior estimates.
The variances 15u are taken from van Valkenhoef (2012) et al., where u is the largest maximum likelihood estimator of treatment differences on the linear scale in single trials [15]. Note that t denotes Student’s t distribution with parameters: location, variance and degrees of freedom.
Model assessment
After the NMA model has been run, guidelines recommend that one assesses the convergence and fit of the model [10,11,12]. In BUGSnet, convergence can be assessed using trace plots and other convergence diagnostics produced by the nma.diag() function. Lastly, the fit of the model and the identification of potential outliers can be carried out using the nma.fit() function which will produce a plot of the leverage values and also display the corresponding effective number of parameters, total residual deviance, and deviance information criterion (DIC). These latter values can be used to help determine or justify model choice when considering two or more competing models (e.g. between a fixed or randomeffects model) and to help identify data points that contribute heavily to the DIC and/or that are influential.
Consistency
A fundamental assumption of an NMA is the assumption of transitivity [2]. Under this assumption, one assumes that one can estimate the difference in the effect of two treatments by subtracting the difference in the effects of the two treatments relative to a common comparator as follows: \( {d}_{\left({t}_{i1},{t}_{ik}\right)}={d}_{\left(1,{t}_{ik}\right)}{d}_{\left(1,{t}_{i1}\right)} \)[2]. Aside from exploring clinical heterogeneity of treatment definitions and modifiers within the network using the data.plot() function, one can also detect violations of the assumption of transitivity by examining statistical consistency within the network. Statistical consistency refers to the statistical agreement between indirect and direct evidence within an evidence network [2]. Evidence of inconsistency would indicate a violation of the transitivity assumption. As noted by Efthimiou et al. (2015), statistical consistency can only be explored if there are closed loops within the network [2]. A variety of methods have been proposed to assess consistency within a network metaanalysis [2, 24, 25]. Such methods are often categorized as being “global” or “local” depending upon whether they examine inconsistency within the entire network or within particular segments thereof [2]. BUGSnet currently implements the inconsistency model (or unrelated mean effects model) as described in the NICEDSU TSD 4 [26]. An inconsistency model is an NMA model similar to the consistency models described above but transitivity \( {d}_{\left({t}_{i1},{t}_{ik}\right)}={d}_{\left(1,{t}_{ik}\right)}{d}_{\left(1,{t}_{i1}\right)} \) is not assumed. Instead, independent priors are defined on each of the \( {d}_{\left({t}_{i1},{t}_{ik}\right)} \) ’s. Inconsistency models therefore have more parameters than consistency models, which needs to be weighted against how well they fit the data compared to the consistency model to determine if there is evidence of inconsistency. The inconsistency model can be specified using the type = "inconsistency" option in the nma.model(). To examine inconsistency at the global level, the fit of the inconsistency model can be compared against a model in which consistency is assumed using the nma.fit() function and comparing the DICs. Local inconsistency can be explored on the leverage plots produced by nma.fit() and also using the nma.compare() function which produces a plot comparing the posterior mean deviance of each data point between the consistency and the inconsistency model.
We chose to implement the inconsistency model method for assessing inconsistency in BUGSnet because it easily handles different network structures and multiarm trials, which is not the case with other methods for assessing inconsistency such as the Bucher method [26, 27]. More options for assessing inconsistency at both the global and local levels will be considered in further BUGSnet releases.
Output
We provide several functions for displaying the results of the NMA in both graphical and tabular formats (league tables, league heat plots, SUCRA plots, SUCRA tables, rankograms and forest plots) to satisfy current guidelines. With respect to plotting the magnitude and uncertainty of the treatment effects, users can the use the nma.forest() function to graph the effect estimates from the NMA against a comparator specified by the user. The effect estimates can also be presented within a league table using the nma.league() function. An important presentation feature in BUGSnet, particularly for large league tables, is that the user can specify an option to colour and arrange the league table into a heatmap that highlights the magnitude of the effect estimates. Users can also graphically display the probability of the ranking of each treatment within a surface under the cumulative ranking curve (SUCRA) plot which can be specified within the nma.rank() function. This function can also be used to present treatment ranks in a tabular format, extract SUCRA values and produce a rankogram. All of the plots produced by these three reporting functions are produced with the ggplot2 package. As such, the user can easily customize the plots (e.g. change the background, add a title) by adding layers using the + command. Also, for reporting relative treatment effects the user can specify whether they want to plot the results on the the linear scale (log scale) or the original scale.
When metaregression is conducted, the nma.rank(), nma.forest() and nma.league() functions allow the user to specify for which value of the covariate they wish to present the results. Even though the covariate is centered for metaregression, the user does not have to do any conversion and results are provided on the original noncentered scale. Another function, nma.regplot() outputs a plot of relative treatment effects on the linear scale across the range of covariate values used in the metaregression, as in the NICEDSU TSD 3 [28].
It is sometimes recommended that users present results from the direct evidence where available [29]. To accommodate this, we have also incorporated the pma() function within BUGSnet which will perform pairwise metaanalysis using the meta package in R and automatically output the results into a tabular format [30].
Results
The following is a demonstration of some of the functions contained within BUGSnet (Table 1) and some of the possible outputs. To accomplish this task, we have recreated an analysis of a dichotomous outcome where studies had variable followup times described in the NICEDSU technical support document 2 (referred to as “Data Example 3”) [17]. The BUGSnet code used to produce this analysis is available in the vignette titled survival in the BUGSnet documentation, and appended as a supplement to this article (see Additional file 1). Additional outputs are presented in the vignette as well as a more detailed description of how to conduct and report network metaanalysis, which is only presented here in brief.
The evidence network used in this analysis consists of 22 randomized trials (including multiarm trials) that examined the effects of six antihypertensive treatments on the risk of developing diabetes [31]. The outcome for this data is the number of new diabetes cases observed during the trial period. The data is organized in the long format (i.e. one row per treatment arm), with variables indicating the study ID, the treatment ID, the number of patients, the number of events, and the mean age (and standard deviation) of participants for each treatment arm (see Table 4). The results of our package are concordant with those reported in the TSD as well results obtained with GeMTC (code and outputs provided as supplement to this article (see Additional files 2, 3, 4 and 5) and NetMetaXL.
Data preparation, description of network and homogeneity
After the data was prepared using the data.prep() function, the net.plot() and the net.tab() functions were used to describe the network of studies in a graphical (Fig. 1) and tabular format respectively (Table 5). As previously discussed, the assumptions of network metaanalysis will be violated when an effect modifier is heterogeneously distributed throughout an evidence base [20]. Prior to conducting the network metaanalysis, analysts can use the data.plot() function to examine the distribution of an effect modifier within the network. The determination of whether or not a variable is an effect modifier and whether or not the observed differences in its distribution are clinically meaningful is determined according to expert opinion and prior evidence. To demonstrate this function, we have simulated a patient characteristic that may modify the treatment effect (i.e. the age of participants). To mimic a lack of reporting, we have omitted the standard deviation for a few of the studies. As observed in Fig. 2, the mean age of participants within each treatment arm (the individual points) is similar to the overall mean age of participants within the evidence base (the red dotted line). According to the standard deviation (the +/− error bars), the variability of ages within each treatment arm appear to be similar as well (where available). Based on this analysis, one would conclude that there is no meaningful heterogeneity in the distribution of age. This analysis would be repeated for all potentially important effect modifiers identified a priori by clinical opinion and a review of previous studies. If no heterogeneity is detected, then one may proceed to conducting the network metaanalysis. If heterogeneity is detected, one can attempt to adjust for imbalances by using metaregression (if there are an adequate number of studies) or by using alternative statistical techniques that leverage individual patient data (e.g. matchingadjusted indirect comparison or simulated treatment comparison) [20].
Network metaanalysis
We conducted an NMA on the Diabetes dataset by fitting a generalized linear model with a complementary loglog link function and binomial likelihood function to account for the dichotomous outcome and differing followup times between studies, which was specified through the use of nma.model(). To be consistent with the NICEDSU technical support document, we specified a burnin of 50,000 iterations followed by 100,000 iterations with 10,000 adaptations in the nma.run() function. We compared the fit of both a fixed and randomeffects model. According to a visual examination of the leverage plots and comparison of the DIC values produced by the nma.fit(), the random effects model would be preferred over the fixed effects model for this particular dataset because the DIC value is lower and because there are fewer outliers in the leverage plot (Fig. 3).
Output
We present results from the generalized linear model that we previously fit to the Diabetes dataset. As visualized in the SUCRA plot obtained from nma.rank(), the angiotensinreceptor blockers’ (ARB) curve is consistently above the other treatments’ curves suggesting that it is the most beneficial treatment with respect to the outcome among the treatments included in the Diabetes evidence network (Fig. 4). The effect estimates and credible intervals produced by the foregoing model are displayed in a league heat plot (Fig. 5) obtained using nma.league(). In Fig. 5, one can see that the difference between ARB and other treatments are all statistically significant at the 95% level except for the ACE inhibitor and Placebo treatments.
Consistency
To assess the presence of inconsistency, we fit an NMA model similar to the one previously described but assuming inconsistency. We obtain leverage plots similar to Fig. 3 using the nma.fit() function where we find that the DIC for the consistency model is marginally smaller than for the inconsistency mode. We also use the nma.compare() function to plot the individual data points’ posterior mean deviance contributions for the consistency model vs the inconsistency model (Fig. 6) as recommended in the NICEDSU TSD 4 [26]. Overall, we conclude that there is a lack of evidence to suggest inconsistency within the network.
Discussion
BUGSnet is intended to be used by researchers when assessing the clinical efficacy of multiple treatments within the context of a submission to a journal or a health technology assessment agency. For conducting a contrastbased Bayesian NMA, the two main competing software packages that one may consider are GeMTC [15] and NetMetaXL [16], for which we have discussed limitations in the introduction. With BUGSnet, we aimed to create a single tool that would compete with the reporting capabilities of NetMetaXL and the analytic capabilities of GeMTC. We have also aimed to provide users with enhanced reporting options not included in existing software such as a function to produce graphs that show the distribution of effect modifiers by trial or by treatment arm and an option to print study names and highlight certain treatment comparisons within the network plot. To help facilitate the use of BUGSnet among new users, we have provided three vignettes (with more vignettes forthcoming) in the R help files that walk users through conducting an NMA using BUGSnet by providing detailed R code and interpretations of the statistical output. Despite these benefits, there are limitations of BUGSnet. BUGSnet is currently limited to exclusively analyzing armlevel data. In contrast, GeMTC can be used to conduct an NMA using entirely armlevel or entirely contrastlevel data [22]. Relative to GeMTC, another limitation of BUGSnet is that GeMTC currently provides a broader range of methods of assessing inconsistency such as the nodesplitting method and a broader range of metaregression analyses such as subgroup metaanalysis. Since it is implemented within the R environment, some users may find BUGSnet more difficult to use relative to NetMetaXL, which is implemented within Microsoft Excel. At this point, armbased models [22] have not been implemented in BUGSnet; the R package pcnetmeta allows such analyses, although it does not readily provide a complete suite of outputs like BUGSnet. We plan to address these shortcomings in future iterations of BUGSnet and interested users should check the previously mentioned URL for updates.
Network metaanalysis is a rapidly evolving area of research with new methods constantly being developed [32]. While the work presented within this paper provides the essential tools required to conduct an NMA in accordance with current guidelines, we plan to implement additional functions and features within this package, based on user feedback, to provide enhanced flexibility and to ensure relevance. Some of the preliminary requests for shortterm additions include: 1) additional functions for detecting inconsistency within the network such as the Bucher method [27]; 2) an option to allow the user to conduct an NMA using studylevel effect estimates; 3) allowing for the relaxation of the proportional hazards assumption when analyzing timetoevent outcomes; 4) allowing for subgroup metaregression and the inclusion of more than one covariate into the metaregression model; 5) a function that will automatically generate a report or slide deck presentation of the results that could be saved as a pdf, html or Word.
As detailed in Table 1, the functions contained within BUGSnet can be used to address the items within the PRISMA, ISPORAMCPNPC, and NICEDSU reporting guidelines that are related to the statistical analysis component of an NMA [11, 12, 29]. However, it should be emphasised that there are several nonstatistical issues described within these guidelines that BUGSnet is not meant to address such as the identification of the research question, the specification of the study population and competing interventions, the development of the search strategy, and the assessment of the risk of bias within each study [10,11,12]. Researchers are urged to consult with these guidelines when planning their NMA to ensure that all aspects of the NMA, both statistical and nonstatistical, adhere to current reporting and methodologic standards.
Conclusions
Here we present a new JAGSbased R package for conducting Bayesian NMA called BUGSnet. Relative to existing NMA software, BUGSnet provides an enhanced set of tools for conducting and reporting results according to published bestpractice guidelines to help overcome the lack of quality identified within this body of literature. In addition to these features, we have attempted to provide ample documentation describing the use and implementation of BUGSnet to help promote the understanding and uptake of this software. Lastly, we plan to monitor the literature and to implement new features within BUGSnet based on the NMA analyst community to ensure that the package remains uptodate with the latest advances in this rapidly developing area of research.
Availability and requirements
Project name: BUGSnet
Project home page: https://bugsnetsoftware.github.io/
Operating system(s): Windows 10 v1809 and Mac OS 10.14 (may work on earlier versions but not tested)
Programming language: R
Other requirements: JAGS 4.3.0
License: Creative Commons AttributionNonCommercialShareAlike 4.0 International
Any restriction to use by nonacademics: Contact authors for nonacademic use.
Availability of data and materials
All of the datasets and material contained within the manuscript can be accessed within the BUGSnet package via the BUGSnet homepage: https://bugsnetsoftware.github.io/
Abbreviations
 ISPORAMCPNPA:

International Society for Pharmacoeconomics and Outcomes Research  Academy of Managed Care Pharmacy  National Pharmaceutical Council
 ITC:

Indirect Treatment Comparisons
 JAGS:

Just Another Gibbs Sampler
 NICEDSU:

National Institute for Health and Care Excellence Decision Support Unit
 NMA:

Network MetaAnalysis
 PRISMA:

Preferred Reporting Items for Systematic Reviews and MetaAnalysis
 SUCRA:

Surface Under the Cumulative Ranking Curve
References
Jansen JP, Fleurence R, Devine B, Itzler R, Barrett A, Hawkins N, Lee K, Boersma C, Annemans L, Cappelleri JC. Interpreting indirect treatment comparisons and network metaanalysis for healthcare decision making: report of the ISPOR task force on indirect treatment comparisons good research practices: part 1. Value Health. 2011;14(4):417–28.
Efthimiou O, Debray TP, van Valkenhoef G, Trelle S, Panayidou K, Moons KG, Reitsma JB, Shang A, Salanti G. GetReal in network metaanalysis: a review of the methodology. Res Synth Methods. 2016;7(3):236–63.
Lee AW. Review of mixed treatment comparisons in published systematic reviews shows marked increase since 2009. J Clin Epidemiol. 2014;67(2):138–43.
Bafeta A, Trinquart L, Seror R, Ravaud P. Reporting of results from network metaanalyses: methodological systematic review. BMJ. 2014;348:g1741.
Hutton B, Salanti G, Chaimani A, Caldwell DM, Schmid C, Thorlund K, Mills E, CatalaLopez F, Turner L, Altman DG, et al. The quality of reporting methods and results in network metaanalyses: an overview of reviews and suggestions for improvement. PLoS One. 2014;9(3):e92508.
Song F, Loke YK, Walsh T, Glenny AM, Eastwood AJ, Altman DG. Methodological problems in the use of indirect comparisons for evaluating healthcare interventions: survey of published systematic reviews. BMJ. 2009;338:b1147.
Donegan S, Williamson P, Gamble C, TudurSmith C. Indirect comparisons: a review of reporting and methodological quality. PLoS One. 2010;5(11):e11054.
Kovic B, Zoratti MJ, Michalopoulos S, Silvestre C, Thorlund K, Thabane L. Deficiencies in addressing effect modification in network metaanalyses: a metaepidemiological survey. J Clin Epidemiol. 2017;88:47–56.
Zarin W, Veroniki AA, Nincic V, Vafaei A, Reynen E, Motiwala SS, Antony J, Sullivan SM, Rios P, Daly C, et al. Characteristics and knowledge synthesis approach for 456 network metaanalyses: a scoping review. BMC Med. 2017;15(1):3.
Jansen JP, Trikalinos T, Cappelleri JC, Daw J, Andes S, Eldessouki R, Salanti G. Indirect treatment comparison/network metaanalysis study questionnaire to assess relevance and credibility to inform health care decision making: an ISPORAMCPNPC good practice task force report. Value Health. 2014;17(2):157–73.
Hutton B, Salanti G, Caldwell DM, Chaimani A, Schmid CH, Cameron C, Ioannidis JP, Straus S, Thorlund K, Jansen JP, et al. The PRISMA extension statement for reporting of systematic reviews incorporating network metaanalyses of health care interventions: checklist and explanations. Ann Intern Med. 2015;162(11):777–84.
Ades A, Caldwell DM, Reken S, Welton NJ, Sutton AJ, Dias S. NICE DSU technical support document 7: evidence synthesis of treatment efficacy in decision making: a reviewer’s checklist. Decision Support Unit London UK. 2012;27905719.
Neupane B, Richer D, Bonner AJ, Kibret T, Beyene J. Network metaanalysis using R: a review of currently available automated packages. PLoS One. 2014;9(12):e115065.
Xu C, Niu Y, Wu J, Gu H, Zhang C. Software and package applicating for network metaanalysis: a usagebased comparative study. J Evid Based Med. 2018;11(3):176–83.
van Valkenhoef G, Lu G, de Brock B, Hillege H, Ades AE, Welton NJ. Automating network metaanalysis. Res Synth Methods. 2012;3(4):285–99.
Brown S, Hutton B, Clifford T, Coyle D, Grima D, Wells G, Cameron C. A Microsoftexcelbased tool for running and critically appraising network metaanalysesan overview and application of NetMetaXL. Systematic Reviews. 2014;3:110.
Dias S, Welton NJ, Sutton AJ, Ades A. NICE DSU technical support document 2: a generalised linear modelling framework for pairwise and network metaanalysis of randomised controlled trials; 2011.
Plummer M: JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. In: Proceedings of the 3rd international workshop on distributed statistical computing: 2003: Vienna; 2003.
Plummer M: JAGS version 4.3.0 user manual. https://sourceforge.net/projects/mcmcjags/files/Manuals/. Accessed 30 Aug 2019.
Jansen JP, Naci H. Is network metaanalysis as valid as standard pairwise metaanalysis? It all depends on the distribution of effect modifiers. BMC Med. 2013;11:159.
Higgins JP, Thompson SG. Quantifying heterogeneity in a metaanalysis. Stat Med. 2002;21(11):1539–58.
Dias S, Ades A, Welton NJ, Jansen JP, Sutton AJ. Network metaanalysis for decisionmaking. Hoboken: Wiley; 2018.
Warn DE, Thompson SG, Spiegelhalter DJ. Bayesian random effects metaanalysis of trials with binary outcomes: methods for the absolute risk difference and relative risk scales. Stat Med. 2002;21(11):1601–23.
Veroniki AA, Vasiliadis HS, Higgins JP, Salanti G. Evaluation of inconsistency in networks of interventions. Int J Epidemiol. 2013;42(1):332–45.
Donegan S, Williamson P, D'Alessandro U, Tudur Smith C. Assessing key assumptions of network metaanalysis: a review of methods. Res Synth Methods. 2013;4(4):291–323.
Dias S, Welton NJ, Sutton AJ, Caldwell DM, Lu G, Ades A. NICE DSU technical support document 4: inconsistency in networks of evidence based on randomised controlled trials; 2011.
Bucher HC, Guyatt GH, Griffith LE, Walter SD. The results of direct and indirect treatment comparisons in metaanalysis of randomized controlled trials. J Clin Epidemiol. 1997;50(6):683–91.
Dias S, Sutton AJ, Welton NJ, Ades A. NICE DSU technical support document 3: heterogeneity: subgroups, metaregression, bias and biasadjustment; 2011.
Jaime Caro J, Eddy DM, Kan H, Kaltz C, Patel B, Eldessouki R, Briggs AH. Questionnaire to assess relevance and credibility of modeling studies for informing health care decision making: an ISPORAMCPNPC good practice task force report. Value Health. 2014;17(2):174–82.
Schwarzer G. Meta: an R package for metaanalysis. R News. 2007;7(3):40–5.
Elliott WJ, Meyer PM. Incident diabetes in clinical trials of antihypertensive drugs: a network metaanalysis. Lancet. 2007;369(9557):201–7.
Nikolakopoulou A, Mavridis D, Furukawa TA, Cipriani A, Tricco AC, Straus SE, Siontis GCM, Egger M, Salanti G. Living network metaanalysis compared with pairwise metaanalysis in comparative effectiveness research: empirical study. BMJ. 2018;360:k585.
Acknowledgements
We would like to thank the two reviewers for their feedback which helped improve this manuscript. We would also like to thank the attendees of our 2019 CADTH Symposium workshop and UBC Network MetaAnalysis Workshop for providing valuable feedback regarding the implementation and ease of use of BUGSnet. AB acknowledges the support of a startup grant from the University of Waterloo and DBo acknowledges the support of the Barrie I Strafford Doctoral Scholarship for Interdisciplinary Studies on Aging.
Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.
Author information
Authors and Affiliations
Contributions
All five authors (AB, DBo, PA, JS and DBr) contributed to the conception and design of the work. AB designed the software architecture. AB, DBo and JS coded the software. DBo and AB drafted the first version of the manuscript and all five authors (AB, DBo, PA, JS and DBr) substantially revised it. All five authors (AB, DBo, PA, JS and DBr) have read and approved the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable
Competing interests
AB and DBo are consultants for Lighthouse Outcomes. JS is employed by Lighthouse Outcomes. PA and DBr are partners at Lighthouse Outcomes.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Additional file 1.
Vignette on network metaanalysis of survival data. A vignette detailing how to obtain the outputs in the Results section using BUGSnet version 1.0.2. The vignette includes all the necessary R code as well as additional outputs and explanations that were not presented in this manuscript for the sake of brevity.
Additional file 2.
R code to compare BUGSnet with GeMTC for the analysis of the diabetes data.
Additional file 3.
Diabetes dataset to accompany the R code in Additional file 2.
Additional file 4.
BUGSnet league table obtained from the R code in Additional file 2.
Additional file 5.
GeMTC league table obtained from the R code in Additional file 2.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Béliveau, A., Boyne, D.J., Slater, J. et al. BUGSnet: an R package to facilitate the conduct and reporting of Bayesian network Metaanalyses. BMC Med Res Methodol 19, 196 (2019). https://doi.org/10.1186/s1287401908292
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1287401908292
Keywords
 Network metaanalysis
 Indirect treatment comparison
 Systematic review
 Bayesian inference
 Knowledge synthesis
 Health technology assessment
 Clinical efficacy
 R package
 Reporting guidelines