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  Treatments                           
#>  AASK   3: Beta Blocker | CCB | ACE Inhibitor
#>  ALLHAT 3: Diuretic | CCB | ACE Inhibitor    
#>  ALPINE 2: Diuretic | ARB                    
#>  ANBP-2 2: Diuretic | ACE Inhibitor          
#>  ASCOT  2: Beta Blocker | CCB                
#>  CAPPP  2: Beta Blocker | ACE Inhibitor      
#>  CHARM  2: Placebo | ARB                     
#>  DREAM  2: Placebo | ACE Inhibitor           
#>  EWPH   2: Diuretic | Placebo                
#>  FEVER  2: Placebo | CCB                     
#>  ... 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  1632
#> d[ARB]               -0.40    0.00 0.05     -0.49     -0.42     -0.40     -0.36     -0.31  2224
#> d[CCB]               -0.20    0.00 0.03     -0.26     -0.22     -0.20     -0.17     -0.13  2045
#> d[Diuretic]           0.06    0.00 0.06     -0.05      0.02      0.06      0.09      0.17  1956
#> d[Placebo]           -0.19    0.00 0.05     -0.29     -0.22     -0.19     -0.16     -0.09  1518
#> lp__             -37970.40    0.08 3.64 -37978.41 -37972.67 -37970.06 -37967.82 -37964.13  2082
#>                  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 18 09:38:31 2022.
#> 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.0 0.08     -0.49     -0.38     -0.33     -0.28     -0.18  1553
#> d[ARB]               -0.40     0.0 0.10     -0.59     -0.46     -0.40     -0.34     -0.21  1722
#> d[CCB]               -0.17     0.0 0.07     -0.30     -0.21     -0.17     -0.13     -0.03  1871
#> d[Diuretic]           0.07     0.0 0.09     -0.10      0.01      0.07      0.13      0.24  1964
#> d[Placebo]           -0.22     0.0 0.09     -0.39     -0.27     -0.21     -0.16     -0.05  1198
#> lp__             -37980.94     0.2 6.69 -37994.55 -37985.36 -37980.74 -37976.31 -37968.54  1146
#> tau                   0.13     0.0 0.04      0.05      0.10      0.12      0.16      0.23   963
#>                  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 18 09:38:58 2022.
#> 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.4 (on 48 data points)
#>                pD: 27.2
#>               DIC: 105.6
(dic_RE <- dic(db_fit_RE))
#> Residual deviance: 53.7 (on 48 data points)
#>                pD: 38.1
#>               DIC: 91.8

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     1972     2882    1
#> d[ACE Inhibitor] -0.36 0.05 -0.46 -0.39 -0.36 -0.32 -0.26     5021     3586    1
#> d[ARB]           -0.45 0.06 -0.57 -0.50 -0.45 -0.41 -0.33     3775     3187    1
#> d[CCB]           -0.25 0.05 -0.36 -0.29 -0.25 -0.22 -0.15     3096     3305    1
#> d[Placebo]       -0.25 0.05 -0.35 -0.28 -0.25 -0.21 -0.14     4678     3402    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.24 -0.13 -0.07 -0.01  0.10     2035     2831    1
#> d[ACE Inhibitor] -0.40 0.09 -0.58 -0.46 -0.40 -0.35 -0.24     4754     2851    1
#> d[ARB]           -0.47 0.11 -0.70 -0.55 -0.47 -0.40 -0.26     4105     2997    1
#> d[CCB]           -0.24 0.08 -0.41 -0.30 -0.24 -0.18 -0.08     4296     2924    1
#> d[Placebo]       -0.29 0.09 -0.47 -0.34 -0.29 -0.22 -0.12     4130     2979    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.07 0.01 0.02 0.04 0.08  0.24     4213     3960    1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.00 0.02 0.03 0.06  0.19     4225     4000    1
#> pred[New 1: ARB]           0.04 0.05 0.00 0.02 0.03 0.05  0.17     4218     3939    1
#> pred[New 1: CCB]           0.05 0.06 0.01 0.02 0.04 0.06  0.20     4217     3960    1
#> pred[New 1: Diuretic]      0.07 0.07 0.01 0.02 0.05 0.08  0.25     4220     3960    1
#> pred[New 1: Placebo]       0.05 0.06 0.01 0.02 0.04 0.07  0.21     4227     3939    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.23     3903     3847    1
#> pred[New 1: ACE Inhibitor] 0.04 0.05 0.00 0.02 0.03 0.05  0.17     3932     4054    1
#> pred[New 1: ARB]           0.04 0.04 0.00 0.01 0.03 0.05  0.16     3902     3941    1
#> pred[New 1: CCB]           0.05 0.05 0.01 0.02 0.04 0.06  0.20     3918     4012    1
#> pred[New 1: Diuretic]      0.06 0.07 0.01 0.02 0.05 0.08  0.25     3912     4055    1
#> pred[New 1: Placebo]       0.05 0.05 0.01 0.02 0.03 0.06  0.19     3937     4011    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     6514     3047    1
#> pred[AASK: ACE Inhibitor] 0.12 0.01 0.10 0.12 0.12 0.13  0.15     3808     3107    1
#> pred[AASK: ARB]           0.12 0.01 0.09 0.11 0.12 0.13  0.15     4206     2961    1
#> pred[AASK: CCB]           0.14 0.01 0.12 0.13 0.14 0.15  0.18     5341     2962    1
#> pred[AASK: Diuretic]      0.18 0.02 0.15 0.17 0.18 0.19  0.22     4187     2876    1
#> pred[AASK: Placebo]       0.14 0.02 0.11 0.13 0.14 0.15  0.17     2888     2946    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.05     3233     2336    1
#> pred[ALLHAT: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03  0.04     4773     2561    1
#> pred[ALLHAT: ARB]           0.03 0.00 0.02 0.03 0.03 0.03  0.04     4297     2854    1
#> pred[ALLHAT: CCB]           0.04 0.00 0.03 0.03 0.04 0.04  0.05     4541     2748    1
#> pred[ALLHAT: Diuretic]      0.05 0.01 0.04 0.04 0.05 0.05  0.06     4782     2671    1
#> pred[ALLHAT: Placebo]       0.03 0.00 0.03 0.03 0.03 0.04  0.04     4618     3130    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.03 0.03  0.05     5806     2580    1
#> pred[ALPINE: ACE Inhibitor] 0.02 0.01 0.01 0.01 0.02 0.02  0.04     6430     2914    1
#> pred[ALPINE: ARB]           0.02 0.01 0.01 0.01 0.02 0.02  0.03     6693     2869    1
#> pred[ALPINE: CCB]           0.02 0.01 0.01 0.02 0.02 0.03  0.04     6297     2911    1
#> pred[ALPINE: Diuretic]      0.03 0.01 0.01 0.02 0.03 0.03  0.05     6564     2683    1
#> pred[ALPINE: Placebo]       0.02 0.01 0.01 0.02 0.02 0.03  0.04     7042     2956    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     2939     2579    1
#> pred[ANBP-2: ACE Inhibitor] 0.05 0.01 0.04 0.04 0.05 0.05  0.06     4760     2727    1
#> pred[ANBP-2: ARB]           0.05 0.01 0.03 0.04 0.05 0.05  0.06     4420     2831    1
#> pred[ANBP-2: CCB]           0.06 0.01 0.04 0.05 0.06 0.06  0.08     4110     2921    1
#> pred[ANBP-2: Diuretic]      0.07 0.01 0.06 0.07 0.07 0.08  0.09     5157     2844    1
#> pred[ANBP-2: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     4557     2978    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     5497     2810    1
#> pred[ASCOT: ACE Inhibitor] 0.08 0.01 0.07 0.08 0.08 0.09  0.10     2000     2675    1
#> pred[ASCOT: ARB]           0.08 0.01 0.06 0.07 0.08 0.08  0.09     2310     2692    1
#> pred[ASCOT: CCB]           0.10 0.01 0.08 0.09 0.09 0.10  0.11     2270     2623    1
#> pred[ASCOT: Diuretic]      0.12 0.01 0.10 0.11 0.12 0.13  0.14     2562     2985    1
#> pred[ASCOT: Placebo]       0.09 0.01 0.08 0.09 0.09 0.10  0.11     1510     2116    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     5812     3130    1
#> pred[CAPPP: ACE Inhibitor] 0.05 0.00 0.05 0.05 0.05 0.06  0.06     1712     2483    1
#> pred[CAPPP: ARB]           0.05 0.01 0.04 0.05 0.05 0.05  0.06     1974     3195    1
#> pred[CAPPP: CCB]           0.06 0.00 0.05 0.06 0.06 0.07  0.07     2841     2743    1
#> pred[CAPPP: Diuretic]      0.08 0.01 0.07 0.08 0.08 0.08  0.10     2418     3151    1
#> pred[CAPPP: Placebo]       0.06 0.01 0.05 0.06 0.06 0.06  0.07     1412     2720    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     2405     2394    1
#> pred[CHARM: ACE Inhibitor] 0.07 0.01 0.05 0.06 0.07 0.07  0.09     4325     2954    1
#> pred[CHARM: ARB]           0.06 0.01 0.05 0.06 0.06 0.07  0.08     4666     2687    1
#> pred[CHARM: CCB]           0.08 0.01 0.06 0.07 0.08 0.08  0.10     3976     2679    1
#> pred[CHARM: Diuretic]      0.10 0.01 0.07 0.09 0.10 0.11  0.13     4300     3084    1
#> pred[CHARM: Placebo]       0.07 0.01 0.06 0.07 0.07 0.08  0.09     4782     2807    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     1947     2479    1
#> pred[DREAM: ACE Inhibitor] 0.17 0.02 0.13 0.16 0.17 0.18  0.21     3383     2472    1
#> pred[DREAM: ARB]           0.16 0.02 0.12 0.15 0.16 0.17  0.21     3493     2668    1
#> pred[DREAM: CCB]           0.20 0.03 0.15 0.18 0.19 0.21  0.25     3078     2401    1
#> pred[DREAM: Diuretic]      0.24 0.03 0.19 0.22 0.24 0.26  0.32     3482     2749    1
#> pred[DREAM: Placebo]       0.19 0.02 0.15 0.17 0.19 0.20  0.23     3745     2456    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     3593     3178    1
#> pred[EWPH: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.04 0.05  0.06     5597     3334    1
#> pred[EWPH: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     5101     2980    1
#> pred[EWPH: CCB]           0.05 0.01 0.04 0.05 0.05 0.06  0.08     4967     3164    1
#> pred[EWPH: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.07  0.09     5467     3015    1
#> pred[EWPH: Placebo]       0.05 0.01 0.03 0.04 0.05 0.06  0.07     5769     3116    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     2744     2207    1
#> pred[FEVER: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03  0.04     4417     2336    1
#> pred[FEVER: ARB]           0.03 0.00 0.02 0.03 0.03 0.03  0.04     4315     2790    1
#> pred[FEVER: CCB]           0.04 0.00 0.03 0.03 0.03 0.04  0.05     4233     2899    1
#> pred[FEVER: Diuretic]      0.04 0.01 0.03 0.04 0.04 0.05  0.06     4331     2750    1
#> pred[FEVER: Placebo]       0.03 0.00 0.03 0.03 0.03 0.04  0.04     4735     2487    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     5975     3214    1
#> pred[HAPPHY: ACE Inhibitor] 0.02  0 0.01 0.02 0.02 0.02  0.02     4446     3287    1
#> pred[HAPPHY: ARB]           0.02  0 0.01 0.02 0.02 0.02  0.02     4150     2840    1
#> pred[HAPPHY: CCB]           0.02  0 0.02 0.02 0.02 0.02  0.03     4710     3405    1
#> pred[HAPPHY: Diuretic]      0.03  0 0.02 0.02 0.03 0.03  0.03     4031     3491    1
#> pred[HAPPHY: Placebo]       0.02  0 0.02 0.02 0.02 0.02  0.02     3306     3129    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     2560     2792    1
#> pred[HOPE: ACE Inhibitor] 0.04 0.01 0.03 0.04 0.04 0.05  0.05     4682     3333    1
#> pred[HOPE: ARB]           0.04 0.01 0.03 0.04 0.04 0.04  0.05     4027     3170    1
#> pred[HOPE: CCB]           0.05 0.01 0.04 0.04 0.05 0.05  0.07     3852     3349    1
#> pred[HOPE: Diuretic]      0.06 0.01 0.05 0.06 0.06 0.07  0.08     4030     2994    1
#> pred[HOPE: Placebo]       0.05 0.01 0.04 0.04 0.05 0.05  0.06     4241     3164    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.09     2718     2434 1.00
#> pred[INSIGHT: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05  0.06     3946     2836 1.00
#> pred[INSIGHT: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     4018     2626 1.00
#> pred[INSIGHT: CCB]           0.06 0.01 0.04 0.05 0.05 0.06  0.07     3693     2578 1.00
#> pred[INSIGHT: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.08  0.09     4330     2302 1.01
#> pred[INSIGHT: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     3566     2972 1.00
#> 
#> ----------------------------------------------------------------- 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     8723     2831    1
#> pred[INVEST: ACE Inhibitor] 0.06 0.00 0.05 0.06 0.06 0.06  0.07     1853     2432    1
#> pred[INVEST: ARB]           0.06 0.01 0.05 0.05 0.06 0.06  0.07     1988     2643    1
#> pred[INVEST: CCB]           0.07 0.00 0.06 0.07 0.07 0.07  0.08     2056     2348    1
#> pred[INVEST: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.09  0.11     2221     3139    1
#> pred[INVEST: Placebo]       0.07 0.01 0.06 0.06 0.07 0.07  0.08     1362     2427    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     7576     2824    1
#> pred[LIFE: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06  0.07     2314     2754    1
#> pred[LIFE: ARB]           0.06 0.01 0.05 0.05 0.06 0.06  0.07     2272     2700    1
#> pred[LIFE: CCB]           0.07 0.01 0.06 0.07 0.07 0.07  0.08     3015     2555    1
#> pred[LIFE: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.09  0.10     2831     2996    1
#> pred[LIFE: Placebo]       0.07 0.01 0.05 0.06 0.07 0.07  0.08     1675     2242    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     4318     2854    1
#> pred[MRC-E: ACE Inhibitor] 0.02  0 0.02 0.02 0.02 0.02  0.03     5802     3286    1
#> pred[MRC-E: ARB]           0.02  0 0.01 0.02 0.02 0.02  0.03     5163     3074    1
#> pred[MRC-E: CCB]           0.03  0 0.02 0.02 0.02 0.03  0.03     5156     3256    1
#> pred[MRC-E: Diuretic]      0.03  0 0.02 0.03 0.03 0.03  0.04     4680     3061    1
#> pred[MRC-E: Placebo]       0.02  0 0.02 0.02 0.02 0.03  0.03     4670     2990    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     6519     3188    1
#> pred[NORDIL: ACE Inhibitor] 0.04 0.00 0.03 0.03 0.04 0.04  0.04     2503     2803    1
#> pred[NORDIL: ARB]           0.03 0.00 0.03 0.03 0.03 0.04  0.04     2366     2782    1
#> pred[NORDIL: CCB]           0.04 0.00 0.04 0.04 0.04 0.04  0.05     2778     2755    1
#> pred[NORDIL: Diuretic]      0.05 0.01 0.04 0.05 0.05 0.06  0.06     2976     2707    1
#> pred[NORDIL: Placebo]       0.04 0.00 0.03 0.04 0.04 0.04  0.05     1713     2504    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.11 0.13 0.14 0.15  0.18     2206     2047    1
#> pred[PEACE: ACE Inhibitor] 0.10 0.01 0.08 0.09 0.10 0.11  0.13     3798     2706    1
#> pred[PEACE: ARB]           0.09 0.01 0.07 0.09 0.09 0.10  0.13     4082     2790    1
#> pred[PEACE: CCB]           0.12 0.02 0.09 0.11 0.12 0.13  0.15     3282     2490    1
#> pred[PEACE: Diuretic]      0.15 0.02 0.11 0.13 0.15 0.16  0.19     3893     3095    1
#> pred[PEACE: Placebo]       0.11 0.01 0.09 0.10 0.11 0.12  0.14     4539     2809    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     2419     2622    1
#> pred[SCOPE: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05  0.06     4179     2864    1
#> pred[SCOPE: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     4579     2749    1
#> pred[SCOPE: CCB]           0.06 0.01 0.04 0.05 0.05 0.06  0.07     3331     2736    1
#> pred[SCOPE: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.08  0.09     4211     3110    1
#> pred[SCOPE: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     4779     3006    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.11     2434     2424    1
#> pred[SHEP: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.07  0.08     3746     2495    1
#> pred[SHEP: ARB]           0.06 0.01 0.04 0.05 0.06 0.06  0.08     4053     2318    1
#> pred[SHEP: CCB]           0.07 0.01 0.05 0.07 0.07 0.08  0.10     3352     2687    1
#> pred[SHEP: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.10  0.12     4148     2811    1
#> pred[SHEP: Placebo]       0.07 0.01 0.05 0.06 0.07 0.08  0.09     3836     2571    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.05 0.05 0.05 0.06  0.06     4917     2939    1
#> pred[STOP-2: ACE Inhibitor] 0.04 0.00 0.03 0.04 0.04 0.04  0.05     3180     3004    1
#> pred[STOP-2: ARB]           0.04 0.00 0.03 0.03 0.04 0.04  0.04     3245     2903    1
#> pred[STOP-2: CCB]           0.05 0.00 0.04 0.04 0.05 0.05  0.05     3784     2183    1
#> pred[STOP-2: Diuretic]      0.06 0.01 0.05 0.05 0.06 0.06  0.07     4117     3361    1
#> pred[STOP-2: Placebo]       0.04 0.00 0.03 0.04 0.04 0.05  0.05     2642     2766    1
#> 
#> ------------------------------------------------------------------ Study: VALUE ---- 
#> 
#>                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[VALUE: Beta Blocker]  0.20 0.02 0.15 0.18 0.19 0.21  0.25     2951     2400    1
#> pred[VALUE: ACE Inhibitor] 0.15 0.02 0.11 0.13 0.14 0.16  0.19     3881     2484    1
#> pred[VALUE: ARB]           0.14 0.02 0.10 0.13 0.14 0.15  0.17     4142     2794    1
#> pred[VALUE: CCB]           0.17 0.02 0.13 0.16 0.17 0.18  0.21     4284     2685    1
#> pred[VALUE: Diuretic]      0.21 0.03 0.16 0.19 0.21 0.22  0.27     4070     2816    1
#> pred[VALUE: Placebo]       0.16 0.02 0.12 0.15 0.16 0.17  0.21     4139     2642    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     2421       NA    1
#> rank[ACE Inhibitor] 1.84 0.55    1   2   2   2     3     3752     3569    1
#> rank[ARB]           1.28 0.53    1   1   1   1     3     3603     3312    1
#> rank[CCB]           3.70 0.53    3   3   4   4     4     3538     2767    1
#> rank[Diuretic]      5.80 0.41    5   6   6   6     6     2878       NA    1
#> rank[Placebo]       3.21 0.60    2   3   3   4     4     2714     2959    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.01      0.79       0.2
#> d[ACE Inhibitor]      0.23      0.70      0.06      0.01      0.00       0.0
#> d[ARB]                0.76      0.21      0.03      0.00      0.00       0.0
#> d[CCB]                0.00      0.02      0.28      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.63      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.76      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. http://nicedsu.org.uk/.
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.↩︎