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  2854    1
#> d[ASPAC]          0.02    0.00 0.04     -0.06     -0.01      0.01      0.04      0.09  5508    1
#> d[PTCA]          -0.48    0.00 0.10     -0.67     -0.54     -0.48     -0.41     -0.28  4900    1
#> d[r-PA]          -0.12    0.00 0.06     -0.24     -0.16     -0.12     -0.08     -0.01  3832    1
#> d[SK + t-PA]     -0.05    0.00 0.05     -0.14     -0.08     -0.05     -0.02      0.04  5515    1
#> d[t-PA]           0.00    0.00 0.03     -0.06     -0.02      0.00      0.02      0.06  4812    1
#> d[TNK]           -0.17    0.00 0.08     -0.32     -0.22     -0.17     -0.12     -0.02  4019    1
#> d[UK]            -0.20    0.00 0.22     -0.64     -0.34     -0.20     -0.05      0.24  4953    1
#> lp__         -43042.64    0.14 5.48 -43054.41 -43046.12 -43042.33 -43038.87 -43032.84  1464    1
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Jun 23 15:45:27 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.4 (on 102 data points)
#>                pD: 58.3
#>               DIC: 163.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.26     -0.19     -0.16     -0.12     -0.06
#> d[ASPAC vs. SK]            0.01    0.00 0.04     -0.07     -0.02      0.01      0.03      0.08
#> d[PTCA vs. SK]            -0.67    0.00 0.19     -1.05     -0.79     -0.67     -0.54     -0.31
#> 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.14     -0.08     -0.04     -0.01      0.05
#> d[t-PA vs. SK]             0.00    0.00 0.03     -0.06     -0.03      0.00      0.02      0.05
#> d[UK vs. SK]              -0.38    0.01 0.52     -1.41     -0.74     -0.37     -0.03      0.66
#> d[ASPAC vs. Acc t-PA]      1.40    0.01 0.41      0.63      1.13      1.39      1.65      2.25
#> d[PTCA vs. Acc t-PA]      -0.22    0.00 0.12     -0.45     -0.30     -0.22     -0.14      0.02
#> d[r-PA vs. Acc t-PA]       0.02    0.00 0.06     -0.11     -0.02      0.02      0.06      0.15
#> d[TNK vs. Acc t-PA]        0.00    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.56     -0.10      0.14      0.37      0.83
#> d[t-PA vs. ASPAC]          0.29    0.01 0.35     -0.40      0.05      0.30      0.53      0.98
#> d[t-PA vs. PTCA]           0.54    0.01 0.42     -0.24      0.26      0.53      0.82      1.41
#> d[UK vs. t-PA]            -0.30    0.00 0.35     -0.98     -0.53     -0.29     -0.07      0.38
#> lp__                  -43039.83    0.14 5.80 -43051.91 -43043.40 -43039.54 -43035.74 -43029.42
#>                       n_eff Rhat
#> d[Acc t-PA vs. SK]     6551    1
#> d[ASPAC vs. SK]        4496    1
#> d[PTCA vs. SK]         4819    1
#> d[r-PA vs. SK]         5296    1
#> d[SK + t-PA vs. SK]    5643    1
#> d[t-PA vs. SK]         4541    1
#> d[UK vs. SK]           5805    1
#> d[ASPAC vs. Acc t-PA]  3407    1
#> d[PTCA vs. Acc t-PA]   4676    1
#> d[r-PA vs. Acc t-PA]   5362    1
#> d[TNK vs. Acc t-PA]    4900    1
#> d[UK vs. Acc t-PA]     4778    1
#> d[t-PA vs. ASPAC]      4430    1
#> d[t-PA vs. PTCA]       3684    1
#> d[UK vs. t-PA]         5052    1
#> lp__                   1633    1
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Jun 23 15:45:52 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.4 (on 102 data points)
#>                pD: 58.3
#>               DIC: 163.7
(dic_ume <- dic(thrombo_fit_ume))
#> Residual deviance: 99.8 (on 102 data points)
#>                pD: 66.1
#>               DIC: 165.9

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     2901     3133    1
#> d[ASPAC vs. SK]            0.02 0.04 -0.06 -0.01  0.01  0.04  0.09     5530     3200    1
#> d[PTCA vs. SK]            -0.48 0.10 -0.67 -0.54 -0.48 -0.41 -0.28     4817     3758    1
#> d[r-PA vs. SK]            -0.12 0.06 -0.24 -0.16 -0.12 -0.08 -0.01     4024     3719    1
#> d[SK + t-PA vs. SK]       -0.05 0.05 -0.14 -0.08 -0.05 -0.02  0.04     5652     3118    1
#> d[t-PA vs. SK]             0.00 0.03 -0.06 -0.02  0.00  0.02  0.06     4772     3321    1
#> d[TNK vs. SK]             -0.17 0.08 -0.32 -0.22 -0.17 -0.12 -0.02     4058     3565    1
#> d[UK vs. SK]              -0.20 0.22 -0.64 -0.34 -0.20 -0.05  0.24     4959     3483    1
#> d[ASPAC vs. Acc t-PA]      0.19 0.06  0.09  0.15  0.19  0.23  0.30     3839     3641    1
#> d[PTCA vs. Acc t-PA]      -0.30 0.10 -0.49 -0.37 -0.30 -0.23 -0.10     6144     3368    1
#> d[r-PA vs. Acc t-PA]       0.06 0.06 -0.05  0.02  0.06  0.09  0.16     6513     3400    1
#> d[SK + t-PA vs. Acc t-PA]  0.13 0.05  0.02  0.09  0.13  0.17  0.23     5426     3269    1
#> d[t-PA vs. Acc t-PA]       0.18 0.05  0.08  0.14  0.18  0.21  0.28     3375     3353    1
#> d[TNK vs. Acc t-PA]        0.01 0.06 -0.12 -0.04  0.01  0.05  0.13     5328     3695    1
#> d[UK vs. Acc t-PA]        -0.02 0.22 -0.47 -0.17 -0.02  0.13  0.42     5080     3620    1
#> d[PTCA vs. ASPAC]         -0.49 0.11 -0.70 -0.56 -0.49 -0.42 -0.28     4989     3549    1
#> d[r-PA vs. ASPAC]         -0.14 0.07 -0.27 -0.19 -0.14 -0.09  0.00     4731     3681    1
#> d[SK + t-PA vs. ASPAC]    -0.06 0.06 -0.18 -0.10 -0.06 -0.02  0.05     5803     3238    1
#> d[t-PA vs. ASPAC]         -0.01 0.04 -0.09 -0.04 -0.01  0.01  0.06     6881     3314    1
#> d[TNK vs. ASPAC]          -0.19 0.09 -0.35 -0.25 -0.18 -0.13 -0.02     4370     3458    1
#> d[UK vs. ASPAC]           -0.21 0.22 -0.66 -0.36 -0.22 -0.06  0.23     5018     3331    1
#> d[r-PA vs. PTCA]           0.35 0.11  0.13  0.28  0.35  0.43  0.57     6127     3369    1
#> d[SK + t-PA vs. PTCA]      0.43 0.11  0.22  0.35  0.43  0.50  0.64     6138     3735    1
#> d[t-PA vs. PTCA]           0.48 0.11  0.27  0.41  0.48  0.55  0.69     4833     3746    1
#> d[TNK vs. PTCA]            0.30 0.12  0.08  0.22  0.31  0.38  0.53     6686     2865    1
#> d[UK vs. PTCA]             0.28 0.24 -0.20  0.11  0.27  0.44  0.75     4923     3498    1
#> d[SK + t-PA vs. r-PA]      0.07 0.07 -0.06  0.03  0.07  0.12  0.21     6290     2917    1
#> d[t-PA vs. r-PA]           0.12 0.07  0.00  0.08  0.12  0.17  0.26     4463     3566    1
#> d[TNK vs. r-PA]           -0.05 0.08 -0.21 -0.10 -0.05  0.01  0.12     7307     3006    1
#> d[UK vs. r-PA]            -0.08 0.23 -0.53 -0.23 -0.08  0.08  0.38     5236     3540    1
#> d[t-PA vs. SK + t-PA]      0.05 0.05 -0.05  0.01  0.05  0.09  0.16     5526     3291    1
#> d[TNK vs. SK + t-PA]      -0.12 0.08 -0.28 -0.18 -0.12 -0.07  0.04     5909     3552    1
#> d[UK vs. SK + t-PA]       -0.15 0.23 -0.60 -0.30 -0.15  0.00  0.30     5087     3424    1
#> d[TNK vs. t-PA]           -0.17 0.08 -0.33 -0.23 -0.17 -0.12 -0.01     4328     3467    1
#> d[UK vs. t-PA]            -0.20 0.22 -0.65 -0.35 -0.20 -0.05  0.24     4983     3422    1
#> d[UK vs. TNK]             -0.03 0.23 -0.48 -0.18 -0.03  0.13  0.42     5174     3448    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.46 0.97    6   7   7   8     9     4106       NA    1
#> rank[Acc t-PA]  3.17 0.81    2   3   3   4     5     3554     3650    1
#> rank[ASPAC]     7.95 1.15    5   7   8   9     9     4779       NA    1
#> rank[PTCA]      1.13 0.34    1   1   1   1     2     3346     3336    1
#> rank[r-PA]      4.40 1.17    2   4   4   5     7     5420     3661    1
#> rank[SK + t-PA] 5.99 1.21    4   5   6   6     9     5524       NA    1
#> rank[t-PA]      7.48 1.08    5   7   8   8     9     4501       NA    1
#> rank[TNK]       3.49 1.24    2   3   3   4     6     4909     3177    1
#> rank[UK]        3.94 2.72    1   2   2   6     9     4754       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.38      0.31
#> d[Acc t-PA]       0.00      0.21      0.46      0.29      0.04      0.00      0.00      0.00
#> d[ASPAC]          0.00      0.00      0.00      0.00      0.03      0.10      0.19      0.25
#> d[PTCA]           0.88      0.12      0.00      0.00      0.00      0.00      0.00      0.00
#> d[r-PA]           0.00      0.06      0.14      0.31      0.38      0.08      0.01      0.01
#> d[SK + t-PA]      0.00      0.00      0.01      0.06      0.26      0.45      0.10      0.07
#> d[t-PA]           0.00      0.00      0.00      0.00      0.03      0.15      0.29      0.33
#> d[TNK]            0.00      0.23      0.31      0.26      0.14      0.03      0.01      0.00
#> d[UK]             0.12      0.38      0.07      0.08      0.10      0.05      0.02      0.02
#>              p_rank[9]
#> d[SK]             0.16
#> d[Acc t-PA]       0.00
#> d[ASPAC]          0.43
#> d[PTCA]           0.00
#> d[r-PA]           0.01
#> d[SK + t-PA]      0.05
#> d[t-PA]           0.19
#> d[TNK]            0.01
#> d[UK]             0.15
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.15      0.53      0.84
#> d[Acc t-PA]       0.00      0.21      0.67      0.96      1.00      1.00      1.00      1.00
#> d[ASPAC]          0.00      0.00      0.00      0.00      0.03      0.13      0.32      0.57
#> d[PTCA]           0.88      1.00      1.00      1.00      1.00      1.00      1.00      1.00
#> d[r-PA]           0.00      0.06      0.20      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.78      0.88      0.95
#> d[t-PA]           0.00      0.00      0.00      0.00      0.04      0.19      0.48      0.81
#> d[TNK]            0.00      0.24      0.55      0.81      0.95      0.98      0.99      0.99
#> d[UK]             0.12      0.50      0.57      0.65      0.75      0.80      0.82      0.85
#>              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://www.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.