Example: Diabetes

library(multinma)
options(mc.cores = parallel::detectCores())

This vignette describes the analysis of data on the number of new cases of diabetes in 22 trials of 6 antihypertensive drugs (Elliott and Meyer 2007; Dias et al. 2011). The data are available in this package as diabetes:

head(diabetes)
#>   studyn studyc trtn         trtc   r    n time
#> 1      1  MRC-E    1     Diuretic  43 1081  5.8
#> 2      1  MRC-E    2      Placebo  34 2213  5.8
#> 3      1  MRC-E    3 Beta Blocker  37 1102  5.8
#> 4      2   EWPH    1     Diuretic  29  416  4.7
#> 5      2   EWPH    2      Placebo  20  424  4.7
#> 6      3   SHEP    1     Diuretic 140 1631  3.0

Setting up the network

We begin by setting up the network. We have arm-level count data giving the number of new cases of diabetes (r) out of the total (n) in each arm, so we use the function set_agd_arm(). For computational efficiency, we let “Beta Blocker” be set as the network reference treatment by default. Elliott and Meyer (2007) and Dias et al. (2011) use “Diuretic” as the reference, but it is a simple matter to transform the results after fitting the NMA model.1

db_net <- set_agd_arm(diabetes, 
                      study = studyc,
                      trt = trtc,
                      r = r, 
                      n = n)
db_net
#> A network with 22 AgD studies (arm-based).
#> 
#> ------------------------------------------------------- AgD studies (arm-based) ---- 
#>  Study  Treatment arms                       
#>  AASK   3: Beta Blocker | ACE Inhibitor | CCB
#>  ALLHAT 3: ACE Inhibitor | CCB | Diuretic    
#>  ALPINE 2: ARB | Diuretic                    
#>  ANBP-2 2: ACE Inhibitor | Diuretic          
#>  ASCOT  2: Beta Blocker | CCB                
#>  CAPPP  2: Beta Blocker | ACE Inhibitor      
#>  CHARM  2: ARB | Placebo                     
#>  DREAM  2: ACE Inhibitor | Placebo           
#>  EWPH   2: Diuretic | Placebo                
#>  FEVER  2: CCB | Placebo                     
#>  ... plus 12 more studies
#> 
#>  Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 6
#> Total number of studies: 22
#> Reference treatment is: Beta Blocker
#> Network is connected

We also have details of length of follow-up in years in each trial (time), which we will use as an offset with a cloglog link function to model the data as rates. We do not have to specify this in the function set_agd_arm(): any additional columns in the data (e.g. offsets or covariates, here the column time) will automatically be made available in the network.

Plot the network structure.

plot(db_net, weight_edges = TRUE, weight_nodes = TRUE)

Meta-analysis models

We fit both fixed effect (FE) and random effects (RE) models.

Fixed effect meta-analysis

First, we fit a fixed effect model using the nma() function with trt_effects = "fixed". We use \(\mathrm{N}(0, 100^2)\) prior distributions for the treatment effects \(d_k\) and study-specific intercepts \(\mu_j\). We can examine the range of parameter values implied by these prior distributions with the summary() method:

summary(normal(scale = 100))
#> A Normal prior distribution: location = 0, scale = 100.
#> 50% of the prior density lies between -67.45 and 67.45.
#> 95% of the prior density lies between -196 and 196.

The model is fitted using the nma() function. We specify that a cloglog link will be used with link = "cloglog" (the Binomial likelihood is the default for these data), and specify the log follow-up time offset using the regression formula regression = ~offset(log(time)).

db_fit_FE <- nma(db_net, 
                 trt_effects = "fixed",
                 link = "cloglog",
                 regression = ~offset(log(time)),
                 prior_intercept = normal(scale = 100),
                 prior_trt = normal(scale = 100))
#> Note: No treatment classes specified in network, any interactions in `regression` formula will be separate (independent) for each treatment.
#> Use set_*() argument `trt_class` and nma() argument `class_interactions` to change this.
#> Note: Setting "Beta Blocker" as the network reference treatment.

Basic parameter summaries are given by the print() method:

