This vignette describes the analysis of 10 trials comparing reduced
fat diets to control (non-reduced fat diets) for preventing mortality
-(Hooper et al. 2000; Dias et al. 2011). The data are
+(Hooper et al.
+2000; Dias et al. 2011). The data are
available in this package as dietary_fat:
-
Setting up the network
We begin by setting up the network - here just a pairwise
@@ -383,35 +391,35 @@
Setting up the network
(
r) and the person-years at risk (
E) in each
arm, so we use the function
set_agd_arm(). We set “Control”
as the reference treatment.
-
diet_net <- set_agd_arm(dietary_fat,
- study = studyc,
- trt = trtc,
- r = r,
- E = E,
- trt_ref = "Control",
- sample_size = n)
-diet_net
-#> A network with 10 AgD studies (arm-based).
-#>
-#> ------------------------------------------------------- AgD studies (arm-based) ----
-#> Study Treatment arms
-#> DART 2: Control | Reduced Fat
-#> London Corn/Olive 3: Control | Reduced Fat | Reduced Fat
-#> London Low Fat 2: Control | Reduced Fat
-#> Minnesota Coronary 2: Control | Reduced Fat
-#> MRC Soya 2: Control | Reduced Fat
-#> Oslo Diet-Heart 2: Control | Reduced Fat
-#> STARS 2: Control | Reduced Fat
-#> Sydney Diet-Heart 2: Control | Reduced Fat
-#> Veterans Administration 2: Control | Reduced Fat
-#> Veterans Diet & Skin CA 2: Control | Reduced Fat
-#>
-#> Outcome type: rate
-#> ------------------------------------------------------------------------------------
-#> Total number of treatments: 2
-#> Total number of studies: 10
-#> Reference treatment is: Control
-#> Network is connected
+
diet_net <- set_agd_arm(dietary_fat,
+ study = studyc,
+ trt = trtc,
+ r = r,
+ E = E,
+ trt_ref = "Control",
+ sample_size = n)
+diet_net
+#> A network with 10 AgD studies (arm-based).
+#>
+#> ------------------------------------------------------- AgD studies (arm-based) ----
+#> Study Treatment arms
+#> DART 2: Control | Reduced Fat
+#> London Corn/Olive 3: Control | Reduced Fat | Reduced Fat
+#> London Low Fat 2: Control | Reduced Fat
+#> Minnesota Coronary 2: Control | Reduced Fat
+#> MRC Soya 2: Control | Reduced Fat
+#> Oslo Diet-Heart 2: Control | Reduced Fat
+#> STARS 2: Control | Reduced Fat
+#> Sydney Diet-Heart 2: Control | Reduced Fat
+#> Veterans Administration 2: Control | Reduced Fat
+#> Veterans Diet & Skin CA 2: Control | Reduced Fat
+#>
+#> Outcome type: rate
+#> ------------------------------------------------------------------------------------
+#> Total number of treatments: 2
+#> Total number of studies: 10
+#> Reference treatment is: Control
+#> Network is connected
We also specify the optional sample_size argument,
although it is not strictly necessary here. In this case
sample_size would only be required to produce a network
@@ -433,41 +441,41 @@
Fixed effect meta-analysis
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.
+
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))
+
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.11 -0.04 -0.01 0.03 0.10 3448 1
-#> lp__ 5325.54 0.06 2.28 5320.39 5324.19 5325.82 5327.18 5329.17 1560 1
-#>
-#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:56:03 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).
+
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.11 -0.05 -0.01 0.03 0.09 3426 1
+#> lp__ 5386.35 0.06 2.36 5380.89 5384.99 5386.66 5388.07 5389.89 1528 1
+#>
+#> Samples were drawn using NUTS(diag_e) at Fri Apr 19 13:11:58 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(diet_fit_FE, pars = c("d", "mu"))
+
# 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)
-
+
plot_prior_posterior(diet_fit_FE)
+
Random effects meta-analysis
@@ -479,169 +487,169 @@
Random effects meta-analysis
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.
+
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))
+
diet_fit_RE <- nma(diet_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:
-
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.20 -0.07 -0.02 0.03 0.15 2045 1
-#> lp__ 5340.70 0.13 3.95 5332.26 5338.20 5340.96 5343.47 5347.61 971 1
-#> tau 0.13 0.00 0.11 0.01 0.05 0.10 0.18 0.41 1172 1
-#>
-#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:56:14 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).
+
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.01 0.00 0.09 -0.20 -0.06 -0.01 0.04 0.18 1222 1
+#> lp__ 5379.44 0.12 3.85 5371.02 5377.13 5379.69 5382.13 5386.15 1062 1
+#> tau 0.14 0.00 0.13 0.00 0.05 0.11 0.19 0.47 865 1
+#>
+#> Samples were drawn using NUTS(diag_e) at Fri Apr 19 13:13:08 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(diet_fit_RE, pars = c("d", "mu", "delta"))
+
# 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"))
-
+
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.2 (on 21 data points)
-#> pD: 11
-#> DIC: 33.2
-
(dic_RE <- dic(diet_fit_RE))
-#> Residual deviance: 21.4 (on 21 data points)
-#> pD: 13.6
-#> DIC: 35
+
(dic_FE <- dic(diet_fit_FE))
+#> Residual deviance: 22.1 (on 21 data points)
+#> pD: 10.9
+#> DIC: 33
+
(dic_RE <- dic(diet_fit_RE))
+#> Residual deviance: 21.2 (on 21 data points)
+#> pD: 13.6
+#> DIC: 34.9
Both 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.
-
-
-
-
+
+
+
+
Further results
-
Dias et al. (2011) produce absolute predictions of
-the mortality rates on reduced fat and control diets, assuming a Normal
+
Dias et al. (2011) produce absolute predictions of the
+mortality rates on reduced fat and control diets, assuming a Normal
distribution on the baseline log rate of mortality with mean \(-3\) and precision \(1.77\). We can replicate these results
using the predict() method. The baseline
argument takes a distr() distribution object, with which we
specify the corresponding Normal distribution. We set
type = "response" to produce predicted rates
(type = "link" would produce predicted log rates).
-
pred_FE <- predict(diet_fit_FE,
- baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5),
- type = "response")
-pred_FE
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Control] 0.07 0.06 0.01 0.03 0.05 0.08 0.21 4202 3972 1
-#> pred[Reduced Fat] 0.06 0.06 0.01 0.03 0.05 0.08 0.21 4224 3978 1
-plot(pred_FE)
-
-
pred_RE <- predict(diet_fit_RE,
- baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5),
- type = "response")
-pred_RE
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Control] 0.07 0.06 0.01 0.03 0.05 0.08 0.22 4043 3537 1
-#> pred[Reduced Fat] 0.07 0.06 0.01 0.03 0.05 0.08 0.22 4020 3720 1
-plot(pred_RE)
-
+
pred_FE <- predict(diet_fit_FE,
+ baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5),
+ type = "response")
+pred_FE
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Control] 0.07 0.06 0.01 0.03 0.05 0.08 0.23 3991 3835 1
+#> pred[Reduced Fat] 0.07 0.06 0.01 0.03 0.05 0.08 0.23 3922 3757 1
+plot(pred_FE)
+
+
pred_RE <- predict(diet_fit_RE,
+ baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5),
+ type = "response")
+pred_RE
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Control] 0.06 0.05 0.01 0.03 0.05 0.08 0.21 4081 4077 1
+#> pred[Reduced Fat] 0.06 0.05 0.01 0.03 0.05 0.08 0.21 4057 3912 1
+plot(pred_RE)
+
If the baseline argument is omitted, predicted rates
will be produced for every study in the network based on their estimated
baseline log rate \(\mu_j\):
-
pred_FE_studies <- predict(diet_fit_FE, type = "response")
-pred_FE_studies
-#> ------------------------------------------------------------------- Study: DART ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[DART: Control] 0.06 0 0.05 0.06 0.06 0.06 0.07 5543 2859 1
-#> pred[DART: Reduced Fat] 0.06 0 0.05 0.06 0.06 0.06 0.07 6370 3212 1
-#>
-#> ------------------------------------------------------ Study: London Corn/Olive ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[London Corn/Olive: Control] 0.07 0.03 0.03 0.06 0.07 0.09 0.13 7039 2639 1
-#> pred[London Corn/Olive: Reduced Fat] 0.07 0.02 0.03 0.05 0.07 0.09 0.13 7254 2690 1
-#>
-#> --------------------------------------------------------- Study: London Low Fat ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[London Low Fat: Control] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 8037 2986 1
-#> pred[London Low Fat: Reduced Fat] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 9224 3205 1
-#>
-#> ----------------------------------------------------- Study: Minnesota Coronary ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Minnesota Coronary: Control] 0.05 0 0.05 0.05 0.05 0.06 0.06 4832 3586 1
-#> pred[Minnesota Coronary: Reduced Fat] 0.05 0 0.05 0.05 0.05 0.06 0.06 6395 3547 1
-#>
-#> --------------------------------------------------------------- Study: MRC Soya ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[MRC Soya: Control] 0.04 0.01 0.03 0.04 0.04 0.04 0.05 6442 2770 1
-#> pred[MRC Soya: Reduced Fat] 0.04 0.01 0.03 0.04 0.04 0.04 0.05 7057 2832 1
-#>
-#> -------------------------------------------------------- Study: Oslo Diet-Heart ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Oslo Diet-Heart: Control] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 6028 2675 1
-#> pred[Oslo Diet-Heart: Reduced Fat] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 8013 3279 1
-#>
-#> ------------------------------------------------------------------ Study: STARS ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[STARS: Control] 0.02 0.01 0.01 0.01 0.02 0.03 0.05 7028 2787 1
-#> pred[STARS: Reduced Fat] 0.02 0.01 0.01 0.01 0.02 0.03 0.05 6997 2750 1
-#>
-#> ------------------------------------------------------ Study: Sydney Diet-Heart ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Sydney Diet-Heart: Control] 0.03 0 0.03 0.03 0.03 0.04 0.04 6388 3364 1
-#> pred[Sydney Diet-Heart: Reduced Fat] 0.03 0 0.03 0.03 0.03 0.04 0.04 7326 3336 1
-#>
-#> ------------------------------------------------ Study: Veterans Administration ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS
-#> pred[Veterans Administration: Control] 0.11 0.01 0.1 0.11 0.11 0.12 0.13 4882 3030
-#> pred[Veterans Administration: Reduced Fat] 0.11 0.01 0.1 0.11 0.11 0.12 0.13 7066 3282
-#> Rhat
-#> pred[Veterans Administration: Control] 1
-#> pred[Veterans Administration: Reduced Fat] 1
-#>
-#> ------------------------------------------------ Study: Veterans Diet & Skin CA ----
-#>
-#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS
-#> pred[Veterans Diet & Skin CA: Control] 0.01 0.01 0 0.01 0.01 0.02 0.03 6770 2728
-#> pred[Veterans Diet & Skin CA: Reduced Fat] 0.01 0.01 0 0.01 0.01 0.02 0.03 6895 2652
-#> Rhat
-#> pred[Veterans Diet & Skin CA: Control] 1
-#> pred[Veterans Diet & Skin CA: Reduced Fat] 1
-plot(pred_FE_studies) + ggplot2::facet_grid(Study~., labeller = ggplot2::label_wrap_gen(width = 10))
-
+
pred_FE_studies <- predict(diet_fit_FE, type = "response")
+pred_FE_studies
+#> ------------------------------------------------------------------- Study: DART ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[DART: Control] 0.06 0 0.05 0.06 0.06 0.06 0.07 5601 3098 1
+#> pred[DART: Reduced Fat] 0.06 0 0.05 0.06 0.06 0.06 0.07 7752 3176 1
+#>
+#> ------------------------------------------------------ Study: London Corn/Olive ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[London Corn/Olive: Control] 0.07 0.02 0.03 0.06 0.07 0.09 0.13 7183 2988 1.01
+#> pred[London Corn/Olive: Reduced Fat] 0.07 0.02 0.03 0.06 0.07 0.09 0.13 7121 3035 1.01
+#>
+#> --------------------------------------------------------- Study: London Low Fat ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[London Low Fat: Control] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 6555 2868 1
+#> pred[London Low Fat: Reduced Fat] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 7683 3049 1
+#>
+#> ----------------------------------------------------- Study: Minnesota Coronary ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Minnesota Coronary: Control] 0.05 0 0.05 0.05 0.05 0.06 0.06 4396 3458 1
+#> pred[Minnesota Coronary: Reduced Fat] 0.05 0 0.05 0.05 0.05 0.06 0.06 8137 3772 1
+#>
+#> --------------------------------------------------------------- Study: MRC Soya ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[MRC Soya: Control] 0.04 0.01 0.03 0.04 0.04 0.04 0.05 6224 3055 1
+#> pred[MRC Soya: Reduced Fat] 0.04 0.01 0.03 0.04 0.04 0.04 0.05 7525 3247 1
+#>
+#> -------------------------------------------------------- Study: Oslo Diet-Heart ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Oslo Diet-Heart: Control] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 6372 3078 1
+#> pred[Oslo Diet-Heart: Reduced Fat] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 7681 3100 1
+#>
+#> ------------------------------------------------------------------ Study: STARS ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[STARS: Control] 0.02 0.01 0.01 0.01 0.02 0.03 0.05 7184 2361 1
+#> pred[STARS: Reduced Fat] 0.02 0.01 0.01 0.01 0.02 0.03 0.05 7356 2306 1
+#>
+#> ------------------------------------------------------ Study: Sydney Diet-Heart ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Sydney Diet-Heart: Control] 0.03 0 0.03 0.03 0.03 0.04 0.04 6935 3061 1
+#> pred[Sydney Diet-Heart: Reduced Fat] 0.03 0 0.03 0.03 0.03 0.04 0.04 7792 3483 1
+#>
+#> ------------------------------------------------ Study: Veterans Administration ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS
+#> pred[Veterans Administration: Control] 0.11 0.01 0.1 0.11 0.11 0.12 0.13 5552 3154
+#> pred[Veterans Administration: Reduced Fat] 0.11 0.01 0.1 0.11 0.11 0.12 0.12 7289 3471
+#> Rhat
+#> pred[Veterans Administration: Control] 1
+#> pred[Veterans Administration: Reduced Fat] 1
+#>
+#> ------------------------------------------------ Study: Veterans Diet & Skin CA ----
+#>
+#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS
+#> pred[Veterans Diet & Skin CA: Control] 0.01 0.01 0 0.01 0.01 0.02 0.03 6465 2903
+#> pred[Veterans Diet & Skin CA: Reduced Fat] 0.01 0.01 0 0.01 0.01 0.02 0.03 6548 2765
+#> Rhat
+#> pred[Veterans Diet & Skin CA: Control] 1
+#> pred[Veterans Diet & Skin CA: Reduced Fat] 1
+plot(pred_FE_studies) + ggplot2::facet_grid(Study~., labeller = ggplot2::label_wrap_gen(width = 10))
+
