Example: Thrombolytic treatments

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

This vignette describes the analysis of 50 trials of 8 thrombolytic drugs (streptokinase, SK; alteplase, t-PA; accelerated alteplase, Acc t-PA; streptokinase plus alteplase, SK+tPA; reteplase, r-PA; tenocteplase, TNK; urokinase, UK; anistreptilase, ASPAC) plus per-cutaneous transluminal coronary angioplasty (PTCA) (Boland et al. 2003; Lu and Ades 2006; Dias et al. 2011). The number of deaths in 30 or 35 days following acute myocardial infarction are recorded. The data are available in this package as thrombolytics:

head(thrombolytics)
#>   studyn trtn      trtc    r     n
#> 1      1    1        SK 1472 20251
#> 2      1    3  Acc t-PA  652 10396
#> 3      1    4 SK + t-PA  723 10374
#> 4      2    1        SK    9   130
#> 5      2    2      t-PA    6   123
#> 6      3    1        SK    5    63

Setting up the network

We begin by setting up the network. We have arm-level count data giving the number of deaths (r) out of the total (n) in each arm, so we use the function set_agd_arm(). By default, SK is set as the network reference treatment.

thrombo_net <- set_agd_arm(thrombolytics, 
                           study = studyn,
                           trt = trtc,
                           r = r, 
                           n = n)
thrombo_net
#> A network with 50 AgD studies (arm-based).
#> 
#> ------------------------------------------------------- AgD studies (arm-based) ---- 
#>  Study Treatments                  
#>  1     3: SK | Acc t-PA | SK + t-PA
#>  2     2: SK | t-PA                
#>  3     2: SK | t-PA                
#>  4     2: SK | t-PA                
#>  5     2: SK | t-PA                
#>  6     3: SK | t-PA | ASPAC        
#>  7     2: SK | t-PA                
#>  8     2: SK | t-PA                
#>  9     2: SK | t-PA                
#>  10    2: SK | SK + t-PA           
#>  ... plus 40 more studies
#> 
#>  Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 9
#> Total number of studies: 50
#> Reference treatment is: SK
#> Network is connected

Plot the network structure.

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

Fixed effects NMA

Following TSD 4 (Dias et al. 2011), we fit a fixed effects NMA 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. By default, this will use a Binomial likelihood and a logit link function, auto-detected from the data.

thrombo_fit <- nma(thrombo_net, 
                   trt_effects = "fixed",
                   prior_intercept = normal(scale = 100),
                   prior_trt = normal(scale = 100))
#> Note: Setting "SK" as the network reference treatment.

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

thrombo_fit
#> A fixed effects NMA with a binomial likelihood (logit link).
#> 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 Rhat
#> d[Acc t-PA]      -0.18    0.00 0.04     -0.26     -0.21     -0.18     -0.15     -0.09  2637    1
#> d[ASPAC]          0.02    0.00 0.04     -0.05     -0.01      0.02      0.04      0.09  4766    1
#> d[PTCA]          -0.47    0.00 0.10     -0.67     -0.54     -0.47     -0.40     -0.28  4097    1
#> d[r-PA]          -0.12    0.00 0.06     -0.24     -0.16     -0.13     -0.08     -0.01  3645    1
#> d[SK + t-PA]     -0.05    0.00 0.05     -0.14     -0.08     -0.05     -0.02      0.04  4996    1
#> d[t-PA]           0.00    0.00 0.03     -0.05     -0.02      0.00      0.02      0.06  5186    1
#> d[TNK]           -0.17    0.00 0.08     -0.32     -0.22     -0.17     -0.12     -0.02  3536    1
#> d[UK]            -0.20    0.00 0.22     -0.64     -0.35     -0.20     -0.06      0.22  4572    1
#> lp__         -43042.88    0.14 5.42 -43054.30 -43046.35 -43042.59 -43039.14 -43032.96  1564    1
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Dec 01 17:43:40 2020.
#> 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(thrombo_fit, pars = c("d", "mu"))

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