db_fit_FE
#> A fixed effects NMA with a binomial likelihood (cloglog link).
#> Regression model: ~offset(log(time)).
#> Inference for Stan model: binomial_1par.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
#> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> 
#>                       mean se_mean   sd      2.5%       25%       50%       75%     97.5% n_eff
#> d[ACE Inhibitor]     -0.30    0.00 0.05     -0.39     -0.33     -0.30     -0.27     -0.21  1227
#> d[ARB]               -0.39    0.00 0.05     -0.49     -0.43     -0.39     -0.36     -0.30  2010
#> d[CCB]               -0.20    0.00 0.03     -0.26     -0.22     -0.20     -0.17     -0.13  1671
#> d[Diuretic]           0.06    0.00 0.06     -0.05      0.02      0.06      0.09      0.17  1523
#> d[Placebo]           -0.19    0.00 0.05     -0.29     -0.22     -0.19     -0.16     -0.09  1179
#> lp__             -38119.60    0.09 3.70 -38127.60 -38121.85 -38119.25 -38117.00 -38113.27  1748
#>                  Rhat
#> d[ACE Inhibitor]    1
#> d[ARB]              1
#> d[CCB]              1
#> d[Diuretic]         1
#> d[Placebo]          1
#> lp__                1
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Jan  9 17:56:28 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).

By default, summaries of the study-specific intercepts \(\mu_j\) are hidden, but could be examined by changing the pars argument:

# Not run
print(db_fit_FE, pars = c("d", "mu"))

The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:

plot_prior_posterior(db_fit_FE)

Random effects meta-analysis

We now fit a random effects model using the nma() function with trt_effects = "random". Again, we use \(\mathrm{N}(0, 100^2)\) prior distributions for the treatment effects \(d_k\) and study-specific intercepts \(\mu_j\), and we additionally use a \(\textrm{half-N}(5^2)\) prior for the heterogeneity standard deviation \(\tau\). We can examine the range of parameter values implied by these prior distributions with the summary() method:

summary(normal(scale = 100))
#> A Normal prior distribution: location = 0, scale = 100.
#> 50% of the prior density lies between -67.45 and 67.45.
#> 95% of the prior density lies between -196 and 196.
summary(half_normal(scale = 5))
#> A half-Normal prior distribution: location = 0, scale = 5.
#> 50% of the prior density lies between 0 and 3.37.
#> 95% of the prior density lies between 0 and 9.8.

Fitting the RE model

db_fit_RE <- nma(db_net, 
                 trt_effects = "random",
                 link = "cloglog",
                 regression = ~offset(log(time)),
                 prior_intercept = normal(scale = 10),
                 prior_trt = normal(scale = 10),
                 prior_het = half_normal(scale = 5),
                 init_r = 0.5)
#> Note: No treatment classes specified in network, any interactions in `regression` formula will be separate (independent) for each treatment.
#> Use set_*() argument `trt_class` and nma() argument `class_interactions` to change this.
#> Note: Setting "Beta Blocker" as the network reference treatment.

Basic parameter summaries are given by the print() method:

db_fit_RE
#> A random effects NMA with a binomial likelihood (cloglog link).
#> Regression model: ~offset(log(time)).
#> Inference for Stan model: binomial_1par.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
#> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> 
#>                       mean se_mean   sd      2.5%       25%       50%       75%     97.5% n_eff
#> d[ACE Inhibitor]     -0.33    0.00 0.08     -0.49     -0.38     -0.33     -0.28     -0.18  1438
#> d[ARB]               -0.40    0.00 0.10     -0.59     -0.46     -0.40     -0.34     -0.21  1606
#> d[CCB]               -0.17    0.00 0.07     -0.30     -0.21     -0.17     -0.13     -0.03  1802
#> d[Diuretic]           0.07    0.00 0.09     -0.10      0.01      0.07      0.13      0.25  1737
#> d[Placebo]           -0.21    0.00 0.09     -0.40     -0.27     -0.21     -0.16     -0.04  1237
#> lp__             -38070.21    0.22 6.87 -38084.72 -38074.58 -38069.88 -38065.62 -38057.55   967
#> tau                   0.13    0.00 0.05      0.06      0.10      0.13      0.16      0.24   845
#>                  Rhat
#> d[ACE Inhibitor] 1.00
#> d[ARB]           1.00
#> d[CCB]           1.00
#> d[Diuretic]      1.00
#> d[Placebo]       1.00
#> lp__             1.00
#> tau              1.01
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Jan  9 17:56:51 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).

By default, summaries of the study-specific intercepts \(\mu_j\) and study-specific relative effects \(\delta_{jk}\) are hidden, but could be examined by changing the pars argument:

# Not run
print(db_fit_RE, pars = c("d", "mu", "delta"))

The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:

plot_prior_posterior(db_fit_RE, prior = c("trt", "het"))

Model comparison

Model fit can be checked using the dic() function:

(dic_FE <- dic(db_fit_FE))
#> Residual deviance: 78.5 (on 48 data points)
#>                pD: 27.3
#>               DIC: 105.8
(dic_RE <- dic(db_fit_RE))
#> Residual deviance: 53.4 (on 48 data points)
#>                pD: 38.1
#>               DIC: 91.5

The FE model is a very poor fit to the data, with a residual deviance much higher than the number of data points. The RE model fits the data better, and has a much lower DIC; we prefer the RE model.

We can also examine the residual deviance contributions with the corresponding plot() method.

plot(dic_FE)

plot(dic_RE)

Further results

For comparison with Elliott and Meyer (2007) and Dias et al. (2011), we can produce relative effects against “Diuretic” using the relative_effects() function with trt_ref = "Diuretic":

(db_releff_FE <- relative_effects(db_fit_FE, trt_ref = "Diuretic"))
#>                   mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Beta Blocker]  -0.06 0.06 -0.17 -0.09 -0.06 -0.02  0.05     1534     2566    1
#> d[ACE Inhibitor] -0.36 0.05 -0.47 -0.39 -0.36 -0.32 -0.25     4582     3450    1
#> d[ARB]           -0.45 0.06 -0.58 -0.49 -0.45 -0.41 -0.33     3181     3052    1
#> d[CCB]           -0.25 0.05 -0.36 -0.29 -0.25 -0.22 -0.15     2595     3234    1
#> d[Placebo]       -0.25 0.06 -0.36 -0.28 -0.25 -0.21 -0.14     3992     2883    1
plot(db_releff_FE, ref_line = 0)

(db_releff_RE <- relative_effects(db_fit_RE, trt_ref = "Diuretic"))
#>                   mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Beta Blocker]  -0.07 0.09 -0.25 -0.13 -0.07 -0.01  0.10     1776     2284    1
#> d[ACE Inhibitor] -0.40 0.09 -0.58 -0.46 -0.40 -0.34 -0.24     4083     3200    1
#> d[ARB]           -0.47 0.11 -0.69 -0.54 -0.47 -0.40 -0.26     3293     2854    1
#> d[CCB]           -0.24 0.09 -0.41 -0.30 -0.24 -0.18 -0.07     3577     3316    1
#> d[Placebo]       -0.29 0.09 -0.47 -0.34 -0.28 -0.23 -0.12     3835     2964    1
plot(db_releff_RE, ref_line = 0)

Dias et al. (2011) produce absolute predictions of the probability of developing diabetes after three years, assuming a Normal distribution on the baseline cloglog probability of developing diabetes on diuretic treatment with mean \(-4.2\) and precision \(1.11\). We can replicate these results using the predict() method. We specify a data frame of newdata, containing the time offset(s) at which to produce predictions (here only 3 years). The baseline argument takes a distr() distribution object with which we specify the corresponding Normal distribution on the baseline cloglog probability, and we set trt_ref = "Diuretic" to indicate that the baseline distribution corresponds to “Diuretic” rather than the network reference “Beta Blocker”. We set type = "response" to produce predicted event probabilities (type = "link" would produce predicted cloglog probabilities).

db_pred_FE <- predict(db_fit_FE, 
                      newdata = data.frame(time = 3),
                      baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5), 
                      trt_ref = "Diuretic",
                      type = "response")
db_pred_FE
#> ------------------------------------------------------------------ Study: New 1 ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[New 1: Beta Blocker]  0.06 0.06 0.01 0.02 0.04 0.07  0.23     3871     3891    1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.00 0.02 0.03 0.06  0.18     3895     3852    1
#> pred[New 1: ARB]           0.04 0.04 0.00 0.01 0.03 0.05  0.16     3899     3889    1
#> pred[New 1: CCB]           0.05 0.05 0.01 0.02 0.03 0.06  0.20     3887     3852    1
#> pred[New 1: Diuretic]      0.06 0.06 0.01 0.02 0.04 0.08  0.24     3905     3891    1
#> pred[New 1: Placebo]       0.05 0.05 0.01 0.02 0.03 0.06  0.20     3904     3813    1
plot(db_pred_FE)

db_pred_RE <- predict(db_fit_RE, 
                      newdata = data.frame(time = 3),
                      baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5), 
                      trt_ref = "Diuretic",
                      type = "response")
