Example: Smoking cessation

This vignette describes the analysis of smoking cessation data (Hasselblad 1998), replicating the analysis in NICE Technical Support Document 4 (Dias et al. 2011). The data are available in this package as smoking:

Setting up the network

We begin by setting up the network. We have arm-level count data giving the number quitting smoking (r) out of the total (n) in each arm, so we use the function set_agd_arm(). Treatment “No intervention” is set as the network reference treatment.

smknet <- set_agd_arm(smoking, 
                      study = studyn,
                      trt = trtc,
                      r = r, 
                      n = n,
                      trt_ref = "No intervention")
smknet
#> A network with 24 AgD studies (arm-based).
#> 
#> ------------------------------------------------------- AgD studies (arm-based) ---- 
#>  Study Treatment arms                                                 
#>  1     3: No intervention | Group counselling | Individual counselling
#>  2     3: Group counselling | Individual counselling | Self-help      
#>  3     2: No intervention | Individual counselling                    
#>  4     2: No intervention | Individual counselling                    
#>  5     2: No intervention | Individual counselling                    
#>  6     2: No intervention | Individual counselling                    
#>  7     2: No intervention | Individual counselling                    
#>  8     2: No intervention | Individual counselling                    
#>  9     2: No intervention | Individual counselling                    
#>  10    2: No intervention | Self-help                                 
#>  ... plus 14 more studies
#> 
#>  Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 4
#> Total number of studies: 24
#> Reference treatment is: No intervention
#> Network is connected

Plot the network structure.

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

Random effects NMA

Following TSD 4, we fit a random effects NMA model, using the nma() function with trt_effects = "random". We use \(\mathrm{N}(0, 100^2)\) prior distributions for the treatment effects \(d_k\) and study-specific intercepts \(\mu_j\), and a \(\textrm{half-N}(5^2)\) prior distribution for the between-study 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.

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.

smkfit <- nma(smknet, 
              trt_effects = "random",
              prior_intercept = normal(scale = 100),
              prior_trt = normal(scale = 100),
              prior_het = normal(scale = 5))

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

smkfit
#> A random 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
#> d[Group counselling]          1.09    0.01 0.44     0.25     0.80     1.08     1.36     2.01  2016
#> d[Individual counselling]     0.84    0.01 0.24     0.39     0.69     0.83     0.99     1.33  1399
#> d[Self-help]                  0.50    0.01 0.40    -0.25     0.24     0.49     0.75     1.33  1867
#> lp__                      -5768.30    0.19 6.28 -5781.39 -5772.40 -5767.95 -5763.84 -5757.07  1133
#> tau                           0.83    0.01 0.18     0.55     0.71     0.81     0.94     1.25  1150
#>                           Rhat
#> d[Group counselling]         1
#> d[Individual counselling]    1
#> d[Self-help]                 1
#> lp__                         1
#> tau                          1
#> 
#> Samples were drawn using NUTS(diag_e) at Thu Feb 24 08:50:02 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(smkfit, pars = c("d", "tau", "mu", "delta"))

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

plot_prior_posterior(smkfit)

By default, this displays all model parameters given prior distributions (in this case \(d_k\), \(\mu_j\), and \(\tau\)), but this may be changed using the prior argument:

plot_prior_posterior(smkfit, prior = "het")

Model fit can be checked using the dic() function

(dic_consistency <- dic(smkfit))
#> Residual deviance: 54.4 (on 50 data points)
#>                pD: 44.1
#>               DIC: 98.6

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

plot(dic_consistency)

Overall model fit seems to be adequate, with almost all points showing good fit (mean residual deviance contribution of 1). The only two points with higher residual deviance (i.e. worse fit) correspond to the two zero counts in the data:

smoking[smoking$r == 0, ]
#>    studyn trtn            trtc r  n
#> 13      6    1 No intervention 0 33
#> 31     15    1 No intervention 0 20

Checking for inconsistency

Note: The results of the inconsistency models here are slightly different to those of Dias et al. (2010, 2011), although the overall conclusions are the same. This is due to the presence of multi-arm trials and a different ordering of treatments, meaning that inconsistency is parameterised differently within the multi-arm trials. The same results as Dias et al. are obtained if the network is instead set up with trtn as the treatment variable.

Unrelated mean effects

We first 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".

smkfit_ume <- nma(smknet, 
                  consistency = "ume",
                  trt_effects = "random",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100),
                  prior_het = normal(scale = 5))
