Class effect NMA analysis using MBNMAdose

Although MBNMAdose is intended to be used for dose-response Model-Based Network Meta-Analysis (MBNMA), it can also be adapted to perform standard Network Meta-Analysis (NMA), and this allows users to take advantage of some of the additional features of MBNMAdose, such as modelling class effects, for use in standard NMA. As well as fitting class effect models, MBNMAdose also allows for nodes-splitting to check for consistency in these models.

To illustrate how this can be done we will use a dataset of inhaled medications for Chronic Obstructive Pulmonary Disease (COPD) from the netmeta package:

library(netmeta)
#> Loading required package: meta
#> Loading required package: metabook
#> Loading 'meta' package (version 8.5-0).
#> Type 'help(meta)' for a brief overview.
#> Loading 'netmeta' package (version 3.5-0).
#> Type 'help("netmeta-package")' for a brief overview.
data("Dong2013")

# Rename column names to match those used in MBNMAdose
Dong2013 <- Dong2013 %>%
    rename(studyID = id, r = death, n = randomized)

Performing standard NMA using mbnma.run()

The simplest use is in a network that includes a placebo treatment. In this dataset we do not have any dose-response information, so there is no value in performing a MBNMA. However, if we assume that every active treatment in the network is a separate “agent” with a dose of 1, and that the Placebo treatment has a dose of 0, then we can use the data in MBNMAdose, and by modelling a linear dose-response function we estimate parameters that are identical to a standard NMA model.

# Define agents and assign a dose of 1 to all agents
Dong2013 <- Dong2013 %>%
    dplyr::rename(agent = treatment) %>%
    dplyr::mutate(dose = dplyr::case_when(agent == "Placebo" ~ 0, agent != "Placebo" ~
        1))

Note that if there is an intervention within the dataset that has been administered at multiple doses, you can force the dataset to be analysed either as a “split” network (in which different doses are assumed to have independent effects) by assigning each of them a separate agent name (e.g. “warfarinlow”, “warfarinhigh”), or as a “lumped” network (in which different doses are assumed to have the same effect) by simply assigning both doses a dose of 1. Further details of “lumping” and “splitting” and the implications of these assumptions can be found in [@pedder2021cons].

Once we have reassigned the doses within the dataset we can create an "mbnma.network" object and create a network plot:

network <- mbnma.network(Dong2013)
summary(network)
#> Description: Network 
#> Number of studies: 41 
#> Number of treatments: 6 
#> Number of agents: 6 
#> Median (min, max) doses per agent (incl placebo): 2 (2, 2)
#> Agent-level network is CONNECTED
#> Ttreatment-level network is CONNECTED
plot(network)

We can then use a linear dose-response MBNMA to analyse the data. The coefficients for the linear slope of the dose-response function are mathematically equivalent to the basic treatment effect parameters estimated in a standard NMA. Note that the results are equivalent in both models (allowing for Monte-Carlo error from the MCMC sampling). The only difference is that the placebo parameter (beta.1[1]) that is equal to zero is not given in the output.