plot_prior_posterior(thrombo_fit, prior = "trt")

Model fit can be checked using the dic() function

(dic_consistency <- dic(thrombo_fit))
#> Residual deviance: 105.9 (on 102 data points)
#>                pD: 58.7
#>               DIC: 164.7

and the residual deviance contributions examined with the corresponding plot() method.

plot(dic_consistency)

There are a number of points which are not very well fit by the model, having posterior mean residual deviance contributions greater than 1.

Checking for inconsistency

We fit an unrelated mean effects (UME) model (Dias et al. 2011) to assess the consistency assumption. Again, we use the function nma(), but now with the argument consistency = "ume".

thrombo_fit_ume <- nma(thrombo_net, 
                  consistency = "ume",
                  trt_effects = "fixed",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100))
#> Note: Setting "SK" as the network reference treatment.
thrombo_fit_ume
#> A fixed effects NMA with a binomial likelihood (logit link).
#> An inconsistency model ('ume') was fitted.
#> 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%
#> d[Acc t-PA vs. SK]        -0.16    0.00 0.05     -0.25     -0.19     -0.16     -0.12     -0.06
#> d[ASPAC vs. SK]            0.01    0.00 0.04     -0.06     -0.02      0.01      0.03      0.07
#> d[PTCA vs. SK]            -0.67    0.00 0.19     -1.05     -0.79     -0.67     -0.54     -0.32
#> d[r-PA vs. SK]            -0.06    0.00 0.09     -0.24     -0.12     -0.06      0.00      0.12
#> d[SK + t-PA vs. SK]       -0.04    0.00 0.05     -0.13     -0.07     -0.04     -0.01      0.05
#> d[t-PA vs. SK]             0.00    0.00 0.03     -0.06     -0.02      0.00      0.02      0.06
#> d[UK vs. SK]              -0.37    0.01 0.52     -1.43     -0.73     -0.36     -0.02      0.63
#> d[ASPAC vs. Acc t-PA]      1.40    0.01 0.41      0.65      1.11      1.38      1.67      2.24
#> d[PTCA vs. Acc t-PA]      -0.21    0.00 0.12     -0.44     -0.29     -0.21     -0.13      0.02
#> d[r-PA vs. Acc t-PA]       0.02    0.00 0.07     -0.12     -0.03      0.02      0.06      0.15
#> d[TNK vs. Acc t-PA]        0.01    0.00 0.06     -0.12     -0.04      0.01      0.05      0.13
#> d[UK vs. Acc t-PA]         0.14    0.01 0.36     -0.57     -0.11      0.13      0.37      0.88
#> d[t-PA vs. ASPAC]          0.29    0.01 0.35     -0.38      0.05      0.28      0.52      0.97
#> d[t-PA vs. PTCA]           0.55    0.01 0.41     -0.26      0.27      0.54      0.81      1.40
#> d[UK vs. t-PA]            -0.30    0.00 0.35     -1.01     -0.53     -0.29     -0.06      0.38
#> lp__                  -43039.70    0.15 5.73 -43051.69 -43043.48 -43039.42 -43035.65 -43029.53
#>                       n_eff Rhat
#> d[Acc t-PA vs. SK]     6165    1
#> d[ASPAC vs. SK]        5091    1
#> d[PTCA vs. SK]         5292    1
#> d[r-PA vs. SK]         6223    1
#> d[SK + t-PA vs. SK]    5920    1
#> d[t-PA vs. SK]         4464    1
#> d[UK vs. SK]           5527    1
#> d[ASPAC vs. Acc t-PA]  3611    1
#> d[PTCA vs. Acc t-PA]   5140    1
#> d[r-PA vs. Acc t-PA]   5807    1
#> d[TNK vs. Acc t-PA]    5304    1
#> d[UK vs. Acc t-PA]     4518    1
#> d[t-PA vs. ASPAC]      4716    1
#> d[t-PA vs. PTCA]       3786    1
#> d[UK vs. t-PA]         5364    1
#> lp__                   1515    1
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Dec 01 17:44:08 2020.
#> 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).