smkfit_ume
#> A random 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%
#> d[Group counselling vs. No intervention]            1.15    0.02 0.80    -0.38     0.62     1.12
#> d[Individual counselling vs. No intervention]       0.90    0.01 0.27     0.38     0.72     0.90
#> d[Self-help vs. No intervention]                    0.32    0.01 0.59    -0.88    -0.05     0.33
#> d[Individual counselling vs. Group counselling]    -0.29    0.01 0.63    -1.51    -0.70    -0.30
#> d[Self-help vs. Group counselling]                 -0.64    0.01 0.72    -2.10    -1.11    -0.63
#> d[Self-help vs. Individual counselling]             0.11    0.02 1.07    -2.09    -0.55     0.13
#> lp__                                            -5765.19    0.23 6.32 -5778.33 -5769.26 -5765.01
#> tau                                                 0.94    0.01 0.23     0.59     0.78     0.91
#>                                                      75%    97.5% n_eff Rhat
#> d[Group counselling vs. No intervention]            1.64     2.84  2685 1.00
#> d[Individual counselling vs. No intervention]       1.07     1.45  1153 1.01
#> d[Self-help vs. No intervention]                    0.71     1.47  2131 1.00
#> d[Individual counselling vs. Group counselling]     0.13     0.99  2420 1.00
#> d[Self-help vs. Group counselling]                 -0.17     0.80  2537 1.00
#> d[Self-help vs. Individual counselling]             0.81     2.21  4019 1.00
#> lp__                                            -5760.75 -5753.64   728 1.00
#> tau                                                 1.07     1.48  1036 1.00
#> 
#> Samples were drawn using NUTS(diag_e) at Thu Feb 24 08:50:25 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).

Comparing the model fit statistics

dic_consistency
#> Residual deviance: 54.4 (on 50 data points)
#>                pD: 44.1
#>               DIC: 98.6
(dic_ume <- dic(smkfit_ume))
#> Residual deviance: 53.7 (on 50 data points)
#>                pD: 45.1
#>               DIC: 98.8

We see that there is little to choose between the two 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, point_alpha = 0.5, interval_alpha = 0.2)

All points lie roughly on the line of equality, so there is no evidence for inconsistency here.

Node-splitting

Another method for assessing inconsistency is node-splitting (Dias et al. 2011, 2010). Whereas the UME model assesses inconsistency globally, node-splitting assesses inconsistency locally for each potentially inconsistent comparison (those with both direct and indirect evidence) in turn.

Node-splitting can be performed using the nma() function with the argument consistency = "nodesplit". By default, all possible comparisons will be split (as determined by the get_nodesplits() function). Alternatively, a specific comparison or comparisons to split can be provided to the nodesplit argument.

smk_nodesplit <- nma(smknet, 
                     consistency = "nodesplit",
                     trt_effects = "random",
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_het = normal(scale = 5))
#> Fitting model 1 of 7, node-split: Group counselling vs. No intervention
#> Fitting model 2 of 7, node-split: Individual counselling vs. No intervention
#> Fitting model 3 of 7, node-split: Self-help vs. No intervention
#> Fitting model 4 of 7, node-split: Individual counselling vs. Group counselling
#> Fitting model 5 of 7, node-split: Self-help vs. Group counselling
#> Fitting model 6 of 7, node-split: Self-help vs. Individual counselling
#> Fitting model 7 of 7, consistency model

The summary() method summarises the node-splitting results, displaying the direct and indirect estimates \(d_\mathrm{dir}\) and \(d_\mathrm{ind}\) from each node-split model, the network estimate \(d_\mathrm{net}\) from the consistency model, the inconsistency factor \(\omega = d_\mathrm{dir} - d_\mathrm{ind}\), and a Bayesian \(p\)-value for inconsistency on each comparison. Since random effects models are fitted, the heterogeneity standard deviation \(\tau\) under each node-split model and under the consistency model is also displayed. The DIC model fit statistics are also provided.