db_pred_RE
#> ------------------------------------------------------------------ Study: New 1 ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[New 1: Beta Blocker]  0.06 0.06 0.01 0.02 0.04 0.08  0.25     3784     3932    1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.00 0.02 0.03 0.06  0.18     3789     3853    1
#> pred[New 1: ARB]           0.04 0.05 0.00 0.01 0.03 0.05  0.17     3784     3892    1
#> pred[New 1: CCB]           0.05 0.06 0.01 0.02 0.03 0.06  0.21     3779     3930    1
#> pred[New 1: Diuretic]      0.07 0.07 0.01 0.02 0.04 0.08  0.26     3776     3852    1
#> pred[New 1: Placebo]       0.05 0.05 0.01 0.02 0.03 0.06  0.21     3772     3852    1
plot(db_pred_RE)

If the baseline and newdata arguments are omitted, predicted probabilities will be produced for every study in the network based on their follow-up times and estimated baseline cloglog probabilities \(\mu_j\):

db_pred_RE_studies <- predict(db_fit_RE, type = "response")
db_pred_RE_studies
#> ------------------------------------------------------------------- Study: AASK ---- 
#> 
#>                           mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[AASK: Beta Blocker]  0.17 0.02 0.14 0.16 0.17 0.18  0.20     4723     3246    1
#> pred[AASK: ACE Inhibitor] 0.13 0.01 0.10 0.12 0.12 0.13  0.15     3772     3061    1
#> pred[AASK: ARB]           0.12 0.01 0.09 0.11 0.12 0.13  0.15     3487     2967    1
#> pred[AASK: CCB]           0.15 0.02 0.12 0.14 0.14 0.16  0.18     4204     3078    1
#> pred[AASK: Diuretic]      0.18 0.02 0.14 0.17 0.18 0.19  0.22     3476     3153    1
#> pred[AASK: Placebo]       0.14 0.02 0.11 0.13 0.14 0.15  0.17     2844     3267    1
#> 
#> ----------------------------------------------------------------- Study: ALLHAT ---- 
#> 
#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ALLHAT: Beta Blocker]  0.04 0.01 0.03 0.04 0.04 0.05  0.06     2494     2348    1
#> pred[ALLHAT: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03  0.04     3888     2440    1
#> pred[ALLHAT: ARB]           0.03 0.00 0.02 0.03 0.03 0.03  0.04     3598     2730    1
#> pred[ALLHAT: CCB]           0.04 0.00 0.03 0.03 0.04 0.04  0.05     3558     2108    1
#> pred[ALLHAT: Diuretic]      0.05 0.01 0.04 0.04 0.05 0.05  0.06     3860     2665    1
#> pred[ALLHAT: Placebo]       0.03 0.00 0.03 0.03 0.03 0.04  0.05     3687     2481    1
#> 
#> ----------------------------------------------------------------- Study: ALPINE ---- 
#> 
#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ALPINE: Beta Blocker]  0.03 0.01 0.01 0.02 0.02 0.03  0.05     5764     3098    1
#> pred[ALPINE: ACE Inhibitor] 0.02 0.01 0.01 0.01 0.02 0.02  0.03     6749     3056    1
#> pred[ALPINE: ARB]           0.02 0.01 0.01 0.01 0.02 0.02  0.03     7128     2936    1
#> pred[ALPINE: CCB]           0.02 0.01 0.01 0.02 0.02 0.03  0.04     6932     3156    1
#> pred[ALPINE: Diuretic]      0.03 0.01 0.01 0.02 0.03 0.03  0.05     7259     2950    1
#> pred[ALPINE: Placebo]       0.02 0.01 0.01 0.02 0.02 0.03  0.04     6941     2811    1
#> 
#> ----------------------------------------------------------------- Study: ANBP-2 ---- 
#> 
#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ANBP-2: Beta Blocker]  0.07 0.01 0.05 0.06 0.07 0.07  0.09     2612     2178    1
#> pred[ANBP-2: ACE Inhibitor] 0.05 0.01 0.04 0.04 0.05 0.05  0.06     4351     2710    1
#> pred[ANBP-2: ARB]           0.05 0.01 0.03 0.04 0.05 0.05  0.06     3664     2424    1
#> pred[ANBP-2: CCB]           0.06 0.01 0.04 0.05 0.06 0.06  0.08     3639     2552    1
#> pred[ANBP-2: Diuretic]      0.07 0.01 0.06 0.07 0.07 0.08  0.09     4509     2419    1
#> pred[ANBP-2: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     4125     2600    1
#> 
#> ------------------------------------------------------------------ Study: ASCOT ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ASCOT: Beta Blocker]  0.11 0.00 0.10 0.11 0.11 0.11  0.12     5041     2398    1
#> pred[ASCOT: ACE Inhibitor] 0.08 0.01 0.07 0.08 0.08 0.09  0.10     2056     2404    1
#> pred[ASCOT: ARB]           0.08 0.01 0.06 0.07 0.08 0.08  0.09     1963     2425    1
#> pred[ASCOT: CCB]           0.10 0.01 0.08 0.09 0.10 0.10  0.11     2089     2170    1
#> pred[ASCOT: Diuretic]      0.12 0.01 0.10 0.11 0.12 0.13  0.14     2014     2428    1
#> pred[ASCOT: Placebo]       0.09 0.01 0.08 0.09 0.09 0.10  0.11     1512     2075    1
#> 
#> ------------------------------------------------------------------ Study: CAPPP ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[CAPPP: Beta Blocker]  0.07 0.00 0.07 0.07 0.07 0.08  0.08     4890     3103    1
#> pred[CAPPP: ACE Inhibitor] 0.05 0.00 0.05 0.05 0.05 0.06  0.06     1680     2220    1
#> pred[CAPPP: ARB]           0.05 0.01 0.04 0.05 0.05 0.05  0.06     2092     2140    1
#> pred[CAPPP: CCB]           0.06 0.00 0.05 0.06 0.06 0.07  0.07     2585     2575    1
#> pred[CAPPP: Diuretic]      0.08 0.01 0.07 0.08 0.08 0.09  0.10     2135     2501    1
#> pred[CAPPP: Placebo]       0.06 0.01 0.05 0.06 0.06 0.06  0.07     1449     1855    1
#> 
#> ------------------------------------------------------------------ Study: CHARM ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[CHARM: Beta Blocker]  0.09 0.01 0.07 0.08 0.09 0.10  0.12     2361     2250    1
#> pred[CHARM: ACE Inhibitor] 0.07 0.01 0.05 0.06 0.07 0.07  0.09     4092     2726    1
#> pred[CHARM: ARB]           0.06 0.01 0.05 0.06 0.06 0.07  0.08     4581     3002    1
#> pred[CHARM: CCB]           0.08 0.01 0.06 0.07 0.08 0.08  0.10     3558     2717    1
#> pred[CHARM: Diuretic]      0.10 0.01 0.07 0.09 0.10 0.10  0.13     4112     2808    1
#> pred[CHARM: Placebo]       0.07 0.01 0.06 0.07 0.07 0.08  0.10     4625     2875    1
#> 
#> ------------------------------------------------------------------ Study: DREAM ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[DREAM: Beta Blocker]  0.23 0.03 0.18 0.21 0.23 0.24  0.29     2492     2621    1
#> pred[DREAM: ACE Inhibitor] 0.17 0.02 0.13 0.16 0.17 0.18  0.21     4212     3075    1
#> pred[DREAM: ARB]           0.16 0.02 0.12 0.14 0.16 0.17  0.21     3712     2893    1
#> pred[DREAM: CCB]           0.20 0.03 0.15 0.18 0.19 0.21  0.25     3651     2834    1
#> pred[DREAM: Diuretic]      0.24 0.03 0.19 0.22 0.24 0.26  0.31     4067     2569    1
#> pred[DREAM: Placebo]       0.19 0.02 0.15 0.17 0.19 0.20  0.24     4509     2798    1
#> 
#> ------------------------------------------------------------------- Study: EWPH ---- 
#> 
#>                           mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[EWPH: Beta Blocker]  0.06 0.01 0.04 0.05 0.06 0.07  0.09     3506     2967    1
#> pred[EWPH: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.04 0.05  0.07     4771     3198    1
#> pred[EWPH: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     4052     3280    1
#> pred[EWPH: CCB]           0.05 0.01 0.04 0.05 0.05 0.06  0.08     4108     3243    1
#> pred[EWPH: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.07  0.10     4852     3138    1
#> pred[EWPH: Placebo]       0.05 0.01 0.03 0.04 0.05 0.06  0.07     4749     2947    1
#> 
#> ------------------------------------------------------------------ Study: FEVER ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[FEVER: Beta Blocker]  0.04 0.01 0.03 0.04 0.04 0.04  0.05     2565     2642    1
#> pred[FEVER: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03  0.04     3706     2837    1
#> pred[FEVER: ARB]           0.03 0.00 0.02 0.02 0.03 0.03  0.04     3737     2459    1
#> pred[FEVER: CCB]           0.03 0.00 0.03 0.03 0.03 0.04  0.05     3697     2988    1
#> pred[FEVER: Diuretic]      0.04 0.01 0.03 0.04 0.04 0.05  0.06     3942     2928    1
#> pred[FEVER: Placebo]       0.03 0.00 0.02 0.03 0.03 0.04  0.04     3795     2544    1
#> 
#> ----------------------------------------------------------------- Study: HAPPHY ---- 
#> 
#>                             mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[HAPPHY: Beta Blocker]  0.02  0 0.02 0.02 0.02 0.03  0.03     4877     3398    1
#> pred[HAPPHY: ACE Inhibitor] 0.02  0 0.01 0.02 0.02 0.02  0.02     3655     2613    1
#> pred[HAPPHY: ARB]           0.02  0 0.01 0.02 0.02 0.02  0.02     3736     3039    1
#> pred[HAPPHY: CCB]           0.02  0 0.02 0.02 0.02 0.02  0.03     3830     2783    1
#> pred[HAPPHY: Diuretic]      0.03  0 0.02 0.02 0.03 0.03  0.03     3133     2826    1
#> pred[HAPPHY: Placebo]       0.02  0 0.02 0.02 0.02 0.02  0.02     2870     3062    1
#> 
#> ------------------------------------------------------------------- Study: HOPE ---- 
#> 
#>                           mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[HOPE: Beta Blocker]  0.06 0.01 0.04 0.05 0.06 0.06  0.08     2246     2122    1
#> pred[HOPE: ACE Inhibitor] 0.04 0.01 0.03 0.04 0.04 0.05  0.06     3893     3079    1
#> pred[HOPE: ARB]           0.04 0.01 0.03 0.04 0.04 0.04  0.05     3666     2986    1
#> pred[HOPE: CCB]           0.05 0.01 0.04 0.04 0.05 0.05  0.07     3231     2836    1
#> pred[HOPE: Diuretic]      0.06 0.01 0.05 0.06 0.06 0.07  0.08     4129     2995    1
#> pred[HOPE: Placebo]       0.05 0.01 0.04 0.04 0.05 0.05  0.06     4106     2752    1
#> 
#> ---------------------------------------------------------------- Study: INSIGHT ---- 
#> 
#>                              mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[INSIGHT: Beta Blocker]  0.07 0.01 0.05 0.06 0.06 0.07  0.08     2646     2628    1
#> pred[INSIGHT: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05  0.06     3590     2713    1
#> pred[INSIGHT: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     3423     2317    1
#> pred[INSIGHT: CCB]           0.06 0.01 0.04 0.05 0.05 0.06  0.07     3740     3088    1
#> pred[INSIGHT: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.07  0.09     4217     2346    1
#> pred[INSIGHT: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     3560     2604    1
#> 
#> ----------------------------------------------------------------- Study: INVEST ---- 
#> 
#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[INVEST: Beta Blocker]  0.08 0.00 0.08 0.08 0.08 0.08  0.09     6040     2876    1
#> pred[INVEST: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06  0.07     2019     2434    1
#> pred[INVEST: ARB]           0.06 0.01 0.05 0.05 0.06 0.06  0.07     1978     2354    1
#> pred[INVEST: CCB]           0.07 0.00 0.06 0.07 0.07 0.07  0.08     2229     2585    1
#> pred[INVEST: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.09  0.11     2069     2556    1
#> pred[INVEST: Placebo]       0.07 0.01 0.06 0.06 0.07 0.07  0.08     1500     1955    1
#> 
#> ------------------------------------------------------------------- Study: LIFE ---- 
#> 
#>                           mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[LIFE: Beta Blocker]  0.08 0.00 0.07 0.08 0.08 0.08  0.09     7078     2746    1
#> pred[LIFE: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06  0.07     2268     2786    1
#> pred[LIFE: ARB]           0.06 0.01 0.05 0.05 0.06 0.06  0.07     1975     2442    1
#> pred[LIFE: CCB]           0.07 0.01 0.06 0.07 0.07 0.07  0.08     2941     2434    1
#> pred[LIFE: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.09  0.11     2351     2509    1
#> pred[LIFE: Placebo]       0.07 0.01 0.05 0.06 0.07 0.07  0.08     1728     2165    1
#> 
#> ------------------------------------------------------------------ Study: MRC-E ---- 
#> 
#>                            mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[MRC-E: Beta Blocker]  0.03  0 0.02 0.03 0.03 0.03  0.04     3514     3056    1
#> pred[MRC-E: ACE Inhibitor] 0.02  0 0.02 0.02 0.02 0.02  0.03     4767     3028    1
#> pred[MRC-E: ARB]           0.02  0 0.01 0.02 0.02 0.02  0.03     4670     3162    1
#> pred[MRC-E: CCB]           0.03  0 0.02 0.02 0.02 0.03  0.03     4412     2808    1
#> pred[MRC-E: Diuretic]      0.03  0 0.02 0.03 0.03 0.03  0.04     4325     2982    1
#> pred[MRC-E: Placebo]       0.02  0 0.02 0.02 0.02 0.03  0.03     4216     3534    1
#> 
#> ----------------------------------------------------------------- Study: NORDIL ---- 
#> 
#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[NORDIL: Beta Blocker]  0.05 0.00 0.04 0.05 0.05 0.05  0.06     6403     3275    1
#> pred[NORDIL: ACE Inhibitor] 0.04 0.00 0.03 0.03 0.04 0.04  0.04     2326     2387    1
#> pred[NORDIL: ARB]           0.03 0.00 0.03 0.03 0.03 0.04  0.04     2037     2491    1
#> pred[NORDIL: CCB]           0.04 0.00 0.04 0.04 0.04 0.04  0.05     2673     2639    1
#> pred[NORDIL: Diuretic]      0.05 0.01 0.04 0.05 0.05 0.06  0.06     2339     2475    1
#> pred[NORDIL: Placebo]       0.04 0.00 0.03 0.04 0.04 0.04  0.05     1833     2091    1
#> 
#> ------------------------------------------------------------------ Study: PEACE ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[PEACE: Beta Blocker]  0.14 0.02 0.10 0.13 0.14 0.15  0.18     2153     2171    1
#> pred[PEACE: ACE Inhibitor] 0.10 0.01 0.08 0.09 0.10 0.11  0.13     3960     2836    1
#> pred[PEACE: ARB]           0.09 0.01 0.07 0.09 0.09 0.10  0.13     3461     2694    1
#> pred[PEACE: CCB]           0.12 0.02 0.09 0.11 0.12 0.13  0.15     3122     2631    1
#> pred[PEACE: Diuretic]      0.15 0.02 0.11 0.13 0.15 0.16  0.19     3703     2711    1
#> pred[PEACE: Placebo]       0.11 0.01 0.09 0.10 0.11 0.12  0.14     4168     2706    1
#> 
#> ------------------------------------------------------------------ Study: SCOPE ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[SCOPE: Beta Blocker]  0.06 0.01 0.05 0.06 0.06 0.07  0.09     2151     2008    1
#> pred[SCOPE: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05  0.06     3607     2638    1
#> pred[SCOPE: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     3985     2702    1
#> pred[SCOPE: CCB]           0.06 0.01 0.04 0.05 0.05 0.06  0.07     2959     2261    1
#> pred[SCOPE: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.08  0.09     3938     2644    1
#> pred[SCOPE: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     4557     2312    1
#> 
#> ------------------------------------------------------------------- Study: SHEP ---- 
#> 
#>                           mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[SHEP: Beta Blocker]  0.09 0.01 0.06 0.08 0.09 0.09  0.12     2162     2221    1
#> pred[SHEP: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.07  0.08     4209     2968    1
#> pred[SHEP: ARB]           0.06 0.01 0.04 0.05 0.06 0.06  0.08     3542     2638    1
#> pred[SHEP: CCB]           0.07 0.01 0.05 0.07 0.07 0.08  0.10     3452     2651    1
#> pred[SHEP: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.10  0.12     4340     2571    1
#> pred[SHEP: Placebo]       0.07 0.01 0.05 0.06 0.07 0.08  0.09     4156     2845    1
#> 
#> ----------------------------------------------------------------- Study: STOP-2 ---- 
#> 
#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[STOP-2: Beta Blocker]  0.05 0.00 0.04 0.05 0.05 0.06  0.06     3564     2894    1
#> pred[STOP-2: ACE Inhibitor] 0.04 0.00 0.03 0.04 0.04 0.04  0.05     2323     2639    1
#> pred[STOP-2: ARB]           0.04 0.00 0.03 0.03 0.04 0.04  0.04     2316     2854    1
#> pred[STOP-2: CCB]           0.05 0.00 0.04 0.04 0.05 0.05  0.05     3358     3095    1
#> pred[STOP-2: Diuretic]      0.06 0.01 0.05 0.05 0.06 0.06  0.07     2722     2768    1
#> pred[STOP-2: Placebo]       0.04 0.00 0.03 0.04 0.04 0.05  0.05     1929     1935    1
#> 
#> ------------------------------------------------------------------ Study: VALUE ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[VALUE: Beta Blocker]  0.20 0.03 0.15 0.18 0.19 0.21  0.25     2623     2133    1
#> pred[VALUE: ACE Inhibitor] 0.15 0.02 0.11 0.13 0.14 0.16  0.19     3370     2479    1
#> pred[VALUE: ARB]           0.14 0.02 0.10 0.13 0.14 0.15  0.17     3700     2326    1
#> pred[VALUE: CCB]           0.17 0.02 0.13 0.16 0.17 0.18  0.21     3514     2400    1
#> pred[VALUE: Diuretic]      0.21 0.03 0.16 0.19 0.21 0.22  0.27     3360     2698    1
#> pred[VALUE: Placebo]       0.16 0.02 0.12 0.15 0.16 0.17  0.21     3284     2180    1
plot(db_pred_RE_studies)

