 Software
 Open access
 Published:
crossnma: An R package to synthesize crossdesign evidence and crossformat data using network metaanalysis and network metaregression
BMC Medical Research Methodology volume 24, Article number: 169 (2024)
Abstract
Background
Although aggregate data (AD) from randomised clinical trials (RCTs) are used in the majority of network metaanalyses (NMAs), other study designs (e.g., cohort studies and other nonrandomised studies, NRS) can be informative about relative treatment effects. The individual participant data (IPD) of the study, when available, are preferred to AD for adjusting for important participant characteristics and to better handle heterogeneity and inconsistency in the network.
Results
We developed the R package crossnma to perform crossformat (IPD and AD) and crossdesign (RCT and NRS) NMA and network metaregression (NMR). The models are implemented as Bayesian threelevel hierarchical models using Just Another Gibbs Sampler (JAGS) software within the R environment. The R package crossnma includes functions to automatically create the JAGS model, reformat the data (based on user input), assess convergence and summarize the results. We demonstrate the workflow within crossnma by using a network of six trials comparing four treatments.
Conclusions
The R package crossnma enables the user to perform NMA and NMR with different data types in a Bayesian framework and facilitates the inclusion of all types of evidence recognising differences in risk of bias.
Background
Background to network metaanalysis
Studies that estimate treatment effects are identified, evaluated, and synthesized in systematic reviews to obtain evidence that answers treatmentrelated questions [1]. Systematic reviews may include a pairwise metaanalysis (PMA) [2, 3], which is a statistical summary of findings from multiple studies comparing two interventions. The PMA extends to network metaanalysis (NMA) to compare multiple competing interventions, providing estimates for the relative effects of each pair of competing treatments [4]. The final NMA estimates are a combination of direct estimates derived from combining study findings and indirect estimates obtained using one or more intermediate comparators under the consistency assumption.
Most NMAs use aggregate data (AD) obtained from published studies. To explore betweenstudy heterogeneity and betweencomparison inconsistency in NMA, we need to study the role of important patient and study characteristics, typically in subgroup analyses or metaregressions [5]. As relationships at study level often fail to reflect associations at individual patient level, aggregated data are not suitable to explore the role of patientlevel characteristics in modifying the treatment effects [6]. However, retrieving individual participant data (IPD) is a difficult and timeconsuming endeavour. The most common scenario is to obtain IPD of some of the included studies, and then combine IPD and AD in a single model.
The vast majority of published NMAs synthesise data from randomised clinical trials (RCTs). Although RCTs are by design less prone to selection bias than nonrandomised designs, biases can still arise from in their conduct [7,8,9] and reported findings [10] or from conflict of interest that can distort the body of evidence [11]. Generalisability of their findings to inform clinical practice is challenged by several of their features: RCTs include patients not necessarily representative of those encountered at the point of care, they are more likely to use placebo or other legacy treatments which are not an option in practice [10,11,12], and they often do not provide data on longterm benefit and safety of interventions. Pragmatic trials offer an alternative, but they can also be costly and difficult to conduct to study rare conditions, so that randomized evidence can be sparse in some health research fields. For these reasons, researchers on evidence synthesis consider sometimes evidence from nonrandomised studies (NRS) despite the risk of confounding and several other biases inherent in their design. To assess the risk of bias (RoB) in the results of a study, tools have been created: the ROBINSI for NRS [13] and the RoB 2 for RCT [14].
To handle different types of studies and data types, we recently introduced a suite of four Bayesian NMA and network metaregression (NMR) models [15] extended a previously described threelevel hierarchical model that combines IPD and AD [16,17,18] and included methods to combine these data when they come from RCT and NRS. The four models can be broadly described as an unadjusted (“naïve”) NMA model, where differences in bias between different study designs are ignored; a model that synthesises only RCT data with prior distributions for the relative treatment effects constructed from NRS evidence and potentially discounted according to RoB; and two biasadjusted models where the relative treatment effect of each study is adjusted according to the underlying risk of bias.
We have implemented the different crossNMA models in a new R package called crossnma (crossdesign and crossformat network metaanalysis and network metaregression). The package enables researchers to perform NMA and NMR on data that are available in different formats as IPD, AD or a combination of both, and each format can come from different study designs as RCT, NRS, or a mixture of both.
This work has been done within the HTx Horizon 2020 project. HTx is supported by the European Union, lasting for 5 years from January 2019. The main aim of HTx is to create a framework for the Next Generation Health Technology Assessment (HTA) to support patientcentered, societally oriented, realtime decisionmaking on access to and reimbursement for health technologies throughout Europe.
Synthesis models available in crossnma
Below we provide a brief description of the four NMR models implemented in crossnma; NMA models can be obtained by simply ignoring covariate terms. Of note, the NMR model is described with one covariate for simplicity and ease of explanation. However, our crossnma package is designed to handle up to three covariates. The model used in the analysis is defined in R function crossnma.model(). The notation used is summarized in Table 1. More details are available in the methodological publication [15].
Unadjusted synthesis of data from RCTs and NRS
This analysis combines the RCT and NRS evidence without accounting for the fact that different studies pertain to different risk of bias (argument method.bias = "naive" in crossnma.model()). Data from each study can be available as IPD or AD. The different types of studies and data format are combined using a threelevel hierarchical model.
Level 1 (individual participant level)
When IPD is available, we observe the outcome \({y}_{ijk}\) for participant \(i\) in study \(j\) receiving treatment \(k\). Each study has a reference treatment \(b\). Observed outcomes are assigned an appropriate likelihood distribution with unknown parameters \({\varphi }_{ijk}\). Then, using a link function \(g(.)\), these parameters are linked to the treatment and covariate effects. The distribution and link function in the model depend on the nature of the data and the effect measure we want to estimate. For example, when we observe binary data and want to estimate odds ratios, a Bernoulli likelihood is considered for the outcome and \(g(.)\) is the logit function. The general form of the NMR model is
where \({x}_{ijk}\) and \({\overline{x} }_{j}\) are covariate value and its study mean, respectively. The prognostic effect of the covariate is quantified by \({\beta }_{0j}\); \({\beta }_{1bk}^{B}\) is the betweenstudy interaction effect, which represents the associations between treatment and study’s mean covariate, and \({\beta }_{1bk}^{W}\) quantifies the withinstudy interactions between treatment effect and covariate. When \({x}_{ijk}={\overline{x} }_{j}=0,\) \({u}_{jb}\) can be interpreted as the average \(g(.)\)transformed outcome in the study reference arm \(b\), and \({\delta }_{jbk}\) is a studyspecific relative effect of treatment \(k\) versus \(b\).
Level 2 (studylevel)
In AD studies, we observe the mean outcome per study arm, \({y}_{jk}\), which is also assigned an appropriate likelihood distribution with unknown parameter \({\varphi }_{.jk}\). Then this parameter is linked to the model parameters via the link function:
Level 3 (crossstudies synthesis)
We combine the parameters from different studies assuming either a randomeffects or a commoneffect model. Table 2 summarizes the different assumptions supported by crossnma along with corresponding arguments to set each assumption in the package. It is important to highlight that studies with AD do not contribute to the estimation of \({\beta }_{0j}\) and \({\beta }_{1bk}^{W}\).
By default, we assign minimally informative prior distributions to the main parameters following Valkenhoef et al. [19]: normal distribution \(\mathrm{\rm N}(0,{\left\{15*{ML}_{max}\right\}}^{2} )\) for \({u}_{jb},{\beta }_{0j}, {B}_{0}, {B}_{1}^{W},{B}_{1}^{B},{d}_{Ak}\) and uniform distribution \({\text{Unif}}(0, {ML}_{max})\) for \(\tau ,{\tau }_{{B}_{0}},{\tau }_{B},{\tau }_{W}\). The quantity \({ML}_{max}\) is calculated by first computing the maximum likelihood estimate of the relative treatment effect for each study and then taking the maximum value of these estimates. Users can provide their own prior distribution for the heterogeneity standard deviation parameters \(\tau ,{\tau }_{{B}_{0}},{\tau }_{B},{\tau }_{W}\).
Use NRS priors for basic parameters in RCT model
In this model, we use NMR to estimate the relative treatment effects using only the NRS (argument method.bias = "prior" in crossnma.model(); arguments starting with run.nrs can be used to control this process). We use the mean \({d}_{Ak}^{NRS}\) and variance \({V}_{AK}^{NRS}\) of the NRS posterior distribution of the relative treatment effects to construct prior distributions for the treatment effects and fit the model in the RCT data; \({{d}_{Ak}\sim \mathrm{\rm N}(d}_{Ak}^{NRS},{V}_{AK}^{NRS})\). To limit the impact of NRS on RCT estimates, we can either inflate the prior variance by dividing it by a common inflation factor \(w\) with \(0<w<1\) or shift NRS means by \(\zeta\) (see reference [20] for more discussion on how to set \(\zeta\) or reference [21] to choose a value based on elicited expert opinion). Treatment effects not observed in NRS are given the default priors; \({d}_{Ak}\sim \mathrm{\rm N}(0,{\left\{15*{ML}_{max}\right\}}^{2} )\).
Biasadjusted model 1
In this and the next model, we adjust treatment effects according to each study’s RoB [22]. The level 1 model in Unadjusted synthesis of data from RCTs and NRS section is extended for biasadjusted model 1 (method.bias = "adjust1") is extended as follows (additional terms are printed in bold):
and the level 2 model becomes \(g\left({\varphi }_{.jk}\right)={u}_{jb} + {{\varvec{\delta}}}_{{\varvec{j}}{\varvec{b}}{\varvec{k}}}{{\varvec{\gamma}}}_{1,{\varvec{j}}{\varvec{b}}{\varvec{k}}}^{{{\varvec{R}}}_{{\varvec{j}}}}+{{\varvec{\delta}}}_{{\varvec{j}}{\varvec{b}}{\varvec{k}}}+{{{\varvec{\gamma}}}_{2,{\varvec{j}}{\varvec{b}}{\varvec{k}}}{\varvec{R}}}_{{\varvec{j}}}+{\beta }_{1bk}^{B} {\overline{x} }_{j.}\)
\({R}_{j}\) is sampled from a Bernoulli distribution with bias probability \({\pi }_{j}\)_{,} for each study. Then \({\pi }_{j}\) is assigned a beta distribution; \({\pi }_{j}\sim Beta\left({a}_{1},{a}_{2}\right)\) where the values of \({a}_{1}\) and \({a}_{2}\) reflect the RoB (low, high or unclear) within RCTs or NRS. The ratio \({a}_{1}/{a}_{2}\) controls the skewness of the beta distribution. When the ratio \({a}_{1}/{a}_{2}\) approaches 1, the mean probability of bias gets closer to 1 and the study acquires ’major’ bias adjustment. The default beta priors in the package are: high bias RCT prior.pi.high.rct = "dbeta(10,1)", low bias RCT prior.pi.low.rct = "dbeta(1,10)", high bias NRS prior.pi.high.nrs = "dbeta(30,1)" and low bias NRS prior.pi.low.nrs = "dbeta(1,30)". Alternatively, we can use study characteristics \({z}_{j}\) to predict \({\pi }_{j}\) through a logistic transformation (internally coded); \(logit({\pi }_{j})=e+f{z}_{j}\). When \({z}_{j}\) is a continuous outcome, \({\text{exp}}(e)\) is the odds of bias at \({z}_{j}=0\) and \({\text{exp}}(f)\) is the odds ratio of bias for a one unit increase in \({z}_{j}\). When \(f\) has a positive value, the bias probability increases with increasing values of \({z}_{j}\).
Table 2 shows how we combine the multiplicative \({(\gamma }_{1,jbk})\) and the additive (\({\gamma }_{2,jbk})\) treatmentspecific bias effects across studies using randomeffects or commoneffect models.
Biasadjusted model 2
Another way to account for differences in RoB of the studies is to replace \({\delta }_{jbk}\) with a biasadjusted relative treatment effect \({\theta }_{jbk}\) (method.bias = "adjust2") [23].
The equations become for level 1 and 2 are
and
Then \({\theta }_{jbk}\) is given either a randomeffect bimodal normal distribution; \({\theta }_{jbk}\sim \left(1{\pi }_{j}\right)\mathrm{\rm N}\left({d}_{Ak}{d}_{Ab},{\tau }^{2}\right)+{\pi }_{j}\mathrm{\rm N}\left({d}_{Ak}{d}_{Ab}+{\gamma }_{jbk},{\tau }^{2}+{\tau }_{\gamma }^{2}\right)\) or assumed common across studies; \({\theta }_{jbk}={d}_{Ak}{d}_{Ab}+{\pi }_{j}{\gamma }_{jbk}.\)
The likelihood of the unknown parameter \({\theta }_{jbk}\) is
Here, the bias probability \({\pi }_{j}\) determines the weight of the biasadjusted distribution (second part of the equation) in the overall likelihood \(L\left({\theta }_{jbk}; {d}_{Ak},{d}_{Ab},\tau ,{\tau }_{\gamma },{\pi }_{j}\right)\). The term \({\gamma }_{jbk}\) is the bias effect, as in biasadjusted model 1.
Implementation
We implement the Bayesian models in a new R package called crossnma. The user can install the package with the command install.packages("crossnma") and then load the library into the current R session with library("crossnma").
The Bayesian model is run in the background using Just Another Gibbs Sampler (JAGS) software [24]. Therefore, the JAGS programme must be installed on the user’s local computer (see https://sourceforge.net/projects/mcmcjags/). A vignette with a binary data example is part of crossnma which can be opened using vignette(“crossnma”). Package updates providing new features or fixing bugs will be posted on the package website: https://github.com/htxr/crossnma.
Workflow within crossnma
Figure 1 presents the workflow for conducting analyses within crossnma. Before running crossnma(), we display the network of evidence using netgraph() (which is a generic function in netmeta) to display the network of evidence. To conduct the data synthesis, there are two main steps: use crossnma.model() to produce the JAGS code and reformat the data, then pass the output to crossnma(), which matches the data with the model, runs the analysis and estimates all model parameters. The generic function plot() can produce a trace plot to evaluate the Markov chain Monte Carlo (MCMC) convergence for each model parameter. The functions summary(), league() and heatplot() can be used, with the output of crossnma() as input, to produce a numerical and graphical summaries of the treatment effect estimates. More details on how to use of the functions and their arguments can be found in Supplementary Document S1.
Comparison with some of the available packages
We compare the output of crossnma (version 1.2.0) with BUGSnet (version 1.1.0), gemtc (version 1.0.2), multinma (version 1.1.2) bnma (version 1.5.1) and netmeta (version 2.82) concerning various features in Table 3. Additionally, we assess the performance of some of these packages using a dataset provided by crossnma, which we describe in the next section. As BUGSnet and gemtc can only synthesize aggregate data, we summarize IPDs at the arm level (Supplementary Data S1). Then we perform NMA with a randomeffects model using all packages. Treatment effect estimates (odds ratios) from crossnma and the three other packages do not differ beyond the MC error, however, the BUGSnet estimate of the betweenstudy variance (\(\tau\)) is substantially larger compare to other packages. This can be attributed to the fact that BUGSnet uses an unrealistic default prior distribution \(\sigma \sim Unif(\mathrm{0,1.62})\) where \(\tau =1/{\sigma }^{2}\). In crossnma, we set \(\tau \sim Unif(\mathrm{0,1.62})\). R code and detailed output of our analyses are provided in Supplementary Data S2, Supplementary Tables S1, S2, S3 and S4.
We performed NMR for age using crossnma and multinma. Supplementary Table S5 shows the estimated treatment effects. The disagreement in the estimation is due to the differences in the implemented models and variations in each package's builtin analysis settings. The code for conducting this comparative analysis is provided in Supplementary Data S4.
Working example
In the following, we illustrate NMA and NMR in crossnma. We analyse fictitious data, simulated to mimic real RCTs with IPD and AD [25]. The code for each analysis is available in Supplementary Data S3 (as well as presented in the vignette).
Description of the network
The evidence network consists of four drugs examined in six studies, with aggregate data (two RCTs) or individual participant data (three RCTs and one cohort study). The IPD dataset contains 1944 rows, i.e., study participants. The AD dataset is provided in armlevel format, with each row representing a study arm and the same variable names as the IPD dataset. Below, we present the two datasets: the first few rows of the IPD dataset and the complete set of rows for the AD dataset. We evaluate the treatment effect using a binary outcome of relapse after two years of followup. The relative treatment effects are expressed as odds ratios, where an OR below 1 indicates the treatment is preferable to the reference.
Data synthesis using the four models available in crossnma
We continue with the analysis as the network is connected. We begin with creating a JAGS model and reformatting both datasets using crossnma.model(). Then, as the network is connected, the output of crossnma.model() is passed to crossnma(), which runs NMA with MCMC for 5000 iterations, 2000 burnin, one thinning and two chains (default settings). In this example, we use drug A as a reference treatment except in the analysis in Using nonrandomized studies as a prior in network metaregression section drug D is the reference.
Unadjusted network metaanalysis
First, we synthesize data from RCT and NRS without distinguishing between them (method.bias = 'naive'). Because there are few studies in the network, we expect the heterogeneity parameter to be estimated inefficiently and thus assume a more informative prior to improve estimation, \(\tau \sim N(0, 1/3)\). The data is analyzed using odds ratio as a summary measure sm = "OR"(which is the default for binary outcomes). By choosing trt.effect = 'random' (default), we are assigning a normal distribution to each relative treatment effect to allow the synthesis across studies, Table 2 lists all supported options.
We can also compute the values of Surface Under the Cumulative Ranking (SUCRA) (by enabling the sucra = TRUE option), but it's essential to specify a negative preferred direction for the outcome (using small.values = "desirable"). This setting indicates that lower values of the relative treatment effect signify the treatment's effectiveness. Conversely, if positive values are preferred, you can set small.values = "undesirable".
The network graph in Fig. 2 was generated with the following command:
In the graph, the thickness of the edges corresponds to the number of studies, while the number of studies is displayed on the edges. Additionally, the node size reflects the number of participants who received each intervention.
Next, we fit the NMA model using crossnma() where we can set the number of iterations, burnin, thinning and chains. We run all subsequent models under the same settings. Note that for our example, we run the MCMC with settings different from the default to ensure convergence.
The R command print(jagsfit1, backtransf = FALSE) produces Table 4a which shows summary statistics for the results (Table 4be are also produced using the print() function for other models). The estimated OR of B vs A can be obtained as exp(d.B), and similarly exp(d.C) and exp(d.D) are the ORs of C and D relative to A.
Function print() also produces the SUCRA rank estimates, where treatment D notably excels with the highest score of 0.941, signifying a strong likelihood of achieving favorable outcomes. In contrast, treatment A has the lowest score at 0.007, implying that it is the least likelihood choice for yielding positive outcomes.
Unadjusted network metaregression
Next, we include age in the model as a potential effect modifier. Because we have few studies and little variation between them, we assume that the age effect within and between studies is equal (argument split.regcoef = FALSE).
In addition to relative treatment effects and its heterogeneity, we obtain estimates of the age effect (b_1 is \({B}_{1}^{W}={B}_{1}^{B}\)) and the heterogeneity standard deviation (tau.b_1 is \({\tau }_{B}={\tau }_{W}\)) in the effect of age across studies (Table 4b).
We could add two more covariates to the NMR model using arguments cov2 and cov3.
The MCMC is run under the same set up as in Unadjusted network metaanalysis section. The league table of the estimates of each treatment vs comparator can be shown as follows (we just display the first two lines of the output).
Using nonrandomized studies as a prior in network metaregression
We use the single NRS study to construct priors for a subsequent NMA fitted on RCT data. The prior variance is inflated by 60% to reduce the contribution of NRS on the final estimation (\(w=0.6\)). Table 4c presents the estimates where the reference treatment is drug D as drug A is not examined in NRS.
The heat plot (Fig. 3) summarizes the relative effect with the 95% credible interval of each treatment on the top compared to the treatment on the left. All estimates are computed for participant age 38.
Biasadjusted model 1
To fit this model, data needs to include studylevel RoB data and indicate the direction of bias (which treatment in the study is expected to be favoured). We assume that additive bias effects are equal across studies. We estimate the probability of bias using the year of study publication.
The common bias effect (g (R output)) is estimated to be 0.112, indicating that older studies tend to overestimate the relative treatment effect when compared to new studies (see Table 4d). We note that we obtained very uncertain estimates for the mean bias effect. This is because the dataset includes only three studies at low and three studies at high RoB. In the presence of more studies, estimation improves both in convergence and precision as shown in Hamza et al. [26].
Biasadjusted model 2
This analysis requires data similar to biasadjusted model 1. We use the default beta priors to estimate the bias probabilities. The overall bias effect (g (R output)=\({g}_{m}\) (model description)) is estimated to be 0.016 (Table 4e), implying that studies with a high RoB slightly underestimate treatment effect when compared to studies with a low RoB but this estimate again comes with large uncertainty.
Models convergence
To evaluate the convergence of the MCMC chains of all models, we use the Gelman and Rubin statistic \(\widehat{R}\) and the number of effective sample sizes (n.eff) shown in Table 4 [27]. Except for the common bias effect (g (R output)) in biasadjusted model 1, \(\widehat{R}\) values are approximately 1. The values of n.eff indicates that sufficient independent samples are used to generate the final estimates. We inspect the trace plot (generated by plot(jagsfit4)) in Fig. 4 to further investigate the convergence of g (besides other parameters of the biasadjusted model 1), and we observe a great deviation between samples. This is because the biasadjusted model 1 includes the bias as a dichotomous variable, which requires having sufficient data at both low and high RoB. The dataset we analyze does not contain enough of this data (3 studies at low and 3 at high RoB).
In Fig. 5, we present density plots to visualize the distributions of the variables derived from the MCMC samples. The density plots illustrate that most variables exhibit normal distributions, reflecting symmetrical data with a clear central tendency. However, in the case of \(\tau\), a positiveonly normal distribution is observed, indicating values restricted to the positive range. Notably, the density distribution of \(g\) demonstrates a wide range of values, suggesting insufficient data for precise estimation. This finding underscores the importance of acquiring additional data to improve the accuracy and reliability of the estimation process.
Computational efficiency of crossnma
The five analyses were conducted using the crossnma package (version 1.2.0) in R (version 4.2.3) on a MacBook Pro (13inch, 2019, Two Thunderbolt 3 ports) with a 1.4 GHz QuadCore Intel Core i5 processor and 8 GB 2133 MHz LPDDR3 RAM. All analyses include 4 studies with IPD and 1944 participants and 2 studies with AD and 4 treatment arms, with a total of 4 treatments. The runtime of the analyses done in this article varied between 2 to 5 min, depending on the specific analysis. The longest runtime of 5 min was observed for network metaregressions that included a single covariate (Unadjusted network metaregression and Using nonrandomized studies as a prior in network metaregression sections). The two biasadjustment models (Biasadjusted model 1 and Biasadjusted model 2 sections) took approximately 3 and a half minutes. The shortest runtime of 2 min was observed for NMA (Unadjusted network metaanalysis section) without adjustment for bias. For all analyses we run 100,000 iterations with a burnin of 40,000 and thinning of 5 to ensure convergence.
Discussion
In this paper, we introduce crossnma, an R package that performs Bayesian NMA and NMR using the JAGS software. In crossnma, data can be collected from different study designs, as RCT or NRS, and provided in IPD or AD formats. The functions within the package enable analysis, result representation and convergence evaluation. We provide detailed instructions on how to use crossnma and we demonstrate this with several analytic examples.
Several R packages for performing NMA on aggregate data are available, such as gemtc [19], bnma [28] and BUGSnet [29] in a Bayesian setting, or netmeta [30] under a frequentist framework. However, data is increasingly becoming available in a variety of formats and designs. For example, there is a growing in the number of IPD analyses, and only userwritten code can be used to perform such analyses. The number of reviews that combine NRS and RCT data is rising as well, and unadjusted synthesis is widely used due to its simplicity and ease of implementation [31]. Network metaregression using only aggregated data can be performed with bnma, rnmamod, gemtc and BUGSnet. Using the methodologies presented by Philippo et al. [32] the R package multinma models jointly effects estimated in studies with IPD and AD formats. The package crossnma implemnets another methodology to merge estimates both formats [26]. Furthermore, crossnma can perform sensitivity analyses for studylevel bias in RCT and NRS data [32]. Some functions in crossnma are similar to and inspired by the Bayesian NMA packages gemtc and BUGSnet. In addition to handling both AD and IPD, the crossnma package can be used to account for various levels of study risk of bias and their impact in the results.
In crossnma, we implemented, among others, a model for the synthesis of IPD and AD data previously described in [16, 17, 26] assuming that the relative treatment effect in IPD and AD synthesis models are the same after accounting for effect modifiers and prognostic factors. A method that makes less assumptions and with theoretical advantages has been presented in [32]. Future updates of crossnma could include an option of the model currently implemented in [33].
Regarding the inclusion of NRS data, it is important to note that in situations where RCT data is unavailable for certain comparisons or treatment interventions, researchers may have to rely on NRS data to inform the analysis. However, it is crucial to recognize that NRS data are susceptible to various biases of unknown magnitude, which necessitates careful consideration when utilizing them in the analysis. In the crossnma package we implemented methods that can decrease the impact that such biases may have in the final NMA results.
While crossnma allows effect modifying covariates to work in a range of ways, including a different regression coefficient for each treatment, data might not enable their estimation. Hence, in practice users are more likely to employ the model that assumes the same regression coefficient for all treatments.
Several limitations should be acknowledged regarding the statistical approaches implemented in crossnma. First, incorporating NRS evidence as prior in the analysis of RCTs can be complicated in practice. Collecting expert opinion about the bias in NRS is time consuming and often impractical. The use of priors from NRS should be implemented via a sensitivity analysis using a range of “downweighing” values for the impact of the prior in the results of NMA. Second, the model by Verde includes a parameter for the probability of bias, which is difficult to estimate from the current data, so informative priors are required. To establish these priors as subjectively as possible, trained data extractors are needed to evaluate the risk of bias in each study using established and reproducible tools, like RoB2 and ROBINSI. Third, the identifiability of all model parameters, and in particular those that relate to bias, depend on the available data. Fourth, sensitivity to the choice of prior distribution necessitates conducting thorough sensitivity analysis. While we provide recommendations in our recent paper [26], further research is needed to explore alternative methods and enhance the applicability of biasadjustment techniques in decisionmaking contexts. Fifth, the implemented models for the synthesis of AD and IPD are an approximations of the model implemented in multinma [33], and their performance is unknown. These models may not be readily generalizable for timetoevent outcomes. Future updates to the package will incorporate the models described in [32], thereby overcoming these limitations. Finally, crossnma assumes similar relative treatment effects in IPD and AD, which holds true only for noncollapsible outcomes. However, for noncollapsible outcomes like logit, this assumption introduces aggregation bias [6].
The model to combine IPD and AD implemented in crossnma assumes distinct regression coefficients for interaction terms at the IPD and AD levels. In contrast, the integration approach implemented in multinma do not require ADspecific interaction terms, as these are inherently defined by the integration process. The models implemented in crossnma can be viewed as an approximation to the models by Philippo et al. [32] which have a theoretical advantage. However, application of the latter model requires additional data or assumptions to establish the correlation structure between covariates, which can be challenging in practice. A largescale comparison of these two modelling approaches using realistic scenarios would shed more light to the impact of model misspecification, violation of model assumptions and extend of aggregation bias.
In addition to the foundational assumptions that underlie conventional metaanalysis (e.g., the assumption that treatment effects are generalizable across patients from the included trials) and the assumptions inherent in NMA (e.g., a connected network and consistency of effects), all metaregression models assume the absence of unobserved effect modifiers [16,17,18, 22, 23]. IPD network metaregression models rely on the assumption of conditional constancy of relative effects, which asserts that relative effects remain constant across different populations at specific levels of a set of covariates.
We acknowledge the following crossnma shortcomings. In addition to the current functions that generate league tables and summary statistics, we plan to develop new functions to present results for the following purposes: displaying the distribution of potential effect modifiers by study, treatment, or both; presenting SUCRA scores as plots and tables to enable ranking treatments, and producing a plot of the estimates of relative treatment effects at various covariate values (for NMR model). Our package supports binary and continuous outcomes, analysed in the vast majority of published NMAs [34]. Future updates will include count and timetoevent outcomes. Also, we plan to develop a separate vignette that focuses specifically on continuous outcomes. These additional features and resource will provide users with a more comprehensive understanding of the package's versatility and how it can be applied in various analysis scenarios. In terms of summary measure, crossnma enables expressing relative treatment effects in terms of odds ratio or risk ratio for binary data and mean difference or standardised mean difference for continuous outcomes.
The data in crossnma must be provided at the arm level which may require additional data manipulation. For example, contrastlevel data can be transformed to the armlevel format using R function longarm() from R package meta. A future extension will expand this to provide contrastlevel data directly. Methods for evaluating inconsistency, including node splitting and unrelated mean effect, are not yet implemented in crossnma. We intend to address these issues in upcoming version of crossnma.
The package crossnma should be used in conjunction with the technical article that describe the models, their assumptions, and limitations. The vignette accompanying the crossnma package but mainly the publication by Hamza et al. [26] contains useful information to enable users to set up models that sensibly reflect the nature of their data. Users are advised in particular to pay close attention to the assumptions behind the models, which are described in [26].
Conclusions
The R package crossnma enables the user to perform NMA and NMR with different data types in a Bayesian framework and facilitates the inclusion of all types of evidence accounting for their differences in risk of bias.
Availability and requirements
Project name: crossnma project
Project home page: https://github.com/htxr/crossnma
Operating system(s): Any OS providing R and JAGS
Programming language: R
Other requirements: JAGS 4.3.0
License: GNU GPL2 or higher versions
Any restrictions to use by nonacademics: no restrictions
Availability of data and materials
Project name: crossnma project
Project home page: https://github.com/htxr/crossnma
Operating system(s): Any OS providing R and JAGS
Programming language: R
Other requirements: JAGS 4.3.0
License: GNU GPL2 or higher versions
Any restrictions to use by nonacademics: no restrictions
Abbreviations
 NMA:

Network MetaAnalysis
 NMR:

Network MetaRegression;
 IPD:

Individual Participant Data
 AD:

Aggregate Data
 RCT:

Randomised clinical trials
 NRS:

NonRandomised Studies
 RoB:

Risk of Bias
 JAGS:

Just Another Gibbs Sampler
 MCMC:

Markov chain Monte Carlo
References
Higgins JPT, Thomas J, Chandler J, Cumpston M, Li T, Page M, et al. Cochrane handbook for systematic reviews of interventions. 2nd ed. Chichester: Wiley; 2019.
Borenstein M, Hedges LV, Higgins JPT, Rothstien HR. Introduction to meta‐analysis. Wiley; 2009. Available from: https://onlinelibrary.wiley.com/doi/abs/10.1002/9780470743386.ch1. Cited 2022 Feb 4.
Borenstein M, Hedges L, Higgins J, Rothstien H. A basic introduction to fixedeffect and randomeffects models for metaanalysis. Res Synth Methods. 2010;1:97–111.
Caldwell DM, Ades AE, Higgins JPT. Simultaneous comparison of multiple treatments: combining direct and indirect evidence. BMJ. 2005;331(7521):897–900.
Dias S, Sutton AJ, Welton NJ, Ades AE. Evidence synthesis for decision making 3: heterogeneity–subgroups, metaregression, bias, and biasadjustment. Med Decis Making. 2013;33(5):618–40.
Rothman KJ, Greenland S, Lash TL. Modern epidemiology. 3rd ed. Philadelphia: Wolters Kluwer Health; 2008.
Chalmers TC, Celano P, Sacks HS, Smith H. Bias in treatment assignment in controlled clinical trials. N Engl J Med. 1983;309(22):1358–61.
Schulz KF, Chalmers I, Hayes RJ, Altman DG. Empirical evidence of bias. Dimensions of methodological quality associated with estimates of treatment effects in controlled trials. JAMA. 1995;273(5):408–12.
Wood L, Egger M, Gluud LL, Schulz KF, Jüni P, Altman DG, et al. Empirical evidence of bias in treatment effect estimates in controlled trials with different interventions and outcomes: metaepidemiological study. BMJ. 2008;336(7644):601–5.
Bero L, Oostvogel F, Bacchetti P, Lee K. Factors associated with findings of published trials of drugdrug comparisons: why some statins appear more efficacious than others. PLoS Med. 2007;4(6):e184.
Flacco ME, Manzoli L, Boccia S, Capasso L, Aleksovska K, Rosso A, et al. Headtohead randomized trials are mostly industry sponsored and almost always favor the industry sponsor. J Clin Epidemiol. 2015;68(7):811–20.
Gazendam AM, SlawaskaEng D, Nucci N, Bhatt O, Ghert M. The impact of industry funding on randomized controlled trials of biologic therapies. Medicines (Basel). 2022;9(3):18.
Sterne JA, Hernán MA, Reeves BC, Savović J, Berkman ND, Viswanathan M, et al. ROBINSI: a tool for assessing risk of bias in nonrandomised studies of interventions. BMJ. 2016;355:i4919.
Sterne JAC, Savović J, Page MJ, Elbers RG, Blencowe NS, Boutron I, et al. RoB 2: a revised tool for assessing risk of bias in randomised trials. BMJ. 2019;366:l4898.
Hamza T, Chalkou K, Pellegrini F, Kuhle J, Benkert P, Lorscheider J, et al. Synthesizing crossdesign evidence and crossformat data using network metaregression. 2022. Available from: https://arxiv.org/abs/2203.06350. Cited 2022 Dec 16.
Saramago P, Sutton AJ, Cooper NJ, Manca A. Mixed treatment comparisons using aggregate and individual participant level data. Stat Med. 2012;31(10970258 (Electronic)):3516–36.
Donegan S, Williamson P, D’Alessandro U, Garner P, Smith CT. Combining individual patient data and aggregate data in mixed treatment comparison metaanalysis: Individual patient data may be beneficial if only for a subset of trials. Stat Med. 2013;32(10970258 (Electronic)):914–30.
Jansen JP. Network metaanalysis of individual and aggregate level data. Res Synth Meth. 2012;3(2):177–90.
van Valkenhoef G, Lu G, de Brock B, Hillege H, Ades AE, Welton NJ. Automating network metaanalysis. Res Syn Meth. 2012;3(4):285–99.
Efthimiou O, Mavridis D, Debray TPA, Samara M, Belger M, Siontis GCM, et al. Combining randomized and nonrandomized evidence in network metaanalysis. Stat Med. 2017;36(8):1210–26.
Turner RM, Spiegelhalter DJ, Smith GCS, Thompson SG. Bias modelling in evidence synthesis. J R Stat Soc Ser A Stat Soc. 2009;172(1):21–47.
Dias S, Welton NJ, Marinho VCC, Salanti G, Higgins JPT, Ades AE. Estimation and adjustment of bias in randomized evidence by using mixed treatment comparison metaanalysis. J R Stat Soc Ser A Stat Soc. 2010;173:613–29.
Verde PE. A biascorrected metaanalysis model for combining, studies of different types and quality. Biom J. 2021;63(2):406–22.
Plummer M. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. 2003.
Tramacere I, Del Giovane C, Salanti G, D’Amico R, Filippini G. Immunomodulators and immunosuppressants for relapsing‐remitting multiple sclerosis: a network meta‐analysis. Cochrane Database Syst Rev. 2015(9):CD011381. Available from: https://www.cochranelibrary.com/cdsr/doi/10.1002/14651858.CD011381.pub2/full.
Hamza T, Chalkou K, Pellegrini F, Kuhle J, Benkert P, Lorscheider J, et al. Synthesizing crossdesign evidence and crossformat data using network metaregression. Res Synth Methods. 2023;14(2):283–300.
Gelman A, Carlin JB, Stern HS, Rubin DB. Bayesian data analysis. 2nd ed. Boca Raton: Chapman & Hall/CRC; 2004.
Seo M, Schmid C. bnma: Bayesian network metaanalysis using “JAGS”. 2022. Available from: https://CRAN.Rproject.org/package=bnma.
Béliveau A, Boyne DJ, Slater J, Brenner D, Arora P. BUGSnet: an R package to facilitate the conduct and reporting of Bayesian network metaanalyses. BMC Med Res Methodol. 2019;19(1):196.
Schwarzer G, Carpenter JR, Rücker G. Network metaanalysis. In: Schwarzer G, Carpenter JR, Rücker G, editors. Metaanalysis with R. Cham: Springer International Publishing; 2015. p. 187–216. https://doi.org/10.1007/9783319214160_8. (Use R!). Cited 2022 Dec 7.
Zhang K, Arora P, Sati N, Béliveau A, Troke N, Veroniki AA, et al. Characteristics and methods of incorporating randomized and nonrandomized evidence in network metaanalyses: a scoping review. J Clin Epidemiol. 2019;113:1–10.
Phillippo DM, Dias S, Ades AE, Belger M, Brnabic A, Schacht A, et al. Multilevel network metaregression for populationadjusted treatment comparisons. J R Stat Soc Ser A Stat Soc. 2020;183(3):1189–210.
Phillippo DM. multinma: Bayesian network metaanalysis of individual and aggregate data. R package version 0.5.1. 2023. https://doi.org/10.5281/zenodo.3904454. https://dmphillippo.github.io/multinma/.
Papakonstantinou T, Nikolakopoulou A, Egger M, Salanti G. In network metaanalysis, most of the information comes from indirect evidence: empirical study. J Clin Epidemiol. 2020;1(124):42–9.
Acknowledgements
We would also like to thank Jiu Li for providing valuable feedback regarding the implementation of crossnma.
Funding
TH and GSa are funded by the HTx project which has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement Nº 825162. This dissemination reflects only the author's view and the Commission is not responsible for any use that may be made of the information it contains.
Author information
Authors and Affiliations
Contributions
TH prepared an initial draft of the manuscript and run the analysis. TH, GSa and GSc contributed to the conceptualization, writing, and editing of the manuscript. TH, GSa and GSc read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 2: Supplementary Data S1.
This dataset is formed by merging both Aggregate Data (AD) and Individual Participant Data (IPD) from crossnma. This dataset is employed in the R code found in Supplementary Data S2 to perform a comparative analysis among the packages BUGSnet, gemtc, multinma, and crossnma.
Additional file 4: Supplementary Table S1.
Summary statistics of network metaanaylsis estimates produced by BUGSnet package, applied to the data in Supplementary Data S1.
Additional file 5: Supplementary Table S2.
Summary statistics of network metaanaylsis estimates produced by crossnma package, applied to the data in Supplementary Data S1.
Additional file 6: Supplementary Table S3.
Summary statistics of network metaanaylsis estimates produced by gemtc package, applied to the data in Supplementary Data S1.
Additional file 7: Supplementary Table S4.
Summary statistics of network metaanaylsis estimates produced by multinma package, applied to the data in Supplementary Data S1.
Additional file 8: Supplementary Table S5.
Summary statistics of network metaregression estimates produced by multinma and crossnma packages. These estimates are based on data available in crossnma, utilizing both Individual Participant Data (IPD) and Aggregate Data (AD).
Additional file 11.
Reviewer reports submitted for version 1.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.
About this article
Cite this article
Hamza, T., Schwarzer, G. & Salanti, G. crossnma: An R package to synthesize crossdesign evidence and crossformat data using network metaanalysis and network metaregression. BMC Med Res Methodol 24, 169 (2024). https://doi.org/10.1186/s12874023021300
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12874023021300