summary(smk_nodesplit)
#> Node-splitting models fitted for 6 comparisons.
#> 
#> ------------------------------ Node-split Group counselling vs. No intervention ---- 
#> 
#>                  mean   sd  2.5%   25%   50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d_net            1.12 0.44  0.25  0.82  1.11 1.40  2.00     1757     2243    1
#> d_dir            1.07 0.76 -0.36  0.57  1.06 1.55  2.67     3314     2388    1
#> d_ind            1.15 0.54  0.09  0.80  1.15 1.50  2.21     1983     2400    1
#> omega           -0.08 0.91 -1.86 -0.68 -0.10 0.51  1.75     2530     2017    1
#> tau              0.87 0.21  0.55  0.73  0.85 0.99  1.37     1075     1539    1
#> tau_consistency  0.84 0.18  0.55  0.71  0.81 0.94  1.25     1378     2149    1
#> 
#> Residual deviance: 54 (on 50 data points)
#>                pD: 44.2
#>               DIC: 98.1
#> 
#> Bayesian p-value: 0.92
#> 
#> ------------------------- Node-split Individual counselling vs. No intervention ---- 
#> 
#>                 mean   sd  2.5%   25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d_net           0.84 0.24  0.38  0.69 0.84 1.00  1.32     1319     1943    1
#> d_dir           0.88 0.26  0.40  0.72 0.88 1.04  1.41     2707     3026    1
#> d_ind           0.58 0.67 -0.73  0.14 0.56 1.02  1.94     1886     2407    1
#> omega           0.30 0.70 -1.10 -0.16 0.31 0.76  1.65     1953     2485    1
#> tau             0.86 0.20  0.55  0.72 0.83 0.97  1.33     1374     2128    1
#> tau_consistency 0.84 0.18  0.55  0.71 0.81 0.94  1.25     1378     2149    1
#> 
#> Residual deviance: 54.1 (on 50 data points)
#>                pD: 44
#>               DIC: 98.1
#> 
#> Bayesian p-value: 0.66
#> 
#> -------------------------------------- Node-split Self-help vs. No intervention ---- 
#> 
#>                  mean   sd  2.5%   25%   50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d_net            0.50 0.40 -0.27  0.23  0.49 0.76  1.29     1736     2483    1
#> d_dir            0.34 0.55 -0.73 -0.01  0.33 0.69  1.45     2810     2556    1
#> d_ind            0.69 0.63 -0.60  0.29  0.68 1.10  1.94     1840     2108    1
#> omega           -0.35 0.84 -2.00 -0.90 -0.36 0.19  1.29     1822     1675    1
#> tau              0.87 0.20  0.57  0.73  0.85 0.97  1.34     1282     1926    1
#> tau_consistency  0.84 0.18  0.55  0.71  0.81 0.94  1.25     1378     2149    1
#> 
#> Residual deviance: 53.7 (on 50 data points)
#>                pD: 44.2
#>               DIC: 97.9
#> 
#> Bayesian p-value: 0.64
#> 
#> ----------------------- Node-split Individual counselling vs. Group counselling ---- 
#> 
#>                  mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d_net           -0.27 0.42 -1.11 -0.55 -0.27  0.00  0.55     2342     2520    1
#> d_dir           -0.12 0.50 -1.10 -0.45 -0.11  0.20  0.86     3754     2641    1
#> d_ind           -0.56 0.63 -1.85 -0.97 -0.56 -0.17  0.69     1771     1977    1
#> omega            0.44 0.68 -0.93  0.00  0.42  0.88  1.83     1885     2275    1
#> tau              0.86 0.20  0.56  0.72  0.83  0.97  1.34     1292     1572    1
#> tau_consistency  0.84 0.18  0.55  0.71  0.81  0.94  1.25     1378     2149    1
#> 
#> Residual deviance: 54.1 (on 50 data points)
#>                pD: 44.5
#>               DIC: 98.6
#> 
#> Bayesian p-value: 0.5
#> 
#> ------------------------------------ Node-split Self-help vs. Group counselling ---- 
#> 
#>                  mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d_net           -0.62 0.48 -1.57 -0.93 -0.62 -0.30  0.36     2688     2823    1
#> d_dir           -0.61 0.68 -1.97 -1.06 -0.61 -0.15  0.75     3445     2715    1
#> d_ind           -0.62 0.66 -1.97 -1.04 -0.60 -0.20  0.68     1556     2287    1
#> omega            0.01 0.89 -1.72 -0.59 -0.01  0.60  1.74     1831     2412    1
#> tau              0.87 0.20  0.56  0.73  0.84  0.98  1.34     1277     1936    1
#> tau_consistency  0.84 0.18  0.55  0.71  0.81  0.94  1.25     1378     2149    1
#> 
#> Residual deviance: 54.3 (on 50 data points)
#>                pD: 44.4
#>               DIC: 98.7
#> 
#> Bayesian p-value: 0.99
#> 
#> ------------------------------- Node-split Self-help vs. Individual counselling ---- 
#> 
#>                  mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d_net           -0.35 0.41 -1.18 -0.61 -0.34 -0.09  0.44     2313     2380    1
#> d_dir            0.09 0.64 -1.17 -0.34  0.08  0.50  1.32     3222     3077    1
#> d_ind           -0.62 0.52 -1.65 -0.96 -0.62 -0.30  0.39     1345     2257    1
#> omega            0.71 0.81 -0.90  0.19  0.72  1.23  2.30     1923     2532    1
#> tau              0.85 0.19  0.54  0.71  0.82  0.95  1.28     1206     2340    1
#> tau_consistency  0.84 0.18  0.55  0.71  0.81  0.94  1.25     1378     2149    1
#> 
#> Residual deviance: 54.2 (on 50 data points)
#>                pD: 44.4
#>               DIC: 98.5
#> 
#> Bayesian p-value: 0.37