We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.

(db_ranks <- posterior_ranks(db_fit_RE))
#>                     mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[Beta Blocker]  5.18 0.43    5   5   5   5     6     1971       NA    1
#> rank[ACE Inhibitor] 1.85 0.54    1   2   2   2     3     3345     2913    1
#> rank[ARB]           1.27 0.52    1   1   1   1     3     3185     3079    1
#> rank[CCB]           3.69 0.54    3   3   4   4     4     2977     2623    1
#> rank[Diuretic]      5.80 0.41    5   6   6   6     6     2280       NA    1
#> rank[Placebo]       3.21 0.60    2   3   3   4     4     2826     2278    1
plot(db_ranks)

(db_rankprobs <- posterior_rank_probs(db_fit_RE))
#>                  p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
#> d[Beta Blocker]       0.00      0.00      0.00      0.02      0.78       0.2
#> d[ACE Inhibitor]      0.23      0.71      0.06      0.01      0.00       0.0
#> d[ARB]                0.77      0.20      0.03      0.00      0.00       0.0
#> d[CCB]                0.00      0.02      0.27      0.69      0.01       0.0
#> d[Diuretic]           0.00      0.00      0.00      0.00      0.20       0.8
#> d[Placebo]            0.01      0.07      0.64      0.28      0.01       0.0
plot(db_rankprobs)