nma.linear <- mbnma.run(network, fun = dpoly(degree = 1), n.iter = 50000)
#> `likelihood` not given by user - set to `binomial` based on data provided
#> `link` not given by user - set to `logit` based on assigned value for `likelihood`
print(nma.linear)
#> Inference for Bugs model at "/tmp/Rtmpwi5mW9/file140a5967b3c2", fit using jags,
#>  3 chains, each with 50000 iterations (first 25000 discarded), n.thin = 25
#>  n.sims = 3000 iterations saved. Running time = 6.633 secs
#>           mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
#> beta.1[2]   0.027   0.085  -0.139  -0.029   0.026   0.083   0.193 1.002  1200
#> beta.1[3]  -0.082   0.081  -0.244  -0.136  -0.083  -0.028   0.081 1.001  3000
#> beta.1[4]  -0.231   0.088  -0.399  -0.292  -0.230  -0.171  -0.058 1.002  3000
#> beta.1[5]  -0.086   0.067  -0.217  -0.133  -0.085  -0.041   0.047 1.001  2400
#> beta.1[6]   0.426   0.184   0.070   0.300   0.427   0.547   0.787 1.003   810
#> totresdev 108.696   9.622  92.138 101.668 108.116 114.833 129.171 1.005   440
#> deviance  414.395   9.622 397.837 407.368 413.815 420.533 434.871 1.005   430
#> 
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> 
#> DIC info (using the rule: pV = var(deviance)/2)
#> pV = 46.1 and DIC = 460.5
#> DIC info (using the rule: pD = Dbar-Dhat, computed via 'rjags::dic.samples')
#> pD = 43.7 and DIC = 456.1
#> DIC is an estimate of expected predictive error (lower deviance is better).
nma <- nma.run(network, n.iter = 50000)
print(nma)
#> $jagsresult
#> Inference for Bugs model at "/tmp/Rtmpwi5mW9/file140a19c50a8f", fit using jags,
#>  3 chains, each with 50000 iterations (first 25000 discarded), n.thin = 25
#>  n.sims = 3000 iterations saved. Running time = 6.177 secs
#>           mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
#> d[1]        0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
#> d[2]        0.031   0.083  -0.130  -0.025   0.031   0.086   0.193 1.002  1000
#> d[3]       -0.074   0.082  -0.230  -0.130  -0.073  -0.018   0.085 1.001  3000
#> d[4]       -0.232   0.088  -0.404  -0.292  -0.231  -0.172  -0.059 1.001  3000
#> d[5]       -0.087   0.065  -0.215  -0.130  -0.086  -0.044   0.040 1.001  3000
#> d[6]        0.412   0.187   0.064   0.279   0.409   0.539   0.781 1.010   210
#> totresdev 108.958   9.710  92.217 102.086 107.964 115.112 129.675 1.002  2000
#> deviance  414.658   9.710 397.917 407.786 413.664 420.812 435.375 1.002  1900
#> 
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> 
#> DIC info (using the rule: pV = var(deviance)/2)
#> pV = 47.1 and DIC = 461.8
#> DIC is an estimate of expected predictive error (lower deviance is better).
#> 
#> $trt.labs
#> [1] "Placebo_0"  "ICS_1"      "LABA_1"     "LABA-ICS_1" "TIO-HH_1"  
#> [6] "TIO-SMI_1" 
#> 
#> $UME
#> [1] FALSE
#> 
#> attr(,"class")
#> [1] "nma"

We can also show the equivalence of results using get.relative() to compare relative effects from both models:

rels <- get.relative(nma.linear, nma)

Without a placebo, estimation is very similar, but requires renaming and recoding the network reference intervention to "Placebo". This is not strictly necessary, as MBNMAdose can handle dose-response datasets that do not include placebo (or dose=0), but it will ensure that parameter estimates are equivalent between the NMA and MBNMA models and will make it easier to estimate relative effects. We illustrate this with the Surgical Site Infection dataset.

# Ensure that Suture-absorbable is the network reference
ssi <- ssi_closure %>%
    dplyr::mutate(agent = factor(trt, levels = c("Suture-absorbable", unique(ssi_closure$trt)[-1])))

# Set dose=0 for network reference and dose=1 for all other interventions
ssi.plac <- ssi %>%
    dplyr::mutate(dose = dplyr::case_when(trt == "Suture-absorbable" ~ 0, TRUE ~
        1))

network.plac <- mbnma.network(ssi.plac)
#> Values for `agent` with dose = 0 have been recoded to `Placebo`
#> agent is being recoded to enforce sequential numbering
#> Values for `class` with dose = 0 have been recoded to `Placebo`
#> class is being recoded to enforce sequential numbering
plot(network.plac)

