Example: Statins for cholesterol lowering

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

This vignette describes the analysis of 19 trials comparing statins to placebo or usual care (Dias et al. 2011). The data are available in this package as statins:

head(statins)
#>   studyn    studyc trtn    trtc prevention   r    n
#> 1      1        4S    1 Placebo  Secondary 256 2223
#> 2      1        4S    2  Statin  Secondary 182 2221
#> 3      2 Bestehorn    1 Placebo  Secondary   4  125
#> 4      2 Bestehorn    2  Statin  Secondary   1  129
#> 5      3     Brown    1 Placebo  Secondary   0   52
#> 6      3     Brown    2  Statin  Secondary   1   94

Dias et al. (2011) used these data to demonstrate meta-regression models adjusting for the binary covariate prevention (primary or secondary prevention), which we recreate here.

Setting up the network

We have data giving the number of deaths (r) out of the total (n) in each arm, so we use the function set_agd_arm() to set up the network. We set placebo as the network reference treatment.

statin_net <- set_agd_arm(statins, 
                          study = studyc,
                          trt = trtc,
                          r = r, 
                          n = n,
                          trt_ref = "Placebo")
statin_net
#> A network with 19 AgD studies (arm-based).
#> 
#> ------------------------------------------------------- AgD studies (arm-based) ---- 
#>  Study     Treatment arms     
#>  4S        2: Placebo | Statin
#>  Bestehorn 2: Placebo | Statin
#>  Brown     2: Placebo | Statin
#>  CCAIT     2: Placebo | Statin
#>  Downs     2: Placebo | Statin
#>  EXCEL     2: Placebo | Statin
#>  Furberg   2: Placebo | Statin
#>  Haskell   2: Placebo | Statin
#>  Jones     2: Placebo | Statin
#>  KAPS      2: Placebo | Statin
#>  ... plus 9 more studies
#> 
#>  Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 2
#> Total number of studies: 19
#> Reference treatment is: Placebo
#> Network is connected

The prevention variable in the statins data frame will automatically be available to use in a meta-regression model.

Meta-analysis models

We fit fixed effect (FE) and random effects (RE) models, with a meta-regression on the binary covariate prevention.

Fixed effect meta-regression

We start by fitting a FE model. We use \(\mathrm{N}(0, 100^2)\) prior distributions for the treatment effect \(d_\mathrm{Statin}\), study-specific intercepts \(\mu_j\), and regression coefficient \(\beta\). 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 with the nma() function, with a fixed effect model specified by trt_effects = "fixed". The regression formula ~ .trt:prevention means that interaction of primary/secondary prevention with treatment will be included; the .trt special variable indicates treatment, and prevention is in the original data set.

