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. By
default, this will use a Binomial likelihood and a logit link function,
auto-detected from the data.
blocker_fit_FE <- nma(blocker_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))Basic parameter summaries are given by the print()
method:
blocker_fit_FE
#> 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[Beta Blocker] -0.26 0.00 0.05 -0.36 -0.30 -0.26 -0.23 -0.16 3309 1
#> lp__ -6087.71 0.09 3.55 -6095.39 -6089.78 -6087.27 -6085.21 -6081.71 1531 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:55:39 2024.
#> 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(blocker_fit_FE, pars = c("d", "mu"))The prior and posterior distributions can be compared visually using
the plot_prior_posterior() function:
plot_prior_posterior(blocker_fit_FE, prior = "trt")
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
blocker_fit_RE <- nma(blocker_net,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5))Basic parameter summaries are given by the print()
method:
blocker_fit_RE
#> 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 Rhat
#> d[Beta Blocker] -0.25 0.00 0.06 -0.37 -0.29 -0.25 -0.21 -0.12 3506 1
#> lp__ -6100.48 0.17 5.56 -6112.03 -6104.09 -6100.14 -6096.65 -6090.33 1082 1
#> tau 0.13 0.00 0.08 0.01 0.07 0.12 0.18 0.30 1140 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:55:50 2024.
#> 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(blocker_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(blocker_fit_RE, prior = c("trt", "het"))
Model comparison
Model fit can be checked using the dic() function:
(dic_FE <- dic(blocker_fit_FE))
#> Residual deviance: 47.3 (on 44 data points)
#> pD: 23.7
#> DIC: 71(dic_RE <- dic(blocker_fit_RE))
#> Residual deviance: 42 (on 44 data points)
#> pD: 28.1
#> DIC: 70.1The residual deviance is lower under the RE model, which is to be expected as this model is more flexible. However, this comes with an increased effective number of parameters (note the increase in \(p_D\)). As a result, the DIC of both models is very similar and the FE model may be preferred for parsimony.
We can also examine the residual deviance contributions with the
corresponding plot() method.
plot(dic_FE)
plot(dic_RE)
There are a number of points which are not very well fit by the FE model, having posterior mean residual deviance contributions greater than 1. Study 14 is a particularly poor fit under the FE model, but its residual deviance is reduced (although still high) under the RE model. The evidence should be given further careful examination, and consideration given to other issues such as the potential for effect-modifying covariates (Dias et al. 2011).
