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.04     -0.39     -0.33     -0.30     -0.27     -0.21  1203
#> d[ARB]               -0.40    0.00 0.05     -0.49     -0.43     -0.40     -0.36     -0.31  2150
#> d[CCB]               -0.20    0.00 0.03     -0.26     -0.22     -0.20     -0.17     -0.14  1807
#> d[Diuretic]           0.06    0.00 0.06     -0.06      0.02      0.06      0.10      0.17  1447
#> d[Placebo]           -0.19    0.00 0.05     -0.28     -0.22     -0.19     -0.16     -0.09  1153
#> lp__             -37970.28    0.09 3.70 -37978.56 -37972.56 -37969.93 -37967.62 -37964.18  1766
#>                  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 Thu Feb 24 09:01:36 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.00 0.08     -0.50     -0.38     -0.33     -0.28     -0.18  1754
#> d[ARB]               -0.40    0.00 0.09     -0.59     -0.46     -0.40     -0.34     -0.22  2327
#> d[CCB]               -0.17    0.00 0.06     -0.29     -0.21     -0.17     -0.13     -0.04  2066
#> d[Diuretic]           0.07    0.00 0.09     -0.10      0.01      0.07      0.13      0.25  1894
#> d[Placebo]           -0.22    0.00 0.09     -0.39     -0.27     -0.21     -0.16     -0.05  1695
#> lp__             -37980.96    0.23 6.86 -37994.95 -37985.38 -37980.65 -37976.20 -37968.17   907
#> tau                   0.13    0.00 0.04      0.06      0.10      0.12      0.15      0.23  1000
#>                  Rhat
#> d[ACE Inhibitor]    1
#> d[ARB]              1
#> d[CCB]              1
#> d[Diuretic]         1
#> d[Placebo]          1
#> lp__                1
#> tau                 1
#> 
#> Samples were drawn using NUTS(diag_e) at Thu Feb 24 09:02:40 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.2 (on 48 data points)
#>                pD: 26.9
#>               DIC: 105.1
(dic_RE <- dic(db_fit_RE))
#> Residual deviance: 53.6 (on 48 data points)
#>                pD: 38
#>               DIC: 91.6

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.10 -0.06 -0.02  0.06     1538     2391    1
#> d[ACE Inhibitor] -0.36 0.05 -0.46 -0.40 -0.36 -0.32 -0.25     4445     3216    1
#> d[ARB]           -0.45 0.06 -0.58 -0.50 -0.45 -0.41 -0.33     3389     2994    1
#> d[CCB]           -0.25 0.05 -0.36 -0.29 -0.25 -0.22 -0.15     2988     2915    1
#> d[Placebo]       -0.25 0.06 -0.36 -0.29 -0.25 -0.21 -0.13     4066     3286    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     1929     2299    1
#> d[ACE Inhibitor] -0.40 0.09 -0.57 -0.46 -0.40 -0.35 -0.23     4074     3545    1
#> d[ARB]           -0.47 0.11 -0.70 -0.54 -0.47 -0.40 -0.26     2845     2400    1
#> d[CCB]           -0.24 0.09 -0.41 -0.29 -0.24 -0.18 -0.07     3710     2765    1
#> d[Placebo]       -0.29 0.09 -0.48 -0.35 -0.29 -0.23 -0.11     3224     2989    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.08  0.23     3805     3852    1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.01 0.02 0.03 0.06  0.17     3804     3853    1
#> pred[New 1: ARB]           0.04 0.04 0.00 0.02 0.03 0.05  0.16     3798     3853    1
#> pred[New 1: CCB]           0.05 0.05 0.01 0.02 0.03 0.06  0.19     3801     3852    1
#> pred[New 1: Diuretic]      0.06 0.06 0.01 0.02 0.04 0.08  0.24     3792     3776    1
#> pred[New 1: Placebo]       0.05 0.05 0.01 0.02 0.03 0.07  0.19     3800     3814    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.24     3706     3654    1
#> pred[New 1: ACE Inhibitor] 0.04 0.05 0.00 0.02 0.03 0.06  0.18     3811     3775    1
#> pred[New 1: ARB]           0.04 0.05 0.00 0.02 0.03 0.05  0.17     3797     3654    1
#> pred[New 1: CCB]           0.05 0.05 0.01 0.02 0.03 0.07  0.21     3752     3819    1
#> pred[New 1: Diuretic]      0.07 0.07 0.01 0.02 0.04 0.08  0.26     3812     3781    1
#> pred[New 1: Placebo]       0.05 0.05 0.01 0.02 0.03 0.06  0.20     3800     3800    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     5546     2966    1
#> pred[AASK: ACE Inhibitor] 0.12 0.01 0.10 0.12 0.12 0.13  0.15     4237     3290    1
#> pred[AASK: ARB]           0.12 0.01 0.09 0.11 0.12 0.13  0.15     4376     3195    1
#> pred[AASK: CCB]           0.14 0.01 0.12 0.13 0.14 0.15  0.18     5175     3244    1
#> pred[AASK: Diuretic]      0.18 0.02 0.14 0.17 0.18 0.19  0.22     4109     3286    1
#> pred[AASK: Placebo]       0.14 0.02 0.11 0.13 0.14 0.15  0.17     3933     3479    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     2677     2696    1
#> pred[ALLHAT: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03  0.04     4512     3023    1
#> pred[ALLHAT: ARB]           0.03 0.00 0.02 0.03 0.03 0.03  0.04     3806     2559    1
#> pred[ALLHAT: CCB]           0.04 0.00 0.03 0.03 0.04 0.04  0.05     3661     2641    1
#> pred[ALLHAT: Diuretic]      0.05 0.01 0.04 0.04 0.05 0.05  0.06     4590     3276    1
#> pred[ALLHAT: Placebo]       0.03 0.00 0.03 0.03 0.03 0.04  0.04     3985     3109    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     6753     3493    1
#> pred[ALPINE: ACE Inhibitor] 0.02 0.01 0.01 0.01 0.02 0.02  0.03     7169     3164    1
#> pred[ALPINE: ARB]           0.02 0.01 0.01 0.01 0.02 0.02  0.03     7692     3430    1
#> pred[ALPINE: CCB]           0.02 0.01 0.01 0.02 0.02 0.03  0.04     7629     3415    1
#> pred[ALPINE: Diuretic]      0.03 0.01 0.01 0.02 0.03 0.03  0.05     7574     3275    1
#> pred[ALPINE: Placebo]       0.02 0.01 0.01 0.02 0.02 0.03  0.04     6979     3264    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     2720     2523    1
#> pred[ANBP-2: ACE Inhibitor] 0.05 0.01 0.04 0.04 0.05 0.05  0.06     4477     2923    1
#> pred[ANBP-2: ARB]           0.05 0.01 0.03 0.04 0.05 0.05  0.06     3441     2776    1
#> pred[ANBP-2: CCB]           0.06 0.01 0.04 0.05 0.06 0.06  0.08     4177     2732    1
#> pred[ANBP-2: Diuretic]      0.07 0.01 0.06 0.07 0.07 0.08  0.09     4875     3071    1
#> pred[ANBP-2: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     4397     2813    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     5935     2912    1
#> pred[ASCOT: ACE Inhibitor] 0.08 0.01 0.07 0.08 0.08 0.09  0.10     2114     2646    1
#> pred[ASCOT: ARB]           0.08 0.01 0.06 0.07 0.08 0.08  0.09     2697     2768    1
#> pred[ASCOT: CCB]           0.10 0.01 0.08 0.09 0.09 0.10  0.11     2463     2578    1
#> pred[ASCOT: Diuretic]      0.12 0.01 0.10 0.11 0.12 0.13  0.14     2152     2296    1
#> pred[ASCOT: Placebo]       0.09 0.01 0.08 0.09 0.09 0.10  0.11     2073     2734    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     5978     2988    1
#> pred[CAPPP: ACE Inhibitor] 0.05 0.00 0.05 0.05 0.05 0.06  0.06     2146     2427    1
#> pred[CAPPP: ARB]           0.05 0.01 0.04 0.05 0.05 0.05  0.06     2698     2524    1
#> pred[CAPPP: CCB]           0.06 0.00 0.05 0.06 0.06 0.07  0.07     3199     2737    1
#> pred[CAPPP: Diuretic]      0.08 0.01 0.07 0.08 0.08 0.08  0.10     2656     2546    1
#> pred[CAPPP: Placebo]       0.06 0.01 0.05 0.06 0.06 0.06  0.07     2096     2430    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     2895     2689    1
#> pred[CHARM: ACE Inhibitor] 0.07 0.01 0.05 0.06 0.07 0.07  0.09     3963     2730    1
#> pred[CHARM: ARB]           0.06 0.01 0.05 0.06 0.06 0.07  0.08     4906     2675    1
#> pred[CHARM: CCB]           0.08 0.01 0.06 0.07 0.08 0.08  0.10     3588     3043    1
#> pred[CHARM: Diuretic]      0.10 0.01 0.07 0.09 0.10 0.11  0.13     3338     2436    1
#> pred[CHARM: Placebo]       0.07 0.01 0.06 0.07 0.07 0.08  0.09     4182     2499    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     2464     2692    1
#> pred[DREAM: ACE Inhibitor] 0.17 0.02 0.13 0.16 0.17 0.18  0.21     4081     2921    1
#> pred[DREAM: ARB]           0.16 0.02 0.12 0.15 0.16 0.17  0.21     3918     2872    1
#> pred[DREAM: CCB]           0.20 0.03 0.15 0.18 0.19 0.21  0.25     3306     2635    1
#> pred[DREAM: Diuretic]      0.24 0.03 0.19 0.22 0.24 0.26  0.31     3639     2878    1
#> pred[DREAM: Placebo]       0.19 0.02 0.15 0.17 0.19 0.20  0.23     4311     3245    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     3306     2558    1
#> pred[EWPH: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.04 0.05  0.06     5586     3060    1
#> pred[EWPH: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     4264     2740    1
#> pred[EWPH: CCB]           0.05 0.01 0.04 0.05 0.05 0.06  0.08     4634     3054    1
#> pred[EWPH: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.07  0.09     5009     2801    1
#> pred[EWPH: Placebo]       0.05 0.01 0.03 0.04 0.05 0.06  0.07     5237     3392    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     3065     2808    1
#> pred[FEVER: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03  0.04     4766     2864    1
#> pred[FEVER: ARB]           0.03 0.00 0.02 0.03 0.03 0.03  0.04     4377     2934    1
#> pred[FEVER: CCB]           0.04 0.00 0.03 0.03 0.03 0.04  0.05     4072     2662    1
#> pred[FEVER: Diuretic]      0.04 0.01 0.03 0.04 0.04 0.05  0.06     4503     2714    1
#> pred[FEVER: Placebo]       0.03 0.00 0.03 0.03 0.03 0.04  0.04     4965     2962    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     5323     3142    1
#> pred[HAPPHY: ACE Inhibitor] 0.02  0 0.01 0.02 0.02 0.02  0.02     3973     2821    1
#> pred[HAPPHY: ARB]           0.02  0 0.01 0.02 0.02 0.02  0.02     3969     2964    1
#> pred[HAPPHY: CCB]           0.02  0 0.02 0.02 0.02 0.02  0.03     4643     3141    1
#> pred[HAPPHY: Diuretic]      0.03  0 0.02 0.02 0.03 0.03  0.03     3330     2564    1
#> pred[HAPPHY: Placebo]       0.02  0 0.02 0.02 0.02 0.02  0.03     3752     3132    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     2739     2829    1
#> pred[HOPE: ACE Inhibitor] 0.04 0.01 0.03 0.04 0.04 0.05  0.06     4374     3142    1
#> pred[HOPE: ARB]           0.04 0.01 0.03 0.04 0.04 0.04  0.05     3739     2915    1
#> pred[HOPE: CCB]           0.05 0.01 0.04 0.04 0.05 0.05  0.07     3661     3011    1
#> pred[HOPE: Diuretic]      0.06 0.01 0.05 0.06 0.06 0.07  0.08     4236     3133    1
#> pred[HOPE: Placebo]       0.05 0.01 0.04 0.04 0.05 0.05  0.06     4868     3116    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     2936     2450    1
#> pred[INSIGHT: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05  0.06     4111     2755    1
#> pred[INSIGHT: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     3805     2569    1
#> pred[INSIGHT: CCB]           0.06 0.01 0.04 0.05 0.05 0.06  0.07     4139     2735    1
#> pred[INSIGHT: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.07  0.09     4434     2823    1
#> pred[INSIGHT: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     3736     2907    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     6939     3074    1
#> pred[INVEST: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06  0.07     2064     2287    1
#> pred[INVEST: ARB]           0.06 0.01 0.05 0.05 0.06 0.06  0.07     2583     2476    1
#> pred[INVEST: CCB]           0.07 0.00 0.06 0.07 0.07 0.07  0.08     2381     2715    1
#> pred[INVEST: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.09  0.11     2102     2423    1
#> pred[INVEST: Placebo]       0.07 0.01 0.06 0.06 0.07 0.07  0.08     1949     2244    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     6900     2461    1
#> pred[LIFE: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06  0.07     2487     2486    1
#> pred[LIFE: ARB]           0.06 0.01 0.05 0.05 0.06 0.06  0.07     2840     2477    1
#> pred[LIFE: CCB]           0.07 0.01 0.06 0.07 0.07 0.07  0.08     3480     2794    1
#> pred[LIFE: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.09  0.11     2540     2780    1
#> pred[LIFE: Placebo]       0.07 0.01 0.05 0.06 0.07 0.07  0.08     2344     2645    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     4761     3446    1
#> pred[MRC-E: ACE Inhibitor] 0.02  0 0.02 0.02 0.02 0.02  0.03     5734     3231    1
#> pred[MRC-E: ARB]           0.02  0 0.01 0.02 0.02 0.02  0.03     5443     3059    1
#> pred[MRC-E: CCB]           0.03  0 0.02 0.02 0.02 0.03  0.03     5575     3516    1
#> pred[MRC-E: Diuretic]      0.03  0 0.02 0.03 0.03 0.03  0.04     4679     3556    1
#> pred[MRC-E: Placebo]       0.02  0 0.02 0.02 0.02 0.03  0.03     5027     3382    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     6043     3426    1
#> pred[NORDIL: ACE Inhibitor] 0.04 0.00 0.03 0.03 0.04 0.04  0.04     2496     2784    1
#> pred[NORDIL: ARB]           0.03 0.00 0.03 0.03 0.03 0.04  0.04     3001     2817    1
#> pred[NORDIL: CCB]           0.04 0.00 0.04 0.04 0.04 0.04  0.05     3095     2883    1
#> pred[NORDIL: Diuretic]      0.05 0.01 0.04 0.05 0.05 0.06  0.06     2541     2918    1
#> pred[NORDIL: Placebo]       0.04 0.00 0.03 0.04 0.04 0.04  0.05     2284     2623    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     2516     2232    1
#> pred[PEACE: ACE Inhibitor] 0.10 0.01 0.08 0.09 0.10 0.11  0.13     4039     2685    1
#> pred[PEACE: ARB]           0.09 0.01 0.07 0.09 0.09 0.10  0.12     3813     2530    1
#> pred[PEACE: CCB]           0.12 0.02 0.09 0.11 0.12 0.13  0.15     3482     2585    1
#> pred[PEACE: Diuretic]      0.15 0.02 0.11 0.13 0.15 0.16  0.19     3327     2808    1
#> pred[PEACE: Placebo]       0.11 0.01 0.09 0.10 0.11 0.12  0.14     4851     2812    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     3209     2804    1
#> pred[SCOPE: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05  0.06     4589     3016    1
#> pred[SCOPE: ARB]           0.04 0.01 0.03 0.04 0.04 0.05  0.06     4873     2898    1
#> pred[SCOPE: CCB]           0.06 0.01 0.04 0.05 0.05 0.06  0.07     4042     2874    1
#> pred[SCOPE: Diuretic]      0.07 0.01 0.05 0.06 0.07 0.08  0.09     4106     2989    1
#> pred[SCOPE: Placebo]       0.05 0.01 0.04 0.05 0.05 0.06  0.07     5094     3301    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     2771     2575    1
#> pred[SHEP: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.07  0.08     4406     2748    1
#> pred[SHEP: ARB]           0.06 0.01 0.04 0.05 0.06 0.06  0.08     3727     2634    1
#> pred[SHEP: CCB]           0.07 0.01 0.05 0.07 0.07 0.08  0.10     3764     2748    1
#> pred[SHEP: Diuretic]      0.09 0.01 0.07 0.08 0.09 0.10  0.12     4299     2960    1
#> pred[SHEP: Placebo]       0.07 0.01 0.05 0.06 0.07 0.08  0.09     4533     3044    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     4255     2813    1
#> pred[STOP-2: ACE Inhibitor] 0.04 0.00 0.03 0.04 0.04 0.04  0.05     2933     2398    1
#> pred[STOP-2: ARB]           0.04 0.00 0.03 0.03 0.04 0.04  0.05     3183     2839    1
#> pred[STOP-2: CCB]           0.05 0.00 0.04 0.04 0.05 0.05  0.05     4075     2843    1
#> pred[STOP-2: Diuretic]      0.06 0.01 0.05 0.05 0.06 0.06  0.07     3413     3274    1
#> pred[STOP-2: Placebo]       0.04 0.00 0.03 0.04 0.04 0.05  0.05     2668     2759    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     3063     2312    1
#> pred[VALUE: ACE Inhibitor] 0.14 0.02 0.11 0.13 0.14 0.16  0.19     3813     2126    1
#> pred[VALUE: ARB]           0.14 0.02 0.11 0.13 0.14 0.14  0.17     4298     2402    1
#> pred[VALUE: CCB]           0.17 0.02 0.13 0.15 0.17 0.18  0.21     4328     2460    1
#> pred[VALUE: Diuretic]      0.21 0.03 0.16 0.19 0.21 0.22  0.27     3703     2242    1
#> pred[VALUE: Placebo]       0.16 0.02 0.12 0.15 0.16 0.17  0.21     3736     2536    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.42    5   5   5   5     6     2434       NA    1
#> rank[ACE Inhibitor] 1.84 0.54    1   2   2   2     3     3528     2788    1
#> rank[ARB]           1.27 0.52    1   1   1   1     3     3226     2638    1
#> rank[CCB]           3.70 0.52    3   3   4   4     4     3543     3127    1
#> rank[Diuretic]      5.80 0.41    5   6   6   6     6     2576       NA    1
#> rank[Placebo]       3.20 0.60    2   3   3   4     4     2890     2987    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.07      0.00      0.00       0.0
#> d[ARB]                0.76      0.21      0.02      0.00      0.00       0.0
#> d[CCB]                0.00      0.02      0.27      0.71      0.01       0.0
#> d[Diuretic]           0.00      0.00      0.00      0.00      0.19       0.8
#> d[Placebo]            0.01      0.08      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.01       0.8         1
#> d[ACE Inhibitor]      0.23      0.93      1.00      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.29      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.72      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://nicedsu.sites.sheffield.ac.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.↩︎