statin_fit_FE <- nma(statin_net, 
                     trt_effects = "fixed",
                     regression = ~.trt:prevention,
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_reg = 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.

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

statin_fit_FE
#> A fixed effects NMA with a binomial likelihood (logit link).
#> Regression model: ~.trt:prevention.
#> 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%
#> beta[.trtStatin:preventionSecondary]    -0.21    0.00 0.11    -0.43    -0.28    -0.21    -0.14
#> d[Statin]                               -0.10    0.00 0.10    -0.29    -0.17    -0.10    -0.04
#> lp__                                 -7246.59    0.08 3.28 -7253.97 -7248.57 -7246.24 -7244.23
#>                                         97.5% n_eff Rhat
#> beta[.trtStatin:preventionSecondary]     0.00  2212    1
#> d[Statin]                                0.08  2164    1
#> lp__                                 -7241.06  1738    1
#> 
#> Samples were drawn using NUTS(diag_e) at Thu Feb 24 09:07:41 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(statin_fit_FE, pars = c("d", "beta", "mu"))

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

plot_prior_posterior(statin_fit_FE, prior = c("trt", "reg"))

Random effects meta-regression

We now fit a RE model. We use \(\mathrm{N}(0, 100^2)\) prior distributions for the treatment effect \(d_\mathrm{Statin}\), study-specific intercepts \(\mu_j\), and regression coefficient \(\beta\). We use a \(\textrm{half-N}(0, 5^2)\) prior distribution 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.

Again, the model is fitted with the nma() function, now with trt_effects = "random". We increase adapt_delta to 0.99 to remove a small number of divergent transition errors (the default for RE models is set to 0.95).

statin_fit_RE <- nma(statin_net, 
                     trt_effects = "random",
                     regression = ~.trt:prevention,
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_reg = normal(scale = 100),
                     prior_het = half_normal(scale = 5),
                     adapt_delta = 0.99)
#> 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.

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

statin_fit_RE
#> A random effects NMA with a binomial likelihood (logit link).
#> Regression model: ~.trt:prevention.
#> 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%
#> beta[.trtStatin:preventionSecondary]    -0.30    0.01 0.27    -0.90    -0.43    -0.28    -0.14
#> d[Statin]                               -0.07    0.01 0.21    -0.48    -0.18    -0.08     0.04
#> lp__                                 -7255.76    0.16 5.25 -7266.54 -7259.18 -7255.58 -7252.10
#> tau                                      0.25    0.01 0.21     0.01     0.09     0.20     0.34
#>                                         97.5% n_eff Rhat
#> beta[.trtStatin:preventionSecondary]     0.22  1499    1
#> d[Statin]                                0.39  1452    1
#> lp__                                 -7246.28  1032    1
#> tau                                      0.79   737    1
#> 
#> Samples were drawn using NUTS(diag_e) at Thu Feb 24 09:08: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(statin_fit_RE, pars = c("d", "beta", "mu", "delta"))

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

plot_prior_posterior(statin_fit_RE, prior = c("trt", "reg", "het"))

Model fit and comparison

Model fit can be checked using the dic() function:

(statin_dic_FE <- dic(statin_fit_FE))
#> Residual deviance: 45.8 (on 38 data points)
#>                pD: 21.5
#>               DIC: 67.3
(statin_dic_RE <- dic(statin_fit_RE))
#> Residual deviance: 42.5 (on 38 data points)
#>                pD: 25.1
#>               DIC: 67.6

The DIC is very similar between FE and RE models, so we might choose the FE model based on parsimony. The residual deviance statistics are larger than the number of data points, suggesting that some data points are not fit well.

We can also examine the residual deviance contributions with the corresponding plot() method.

plot(statin_dic_FE)

plot(statin_dic_RE)

There are a number of studies which are not fit well under either model, having posterior mean residual deviance contributions greater than 1, and should be investigated to see if there are further substantive differences between studies.

Further results

We can produce estimates of the relative effect of statins vs. placebo for either primary or secondary prevention, using the relative_effects() function. The newdata argument specifies a data frame containing the levels of the covariate prevention that we are interested in, and the study argument is used to specify a column of newdata for an informative label.

statin_releff_FE <- relative_effects(statin_fit_FE,
                                     newdata = data.frame(prevention = c("Primary", "Secondary")),
                                     study = prevention)

statin_releff_FE
#> ---------------------------------------------------------------- Study: Primary ---- 
#> 
#> Covariate values:
#>  prevention
#>     Primary
#> 
#>                    mean  sd  2.5%   25%  50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Primary: Statin] -0.1 0.1 -0.29 -0.17 -0.1 -0.04  0.08     2200     2354    1
#> 
#> -------------------------------------------------------------- Study: Secondary ---- 
#> 
#> Covariate values:
#>  prevention
#>   Secondary
#> 
#>                       mean   sd  2.5%   25%   50%   75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Secondary: Statin] -0.31 0.05 -0.42 -0.35 -0.31 -0.28 -0.21     4047     3614    1

The plot() method may be used to visually compare these estimates:

plot(statin_releff_FE, 
     ref_line = 0)

Model parameters may be plotted with the corresponding plot() method:

plot(statin_fit_FE, 
     pars = "beta", 
     ref_line = 0,
     stat = "halfeye")

Whilst the 95% Credible Interval includes zero, there is a suggestion that statins are more effective for secondary prevention.

References

Dias, S., A. J. Sutton, N. J. Welton, and A. E. Ades. 2011. “NICE DSU Technical Support Document 3: Heterogeneity: Subgroups, Meta-Regression, Bias and Bias-Adjustment.” National Institute for Health and Care Excellence. https://nicedsu.sites.sheffield.ac.uk.