Comparing the model fit statistics

dic_consistency
#> Residual deviance: 105.9 (on 102 data points)
#>                pD: 58.7
#>               DIC: 164.7
(dic_ume <- dic(thrombo_fit_ume))
#> Residual deviance: 99.6 (on 102 data points)
#>                pD: 65.8
#>               DIC: 165.4

Whilst the UME model fits the data better, having a lower residual deviance, the additional parameters in the UME model mean that the DIC is very similar between both models. However, it is also important to examine the individual contributions to model fit of each data point under the two models (a so-called “dev-dev” plot). Passing two nma_dic objects produced by the dic() function to the plot() method produces this dev-dev plot:

plot(dic_consistency, dic_ume, show_uncertainty = FALSE)

The four points lying in the lower right corner of the plot have much lower posterior mean residual deviance under the UME model, indicating that these data are potentially inconsistent. These points correspond to trials 44 and 45, the only two trials comparing Acc t-PA to ASPAC. The ASPAC vs. Acc t-PA estimates are very different under the consistency model and inconsistency (UME) model, suggesting that these two trials may be systematically different from the others in the network.

Further results

Relative effects for all pairwise contrasts between treatments can be produced using the relative_effects() function, with all_contrasts = TRUE.

(thrombo_releff <- relative_effects(thrombo_fit, all_contrasts = TRUE))
#>                            mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Acc t-PA vs. SK]        -0.18 0.04 -0.26 -0.21 -0.18 -0.15 -0.09     2662     2924    1
#> d[ASPAC vs. SK]            0.02 0.04 -0.05 -0.01  0.02  0.04  0.09     4827     3613    1
#> d[PTCA vs. SK]            -0.47 0.10 -0.67 -0.54 -0.47 -0.40 -0.28     4006     3528    1
#> d[r-PA vs. SK]            -0.12 0.06 -0.24 -0.16 -0.13 -0.08 -0.01     3679     3135    1
#> d[SK + t-PA vs. SK]       -0.05 0.05 -0.14 -0.08 -0.05 -0.02  0.04     5143     3379    1
#> d[t-PA vs. SK]             0.00 0.03 -0.05 -0.02  0.00  0.02  0.06     5023     3454    1
#> d[TNK vs. SK]             -0.17 0.08 -0.32 -0.22 -0.17 -0.12 -0.02     3555     3171    1
#> d[UK vs. SK]              -0.20 0.22 -0.64 -0.35 -0.20 -0.06  0.22     4510     3174    1
#> d[ASPAC vs. Acc t-PA]      0.20 0.05  0.09  0.16  0.20  0.23  0.30     3442     3550    1
#> d[PTCA vs. Acc t-PA]      -0.30 0.10 -0.49 -0.36 -0.30 -0.23 -0.10     5330     2920    1
#> d[r-PA vs. Acc t-PA]       0.05 0.06 -0.05  0.02  0.05  0.09  0.16     5481     3574    1
#> d[SK + t-PA vs. Acc t-PA]  0.13 0.05  0.02  0.09  0.13  0.17  0.24     4619     3591    1
#> d[t-PA vs. Acc t-PA]       0.18 0.05  0.08  0.15  0.18  0.21  0.28     3184     3534    1
#> d[TNK vs. Acc t-PA]        0.01 0.06 -0.12 -0.04  0.01  0.05  0.14     5634     2994    1
#> d[UK vs. Acc t-PA]        -0.02 0.22 -0.46 -0.17 -0.02  0.13  0.40     4611     3140    1
#> d[PTCA vs. ASPAC]         -0.49 0.11 -0.70 -0.57 -0.49 -0.42 -0.29     4214     3607    1
#> d[r-PA vs. ASPAC]         -0.14 0.07 -0.27 -0.19 -0.14 -0.09 -0.01     4018     3505    1
#> d[SK + t-PA vs. ASPAC]    -0.07 0.06 -0.18 -0.11 -0.07 -0.03  0.05     5336     3548    1
#> d[t-PA vs. ASPAC]         -0.01 0.04 -0.08 -0.04 -0.01  0.01  0.06     7427     3154    1
#> d[TNK vs. ASPAC]          -0.19 0.08 -0.35 -0.25 -0.19 -0.13 -0.02     3955     3561    1
#> d[UK vs. ASPAC]           -0.22 0.22 -0.66 -0.37 -0.22 -0.07  0.21     4591     3400    1
#> d[r-PA vs. PTCA]           0.35 0.11  0.13  0.27  0.35  0.42  0.57     4933     3317    1
#> d[SK + t-PA vs. PTCA]      0.42 0.11  0.22  0.35  0.42  0.50  0.63     4904     3477    1
#> d[t-PA vs. PTCA]           0.48 0.10  0.28  0.40  0.48  0.55  0.68     4151     3202    1
#> d[TNK vs. PTCA]            0.30 0.12  0.07  0.22  0.30  0.39  0.54     6017     3118    1
#> d[UK vs. PTCA]             0.27 0.24 -0.21  0.11  0.27  0.44  0.74     4769     3265    1
#> d[SK + t-PA vs. r-PA]      0.07 0.07 -0.06  0.03  0.08  0.12  0.21     5322     3106    1
#> d[t-PA vs. r-PA]           0.13 0.07  0.00  0.08  0.13  0.17  0.26     3999     3326    1
#> d[TNK vs. r-PA]           -0.05 0.08 -0.21 -0.10 -0.05  0.01  0.12     6474     3329    1
#> d[UK vs. r-PA]            -0.08 0.23 -0.52 -0.23 -0.08  0.08  0.36     4894     3441    1
#> d[t-PA vs. SK + t-PA]      0.05 0.05 -0.05  0.02  0.05  0.09  0.16     5094     3741    1
#> d[TNK vs. SK + t-PA]      -0.12 0.08 -0.29 -0.18 -0.12 -0.06  0.04     5838     2754    1
#> d[UK vs. SK + t-PA]       -0.15 0.22 -0.59 -0.30 -0.15  0.00  0.28     4749     3195    1
#> d[TNK vs. t-PA]           -0.17 0.08 -0.34 -0.23 -0.17 -0.12 -0.02     3679     3106    1
#> d[UK vs. t-PA]            -0.21 0.22 -0.65 -0.35 -0.21 -0.06  0.23     4583     3285    1
#> d[UK vs. TNK]             -0.03 0.23 -0.48 -0.19 -0.02  0.12  0.41     4847     3522    1
plot(thrombo_releff, ref_line = 0)

