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 Poisson likelihood with a log link function,
auto-detected from the data.
diet_fit_FE <- nma(diet_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))Basic parameter summaries are given by the print()
method:
diet_fit_FE
#> A fixed effects NMA with a poisson likelihood (log link).
#> Inference for Stan model: poisson.
#> 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[Reduced Fat] -0.01 0.00 0.05 -0.12 -0.04 -0.01 0.03 0.10 3295 1
#> lp__ 5386.23 0.06 2.41 5380.54 5384.77 5386.60 5388.01 5390.01 1699 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue May 23 11:34:52 2023.
#> 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(diet_fit_FE, pars = c("d", "mu"))The prior and posterior distributions can be compared visually using
the plot_prior_posterior() function:
plot_prior_posterior(diet_fit_FE)
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
diet_fit_RE <- nma(diet_net,
trt_effects = "random",
prior_intercept = normal(scale = 10),
prior_trt = normal(scale = 10),
prior_het = half_normal(scale = 5))Basic parameter summaries are given by the print()
method:
diet_fit_RE
#> A random effects NMA with a poisson likelihood (log link).
#> Inference for Stan model: poisson.
#> 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[Reduced Fat] -0.02 0.00 0.09 -0.21 -0.07 -0.02 0.03 0.17 1269 1
#> lp__ 5378.69 0.11 3.83 5370.57 5376.27 5378.85 5381.35 5385.60 1216 1
#> tau 0.13 0.00 0.11 0.00 0.05 0.10 0.18 0.42 1008 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue May 23 11:35:05 2023.
#> 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(diet_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(diet_fit_RE, prior = c("trt", "het"))
Model comparison
Model fit can be checked using the dic() function:
(dic_FE <- dic(diet_fit_FE))
#> Residual deviance: 22.4 (on 21 data points)
#> pD: 11.1
#> DIC: 33.5(dic_RE <- dic(diet_fit_RE))
#> Residual deviance: 21.4 (on 21 data points)
#> pD: 13.5
#> DIC: 34.9Both models appear to fit the data well, as the residual deviance is close to the number of data points. The DIC is very similar between models, so 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)