(db_cumrankprobs <- posterior_rank_probs(db_fit_RE, cumulative = TRUE))
#>                  p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
#> d[Beta Blocker]       0.00      0.00      0.00      0.02       0.8         1
#> d[ACE Inhibitor]      0.23      0.93      0.99      1.00       1.0         1
#> d[ARB]                0.77      0.97      1.00      1.00       1.0         1
#> d[CCB]                0.00      0.02      0.30      0.99       1.0         1
#> d[Diuretic]           0.00      0.00      0.00      0.00       0.2         1
#> d[Placebo]            0.01      0.08      0.71      0.99       1.0         1
plot(db_cumrankprobs)

References

Dias, S., N. J. Welton, A. J. Sutton, and A. E. Ades. 2011. “NICE DSU Technical Support Document 2: A Generalised Linear Modelling Framework for Pair-Wise and Network Meta-Analysis of Randomised Controlled Trials.” National Institute for Health and Care Excellence. https://www.sheffield.ac.uk/nice-dsu.
Elliott, W. J., and P. M. Meyer. 2007. “Incident Diabetes in Clinical Trials of Antihypertensive Drugs: A Network Meta-Analysis.” The Lancet 369 (9557): 201–7. https://doi.org/10.1016/s0140-6736(07)60108-1.

  1. The gain in efficiency here from using “Beta Blocker” as the network reference treatment instead of “Diuretic” is considerable - around 4-8 times, in terms of effective samples per second. The functions in this package will always attempt to choose a default network reference treatment that maximises computational efficiency and stability. If you have chosen an alternative network reference treatment and the model runs very slowly or has low effective sample size, this is a likely cause.↩︎