Treatment rankings, rank probabilities, and cumulative rank probabilities.

(thrombo_ranks <- posterior_ranks(thrombo_fit))
#>                 mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[SK]        7.43 0.96    6   7   7   8     9     3734       NA    1
#> rank[Acc t-PA]  3.17 0.81    2   3   3   4     5     4253     3236    1
#> rank[ASPAC]     8.02 1.11    5   7   8   9     9     4683       NA    1
#> rank[PTCA]      1.13 0.35    1   1   1   1     2     3560     2945    1
#> rank[r-PA]      4.38 1.16    2   4   4   5     7     5039     3589    1
#> rank[SK + t-PA] 5.95 1.20    4   5   6   6     9     5133       NA    1
#> rank[t-PA]      7.51 1.07    5   7   8   8     9     5013       NA    1
#> rank[TNK]       3.49 1.25    2   3   3   4     6     4786     3630    1
#> rank[UK]        3.91 2.67    1   2   3   5     9     4442       NA    1
plot(thrombo_ranks)

(thrombo_rankprobs <- posterior_rank_probs(thrombo_fit))
#>              p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6] p_rank[7] p_rank[8]
#> d[SK]             0.00      0.00      0.00      0.00      0.02      0.14      0.39      0.31
#> d[Acc t-PA]       0.00      0.21      0.46      0.29      0.05      0.00      0.00      0.00
#> d[ASPAC]          0.00      0.00      0.00      0.00      0.02      0.08      0.18      0.26
#> d[PTCA]           0.87      0.13      0.00      0.00      0.00      0.00      0.00      0.00
#> d[r-PA]           0.00      0.06      0.15      0.30      0.38      0.08      0.01      0.01
#> d[SK + t-PA]      0.00      0.00      0.01      0.06      0.26      0.46      0.09      0.06
#> d[t-PA]           0.00      0.00      0.00      0.00      0.03      0.14      0.30      0.33
#> d[TNK]            0.00      0.24      0.30      0.26      0.15      0.03      0.01      0.01
#> d[UK]             0.13      0.36      0.07      0.09      0.10      0.06      0.02      0.02
#>              p_rank[9]
#> d[SK]             0.15
#> d[Acc t-PA]       0.00
#> d[ASPAC]          0.45
#> d[PTCA]           0.00
#> d[r-PA]           0.01
#> d[SK + t-PA]      0.05
#> d[t-PA]           0.20
#> d[TNK]            0.00
#> d[UK]             0.14
plot(thrombo_rankprobs)