The DIC of each inconsistency model is unchanged from the consistency model, no node-splits result in reduced heterogeneity standard deviation \(\tau\) compared to the consistency model, and the Bayesian \(p\)-values are all large. There is no evidence of inconsistency.

We can visually compare the posterior distributions of the direct, indirect, and network estimates using the plot() method. These are all in agreement; the posterior densities of the direct and indirect estimates overlap. Notice that there is not much indirect information for the Individual counselling vs. No intervention comparison, so the network (consistency) estimate is very similar to the direct estimate for this comparison.

plot(smk_nodesplit) +
  ggplot2::theme(legend.position = "bottom", legend.direct = "horizontal")

Further results

Pairwise relative effects, for all pairwise contrasts with all_contrasts = TRUE.

(smk_releff <- relative_effects(smkfit, all_contrasts = TRUE))
#>                                                  mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS
#> d[Group counselling vs. No intervention]         1.09 0.44  0.25  0.80  1.08  1.36  2.01     2052
#> d[Individual counselling vs. No intervention]    0.84 0.24  0.39  0.69  0.83  0.99  1.33     1406
#> d[Self-help vs. No intervention]                 0.50 0.40 -0.25  0.24  0.49  0.75  1.33     1896
#> d[Individual counselling vs. Group counselling] -0.25 0.42 -1.09 -0.52 -0.24  0.02  0.55     2875
#> d[Self-help vs. Group counselling]              -0.59 0.49 -1.59 -0.90 -0.59 -0.27  0.40     2594
#> d[Self-help vs. Individual counselling]         -0.33 0.41 -1.12 -0.61 -0.34 -0.07  0.49     2061
#>                                                 Tail_ESS Rhat
#> d[Group counselling vs. No intervention]            2414    1
#> d[Individual counselling vs. No intervention]       2030    1
#> d[Self-help vs. No intervention]                    2444    1
#> d[Individual counselling vs. Group counselling]     2799    1
#> d[Self-help vs. Group counselling]                  2529    1
#> d[Self-help vs. Individual counselling]             2234    1
plot(smk_releff, ref_line = 0)

Treatment rankings, rank probabilities, and cumulative rank probabilities. We set lower_better = FALSE since a higher log odds of cessation is better (the outcome is positive).

(smk_ranks <- posterior_ranks(smkfit, lower_better = FALSE))
#>                              mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[No intervention]        3.90 0.31    3   4   4   4     4     2139       NA    1
#> rank[Group counselling]      1.39 0.64    1   1   1   2     3     2780     2922    1
#> rank[Individual counselling] 1.93 0.63    1   2   2   2     3     2508     2740    1
#> rank[Self-help]              2.78 0.70    1   3   3   3     4     2144       NA    1
plot(smk_ranks)

(smk_rankprobs <- posterior_rank_probs(smkfit, lower_better = FALSE))
#>                           p_rank[1] p_rank[2] p_rank[3] p_rank[4]
#> d[No intervention]             0.00      0.00      0.10      0.90
#> d[Group counselling]           0.70      0.23      0.07      0.01
#> d[Individual counselling]      0.24      0.59      0.17      0.00
#> d[Self-help]                   0.07      0.18      0.66      0.09
plot(smk_rankprobs)

(smk_cumrankprobs <- posterior_rank_probs(smkfit, lower_better = FALSE, cumulative = TRUE))
#>                           p_rank[1] p_rank[2] p_rank[3] p_rank[4]
#> d[No intervention]             0.00      0.00      0.10         1
#> d[Group counselling]           0.70      0.92      0.99         1
#> d[Individual counselling]      0.24      0.83      1.00         1
#> d[Self-help]                   0.07      0.24      0.91         1
plot(smk_cumrankprobs)