# Note that Suture-absorbable (the comparator) has been renamed to Placebo
# Run linear MBNMA model
nma.linear <- mbnma.run(network.plac, fun = dpoly(degree = 1), n.iter = 50000)
#> `likelihood` not given by user - set to `binomial` based on data provided
#> `link` not given by user - set to `logit` based on assigned value for `likelihood`
summary(nma.linear)
#> ========================================
#> Dose-response MBNMA
#> ========================================
#> 
#> Likelihood: binomial
#> Link function: logit
#> Dose-response function: poly
#> 
#> Pooling method
#> 
#> Method: Common (fixed) effects estimated for relative effects
#> 
#> 
#> beta.1 dose-response parameter results
#> 
#> Pooling: relative effects for each agent
#> 
#> |Agent                               |Parameter  |  Median|    2.5%|   97.5%|
#> |:-----------------------------------|:----------|-------:|-------:|-------:|
#> |Clips                               |beta.1[2]  | -0.0910| -1.0239|  0.8203|
#> |Suture--monofilament                |beta.1[3]  |  0.0333| -0.8322|  0.8193|
#> |Staples                             |beta.1[4]  |  0.0394| -0.5341|  0.6133|
#> |Suture                              |beta.1[5]  | -0.0566| -0.7231|  0.6793|
#> |Tissue adhesive gel-low viscosity   |beta.1[6]  |  0.6430| -0.0962|  1.4208|
#> |Tape/film                           |beta.1[7]  | -0.3082| -1.0852|  0.3651|
#> |Preference of physician             |beta.1[8]  | -0.2220| -1.5997|  1.0016|
#> |Tissue adhesive gel-high viscosity  |beta.1[9]  | -0.1908| -1.9480|  1.4476|
#> |Suture-non-absorbable               |beta.1[10] | -0.1055| -0.3872|  0.1765|
#> |Suture-antimicrobial                |beta.1[11] | -0.2713| -0.4259| -0.1179|
#> |Suture-non-absorbable-monofilament  |beta.1[12] | -0.6423| -1.2468| -0.0030|
#> |Suture-absorbable-monofilament      |beta.1[13] | -0.5101| -1.1469|  0.1512|
#> |Tissue adhesive gel                 |beta.1[14] |  0.9403| -1.6874|  5.1848|
#> |Suture-absorbable-multifilament     |beta.1[15] | -0.7533| -1.3815| -0.0827|
#> |Suture-non-absorbable-multifilament |beta.1[16] | -1.0279| -1.9485| -0.1391|
#> 
#> 
#> Model Fit Statistics
#> Deviance = 1016.7
#> Residual deviance = 287.7
#> Deviance Information Criterion (DIC) = 1140.1

The linear dose-response coefficients can then be interpreted as the relative effect for each intervention versus the network reference ("Suture-absorbable").

Benefits of using mbnma.run() for standard NMA models

Now that we have shown how to specify a standard NMA model within the MBNMA framework in mbnma.run(), we can now use MBNMAdose to implement some more interesting models, such as class effect models and node-splits to assess consistency.

A class effects model can be implemented using the class.effect argument in mbnma.run(), introducing either a "common" or 2 class effect on the single linear dose-response parameter, beta.1:

# Random class effect model
nma.class <- mbnma.run(network.plac, fun = dpoly(degree = 1), class.effect = list(beta.1 = "random"),
    n.iter = 50000)

A "common" class effect assumes that all treatments within a class have the same effect, whilst a "random" class effect assumes that treatment-level effects are randomly distributed around a mean class effect with a standard deviation (SD) that is estimated within the model.

Common and random class effect models can be compared using model fit statistics (e.g. Deviance Information Criterion) to identify which is the most parsimonious model. Note that within MBNMAdose when a random class effect is fitted this makes the assumption that all classes share the same within-class SD. This may not necessarily be valid, but relaxing this cannot currently be done in MBNMAdose and it requires specific JAGS code to be written.