(thrombo_cumrankprobs <- posterior_rank_probs(thrombo_fit, cumulative = TRUE))
#>              p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6] p_rank[7] p_rank[8]
#> d[SK]             0.00      0.00      0.00      0.00      0.02      0.16      0.54      0.85
#> d[Acc t-PA]       0.00      0.21      0.67      0.95      1.00      1.00      1.00      1.00
#> d[ASPAC]          0.00      0.00      0.00      0.00      0.03      0.11      0.29      0.55
#> d[PTCA]           0.87      1.00      1.00      1.00      1.00      1.00      1.00      1.00
#> d[r-PA]           0.00      0.06      0.21      0.51      0.89      0.97      0.98      0.99
#> d[SK + t-PA]      0.00      0.00      0.01      0.07      0.33      0.79      0.88      0.95
#> d[t-PA]           0.00      0.00      0.00      0.00      0.03      0.17      0.47      0.80
#> d[TNK]            0.00      0.24      0.55      0.80      0.95      0.98      0.99      1.00
#> d[UK]             0.13      0.49      0.57      0.65      0.75      0.81      0.83      0.86
#>              p_rank[9]
#> d[SK]                1
#> d[Acc t-PA]          1
#> d[ASPAC]             1
#> d[PTCA]              1
#> d[r-PA]              1
#> d[SK + t-PA]         1
#> d[t-PA]              1
#> d[TNK]               1
#> d[UK]                1
plot(thrombo_cumrankprobs)

References

Boland, A., Y. Dundar, A. Bagust, A. Haycox, R. Hill, R. Mujica Mota, T. Walley, and R. Dickson. 2003. “Early Thrombolysis for the Treatment of Acute Myocardial Infarction: A Systematic Review and Economic Evaluation.” Health Technology Assessment 7 (15). https://doi.org/10.3310/hta7150.

Dias, S., N. J. Welton, A. J. Sutton, D. M. Caldwell, G. Lu, and A. E. Ades. 2011. “NICE DSU Technical Support Document 4: Inconsistency in Networks of Evidence Based on Randomised Controlled Trials.” National Institute for Health and Care Excellence. http://nicedsu.org.uk/.

Lu, G. B., and A. E. Ades. 2006. “Assessing Evidence Inconsistency in Mixed Treatment Comparisons.” Journal of the American Statistical Association 101 (474): 447–59. https://doi.org/10.1198/016214505000001302.