library(multinma)
options(mc.cores = parallel::detectCores())
This vignette describes the analysis of data on the mean off-time
reduction in patients given dopamine agonists as adjunct therapy in
Parkinson’s disease, in a network of 7 trials of 4 active drugs plus
placebo (Dias et al. 2011). The data are
available in this package as parkinsons
:
head(parkinsons)
#> studyn trtn y se n diff se_diff
#> 1 1 1 -1.22 0.504 54 NA 0.504
#> 2 1 3 -1.53 0.439 95 -0.31 0.668
#> 3 2 1 -0.70 0.282 172 NA 0.282
#> 4 2 2 -2.40 0.258 173 -1.70 0.382
#> 5 3 1 -0.30 0.505 76 NA 0.505
#> 6 3 2 -2.60 0.510 71 -2.30 0.718
We consider analysing these data in three separate ways:
y
and corresponding
standard errors se
);diff
and
corresponding standard errors se_diff
);Note: In this case, with Normal likelihoods for both arms and contrasts, we will see that the three analyses give identical results. In general, unless the arm-based likelihood is Normal, results from a model using a contrast-based likelihood will not exactly match those from a model using an arm-based likelihood, since the contrast-based Normal likelihood is only an approximation. Similarity of results depends on the suitability of the Normal approximation, which may not always be appropriate - e.g. with a small number of events or small sample size for a binary outcome. The use of an arm-based likelihood (sometimes called an “exact” likelihood) is therefore preferable where possible in general.
We begin with an analysis of the arm-based data - means and standard errors.
We have arm-level continuous data giving the mean off-time reduction
(y
) and standard error (se
) in each arm. We
use the function set_agd_arm()
to set up the network.
<- set_agd_arm(parkinsons,
arm_net study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n)
arm_net#> A network with 7 AgD studies (arm-based).
#>
#> ------------------------------------------------------- AgD studies (arm-based) ----
#> Study Treatment arms
#> 1 2: 1 | 3
#> 2 2: 1 | 2
#> 3 3: 4 | 1 | 2
#> 4 2: 4 | 3
#> 5 2: 4 | 3
#> 6 2: 4 | 5
#> 7 2: 4 | 5
#>
#> Outcome type: continuous
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 5
#> Total number of studies: 7
#> Reference treatment is: 4
#> Network is connected
We let treatment 4 be set by default as the network reference
treatment, since this results in considerably improved sampling
efficiency over choosing treatment 1 as the network reference. The
sample_size
argument is optional, but enables the nodes to
be weighted by sample size in the network plot.
Plot the network structure.
plot(arm_net, weight_edges = TRUE, weight_nodes = TRUE)
We fit both fixed effect (FE) and random effects (RE) models.
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.
<- nma(arm_net,
arm_fit_FE trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10))
#> Note: Setting "4" as the network reference treatment.
Basic parameter summaries are given by the print()
method:
arm_fit_FE#> A fixed effects NMA with a normal likelihood (identity link).
#> Inference for Stan model: normal.
#> 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[1] 0.52 0.01 0.49 -0.44 0.19 0.51 0.84 1.49 1455 1
#> d[2] -1.29 0.01 0.53 -2.36 -1.63 -1.30 -0.93 -0.24 1447 1
#> d[3] 0.05 0.01 0.32 -0.59 -0.17 0.05 0.26 0.68 2084 1
#> d[5] -0.30 0.00 0.21 -0.70 -0.45 -0.31 -0.17 0.11 3020 1
#> lp__ -58.31 0.06 2.45 -64.18 -59.72 -57.98 -56.53 -54.61 1593 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:57:15 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(arm_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using
the plot_prior_posterior()
function:
plot_prior_posterior(arm_fit_FE)
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
<- nma(arm_net,
arm_fit_RE seed = 379394727,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99)
#> Note: Setting "4" as the network reference treatment.
#> Warning: There were 2 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
We do see a small number of divergent transition errors, which cannot
simply be removed by increasing the value of the
adapt_delta
argument (by default set to 0.95
for RE models). To diagnose, we use the pairs()
method to
investigate where in the posterior distribution these divergences are
happening (indicated by red crosses):
pairs(arm_fit_RE, pars = c("mu[4]", "d[3]", "delta[4: 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation. In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed. We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print()
method:
arm_fit_RE#> A random effects NMA with a normal likelihood (identity link).
#> Inference for Stan model: normal.
#> 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[1] 0.52 0.02 0.62 -0.63 0.14 0.51 0.90 1.78 1522 1
#> d[2] -1.32 0.02 0.68 -2.68 -1.74 -1.32 -0.91 -0.01 1633 1
#> d[3] 0.03 0.01 0.48 -0.89 -0.24 0.04 0.31 0.94 2222 1
#> d[5] -0.31 0.01 0.43 -1.15 -0.50 -0.31 -0.10 0.49 1629 1
#> lp__ -76.05 0.10 3.59 -84.03 -78.19 -75.72 -73.52 -70.04 1358 1
#> tau 0.38 0.01 0.38 0.01 0.12 0.27 0.51 1.37 723 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:57:26 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(arm_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(arm_fit_RE)
Model fit can be checked using the dic()
function:
<- dic(arm_fit_FE))
(arm_dic_FE #> Residual deviance: 13.5 (on 15 data points)
#> pD: 11.2
#> DIC: 24.7
<- dic(arm_fit_RE))
(arm_dic_RE #> Residual deviance: 13.5 (on 15 data points)
#> pD: 12.3
#> DIC: 25.8
Both models fit the data well, having posterior mean residual deviance close to the number of data points. The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the
corresponding plot()
method.
plot(arm_dic_FE)
plot(arm_dic_RE)
For comparison with Dias et al. (2011), we can produce
relative effects against placebo using the
relative_effects()
function with
trt_ref = 1
:
<- relative_effects(arm_fit_FE, trt_ref = 1))
(arm_releff_FE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[4] -0.52 0.49 -1.49 -0.84 -0.51 -0.19 0.44 1462 2068 1
#> d[2] -1.81 0.33 -2.48 -2.03 -1.81 -1.58 -1.16 5387 3162 1
#> d[3] -0.47 0.50 -1.44 -0.82 -0.47 -0.13 0.50 2184 2601 1
#> d[5] -0.82 0.53 -1.87 -1.17 -0.82 -0.47 0.21 1571 2072 1
plot(arm_releff_FE, ref_line = 0)
<- relative_effects(arm_fit_RE, trt_ref = 1))
(arm_releff_RE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[4] -0.52 0.62 -1.78 -0.90 -0.51 -0.14 0.63 1551 1917 1
#> d[2] -1.85 0.50 -2.86 -2.12 -1.84 -1.55 -0.91 4153 2853 1
#> d[3] -0.49 0.62 -1.69 -0.87 -0.49 -0.11 0.68 3010 2412 1
#> d[5] -0.83 0.76 -2.38 -1.25 -0.81 -0.38 0.58 1704 1902 1
plot(arm_releff_RE, ref_line = 0)
Following Dias et al. (2011), we produce absolute predictions
of the mean off-time reduction on each treatment assuming a Normal
distribution for the outcomes on treatment 1 (placebo) with mean \(-0.73\) and precision \(21\). We use the predict()
method, where the baseline
argument takes a
distr()
distribution object with which we specify the
corresponding Normal distribution, and we specify
trt_ref = 1
to indicate that the baseline distribution
corresponds to treatment 1. (Strictly speaking,
type = "response"
is unnecessary here, since the identity
link function was used.)
<- predict(arm_fit_FE,
arm_pred_FE baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
trt_ref = 1)
arm_pred_FE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[4] -1.24 0.53 -2.29 -1.60 -1.25 -0.88 -0.20 1570 2220 1
#> pred[1] -0.73 0.22 -1.16 -0.88 -0.73 -0.58 -0.29 3910 3936 1
#> pred[2] -2.53 0.40 -3.31 -2.78 -2.53 -2.27 -1.75 4820 3588 1
#> pred[3] -1.20 0.54 -2.25 -1.56 -1.20 -0.81 -0.15 2260 2867 1
#> pred[5] -1.55 0.57 -2.69 -1.92 -1.55 -1.18 -0.45 1664 2609 1
plot(arm_pred_FE)
<- predict(arm_fit_RE,
arm_pred_RE baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
trt_ref = 1)
arm_pred_RE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[4] -1.25 0.65 -2.60 -1.66 -1.24 -0.82 -0.06 1653 1893 1
#> pred[1] -0.73 0.22 -1.16 -0.87 -0.73 -0.58 -0.30 3976 3867 1
#> pred[2] -2.58 0.55 -3.64 -2.89 -2.57 -2.25 -1.56 4242 2675 1
#> pred[3] -1.22 0.66 -2.50 -1.64 -1.21 -0.80 0.03 3146 2589 1
#> pred[5] -1.56 0.79 -3.12 -2.00 -1.53 -1.09 -0.14 1760 1765 1
plot(arm_pred_RE)
If the baseline
argument is omitted, predictions of mean
off-time reduction will be produced for every study in the network based
on their estimated baseline response \(\mu_j\):
<- predict(arm_fit_FE, type = "response")
arm_pred_FE_studies
arm_pred_FE_studies#> ---------------------------------------------------------------------- Study: 1 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[1: 4] -1.64 0.46 -2.59 -1.94 -1.64 -1.35 -0.76 1940 2553 1
#> pred[1: 1] -1.13 0.43 -1.97 -1.42 -1.12 -0.84 -0.29 3545 2759 1
#> pred[1: 2] -2.93 0.51 -3.95 -3.29 -2.93 -2.59 -1.96 3165 2850 1
#> pred[1: 3] -1.60 0.40 -2.40 -1.86 -1.59 -1.33 -0.80 3429 3084 1
#> pred[1: 5] -1.95 0.50 -2.96 -2.28 -1.95 -1.61 -0.98 1982 2436 1
#>
#> ---------------------------------------------------------------------- Study: 2 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[2: 4] -1.16 0.52 -2.17 -1.50 -1.15 -0.82 -0.13 1447 2010 1
#> pred[2: 1] -0.64 0.26 -1.14 -0.81 -0.64 -0.47 -0.13 4796 3864 1
#> pred[2: 2] -2.45 0.25 -2.92 -2.62 -2.44 -2.28 -1.97 4864 3265 1
#> pred[2: 3] -1.11 0.54 -2.18 -1.47 -1.11 -0.76 -0.06 2060 2552 1
#> pred[2: 5] -1.46 0.56 -2.57 -1.82 -1.46 -1.10 -0.37 1557 2098 1
#>
#> ---------------------------------------------------------------------- Study: 3 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[3: 4] -1.11 0.42 -1.95 -1.39 -1.10 -0.83 -0.29 1714 2578 1
#> pred[3: 1] -0.59 0.38 -1.33 -0.85 -0.60 -0.34 0.13 4890 3255 1
#> pred[3: 2] -2.40 0.40 -3.19 -2.67 -2.39 -2.13 -1.62 4261 2897 1
#> pred[3: 3] -1.06 0.48 -2.01 -1.38 -1.06 -0.75 -0.12 2691 3161 1
#> pred[3: 5] -1.41 0.47 -2.33 -1.72 -1.42 -1.09 -0.49 1803 2763 1
#>
#> ---------------------------------------------------------------------- Study: 4 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[4: 4] -0.40 0.30 -0.99 -0.60 -0.39 -0.20 0.18 2265 2762 1
#> pred[4: 1] 0.12 0.51 -0.87 -0.24 0.11 0.48 1.09 2127 2799 1
#> pred[4: 2] -1.69 0.56 -2.77 -2.06 -1.69 -1.31 -0.59 1962 2477 1
#> pred[4: 3] -0.35 0.24 -0.83 -0.52 -0.36 -0.18 0.12 4953 3214 1
#> pred[4: 5] -0.70 0.36 -1.40 -0.94 -0.70 -0.46 -0.01 2477 2779 1
#>
#> ---------------------------------------------------------------------- Study: 5 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[5: 4] -0.56 0.35 -1.25 -0.79 -0.54 -0.33 0.12 2514 2543 1
#> pred[5: 1] -0.04 0.55 -1.15 -0.42 -0.03 0.34 1.02 2143 2665 1
#> pred[5: 2] -1.85 0.60 -3.01 -2.23 -1.84 -1.46 -0.69 1966 2181 1
#> pred[5: 3] -0.51 0.29 -1.08 -0.71 -0.51 -0.32 0.06 4987 2803 1
#> pred[5: 5] -0.86 0.40 -1.67 -1.13 -0.86 -0.59 -0.07 2686 2693 1
#>
#> ---------------------------------------------------------------------- Study: 6 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[6: 4] -2.19 0.17 -2.53 -2.31 -2.19 -2.08 -1.86 2909 2629 1
#> pred[6: 1] -1.68 0.51 -2.71 -2.02 -1.68 -1.34 -0.65 1584 2059 1
#> pred[6: 2] -3.48 0.56 -4.58 -3.84 -3.49 -3.11 -2.37 1566 2208 1
#> pred[6: 3] -2.15 0.37 -2.87 -2.39 -2.15 -1.90 -1.43 2340 2707 1
#> pred[6: 5] -2.50 0.17 -2.82 -2.61 -2.50 -2.39 -2.17 5579 3339 1
#>
#> ---------------------------------------------------------------------- Study: 7 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[7: 4] -1.79 0.18 -2.13 -1.91 -1.79 -1.67 -1.45 3591 3141 1
#> pred[7: 1] -1.28 0.52 -2.27 -1.63 -1.28 -0.93 -0.21 1579 2161 1
#> pred[7: 2] -3.08 0.56 -4.17 -3.47 -3.09 -2.71 -1.95 1546 2169 1
#> pred[7: 3] -1.75 0.37 -2.48 -2.00 -1.75 -1.50 -1.05 2259 2571 1
#> pred[7: 5] -2.10 0.21 -2.49 -2.24 -2.10 -1.96 -1.69 4670 3172 1
plot(arm_pred_FE_studies)
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
<- posterior_ranks(arm_fit_FE))
(arm_ranks #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[4] 3.51 0.72 2 3 3 4 5 2091 NA 1
#> rank[1] 4.62 0.79 2 5 5 5 5 1999 NA 1
#> rank[2] 1.06 0.31 1 1 1 1 2 2051 2153 1
#> rank[3] 3.55 0.92 2 3 4 4 5 2948 NA 1
#> rank[5] 2.26 0.67 1 2 2 2 4 2493 2702 1
plot(arm_ranks)
<- posterior_rank_probs(arm_fit_FE))
(arm_rankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5]
#> d[4] 0.00 0.04 0.49 0.38 0.09
#> d[1] 0.00 0.04 0.07 0.12 0.77
#> d[2] 0.95 0.04 0.01 0.00 0.00
#> d[3] 0.00 0.16 0.26 0.44 0.14
#> d[5] 0.05 0.72 0.17 0.05 0.01
plot(arm_rankprobs)
<- posterior_rank_probs(arm_fit_FE, cumulative = TRUE))
(arm_cumrankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5]
#> d[4] 0.00 0.05 0.53 0.91 1
#> d[1] 0.00 0.04 0.11 0.23 1
#> d[2] 0.95 0.99 1.00 1.00 1
#> d[3] 0.00 0.16 0.42 0.86 1
#> d[5] 0.05 0.76 0.94 0.99 1
plot(arm_cumrankprobs)
We now perform an analysis using the contrast-based data (mean differences and standard errors).
With contrast-level data giving the mean difference in off-time
reduction (diff
) and standard error (se_diff
),
we use the function set_agd_contrast()
to set up the
network.
<- set_agd_contrast(parkinsons,
contr_net study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)
contr_net#> A network with 7 AgD studies (contrast-based).
#>
#> -------------------------------------------------- AgD studies (contrast-based) ----
#> Study Treatment arms
#> 1 2: 1 | 3
#> 2 2: 1 | 2
#> 3 3: 4 | 1 | 2
#> 4 2: 4 | 3
#> 5 2: 4 | 3
#> 6 2: 4 | 5
#> 7 2: 4 | 5
#>
#> Outcome type: continuous
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 5
#> Total number of studies: 7
#> Reference treatment is: 4
#> Network is connected
The sample_size
argument is optional, but enables the
nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(contr_net, weight_edges = TRUE, weight_nodes = TRUE)
We fit both fixed effect (FE) and random effects (RE) models.
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\). 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.
<- nma(contr_net,
contr_fit_FE trt_effects = "fixed",
prior_trt = normal(scale = 100))
#> Note: Setting "4" as the network reference treatment.
Basic parameter summaries are given by the print()
method:
contr_fit_FE#> A fixed effects NMA with a normal likelihood (identity link).
#> Inference for Stan model: normal.
#> 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[1] 0.52 0.01 0.48 -0.40 0.19 0.51 0.85 1.50 2071 1
#> d[2] -1.29 0.01 0.52 -2.34 -1.64 -1.30 -0.94 -0.27 2191 1
#> d[3] 0.04 0.01 0.33 -0.59 -0.19 0.04 0.27 0.69 2688 1
#> d[5] -0.30 0.00 0.21 -0.71 -0.44 -0.31 -0.16 0.12 3076 1
#> lp__ -25.27 0.03 1.41 -28.91 -25.98 -24.95 -24.22 -23.49 1720 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:57:46 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).
The prior and posterior distributions can be compared visually using
the plot_prior_posterior()
function:
plot_prior_posterior(contr_fit_FE)
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
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
<- nma(contr_net,
contr_fit_RE seed = 1150676438,
trt_effects = "random",
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99)
#> Note: Setting "4" as the network reference treatment.
#> Warning: There were 1 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
We do see a small number of divergent transition errors, which cannot
simply be removed by increasing the value of the
adapt_delta
argument (by default set to 0.95
for RE models). To diagnose, we use the pairs()
method to
investigate where in the posterior distribution these divergences are
happening (indicated by red crosses):
pairs(contr_fit_RE, pars = c("d[3]", "delta[4: 4 vs. 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation. In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed. We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print()
method:
contr_fit_RE#> A random effects NMA with a normal likelihood (identity link).
#> Inference for Stan model: normal.
#> 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[1] 0.53 0.01 0.62 -0.64 0.15 0.52 0.90 1.77 2374 1
#> d[2] -1.31 0.01 0.70 -2.64 -1.73 -1.31 -0.91 0.08 2315 1
#> d[3] 0.03 0.01 0.48 -0.89 -0.24 0.03 0.30 0.92 2405 1
#> d[5] -0.30 0.01 0.46 -1.27 -0.51 -0.31 -0.10 0.58 1453 1
#> lp__ -32.87 0.09 2.87 -39.20 -34.57 -32.61 -30.85 -27.99 1067 1
#> tau 0.40 0.02 0.43 0.01 0.12 0.28 0.51 1.51 789 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:57:55 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 relative effects \(\delta_{jk}\) are hidden, but could be
examined by changing the pars
argument:
# Not run
print(contr_fit_RE, pars = c("d", "delta"))
The prior and posterior distributions can be compared visually using
the plot_prior_posterior()
function:
plot_prior_posterior(contr_fit_RE)
Model fit can be checked using the dic()
function:
<- dic(contr_fit_FE))
(contr_dic_FE #> Residual deviance: 6.3 (on 8 data points)
#> pD: 4.1
#> DIC: 10.4
<- dic(contr_fit_RE))
(contr_dic_RE #> Residual deviance: 6.5 (on 8 data points)
#> pD: 5.3
#> DIC: 11.8
Both models fit the data well, having posterior mean residual deviance close to the number of data points. The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the
corresponding plot()
method.
plot(contr_dic_FE)
plot(contr_dic_RE)
For comparison with Dias et al. (2011), we can produce
relative effects against placebo using the
relative_effects()
function with
trt_ref = 1
:
<- relative_effects(contr_fit_FE, trt_ref = 1))
(contr_releff_FE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[4] -0.52 0.48 -1.50 -0.85 -0.51 -0.19 0.40 2079 1902 1
#> d[2] -1.81 0.33 -2.47 -2.03 -1.81 -1.59 -1.18 5161 3033 1
#> d[3] -0.48 0.49 -1.43 -0.81 -0.46 -0.16 0.49 2944 2987 1
#> d[5] -0.82 0.52 -1.85 -1.18 -0.81 -0.46 0.19 2210 2363 1
plot(contr_releff_FE, ref_line = 0)
<- relative_effects(contr_fit_RE, trt_ref = 1))
(contr_releff_RE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[4] -0.53 0.62 -1.77 -0.90 -0.52 -0.15 0.64 2734 1906 1
#> d[2] -1.84 0.51 -2.83 -2.14 -1.84 -1.55 -0.81 3720 2532 1
#> d[3] -0.50 0.63 -1.76 -0.86 -0.50 -0.12 0.72 3618 2669 1
#> d[5] -0.83 0.74 -2.34 -1.26 -0.84 -0.39 0.56 2394 1698 1
plot(contr_releff_RE, ref_line = 0)
Following Dias et al. (2011), we produce absolute predictions
of the mean off-time reduction on each treatment assuming a Normal
distribution for the outcomes on treatment 1 (placebo) with mean \(-0.73\) and precision \(21\). We use the predict()
method, where the baseline
argument takes a
distr()
distribution object with which we specify the
corresponding Normal distribution, and we specify
trt_ref = 1
to indicate that the baseline distribution
corresponds to treatment 1. (Strictly speaking,
type = "response"
is unnecessary here, since the identity
link function was used.)
<- predict(contr_fit_FE,
contr_pred_FE baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
trt_ref = 1)
contr_pred_FE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[4] -1.25 0.53 -2.32 -1.60 -1.24 -0.89 -0.22 2373 2538 1
#> pred[1] -0.73 0.22 -1.16 -0.88 -0.73 -0.58 -0.30 3988 3629 1
#> pred[2] -2.54 0.39 -3.34 -2.80 -2.54 -2.28 -1.77 4817 3183 1
#> pred[3] -1.21 0.54 -2.27 -1.56 -1.20 -0.85 -0.15 3205 3270 1
#> pred[5] -1.55 0.57 -2.66 -1.93 -1.55 -1.17 -0.44 2493 2437 1
plot(contr_pred_FE)
<- predict(contr_fit_RE,
contr_pred_RE baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
trt_ref = 1)
contr_pred_RE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[4] -1.26 0.66 -2.54 -1.67 -1.24 -0.86 -0.01 2879 2089 1
#> pred[1] -0.73 0.22 -1.16 -0.88 -0.73 -0.59 -0.31 3715 3592 1
#> pred[2] -2.57 0.56 -3.66 -2.91 -2.58 -2.24 -1.46 3722 2467 1
#> pred[3] -1.23 0.67 -2.57 -1.63 -1.24 -0.81 0.05 3715 2838 1
#> pred[5] -1.57 0.78 -3.12 -2.02 -1.56 -1.09 -0.09 2359 1834 1
plot(contr_pred_RE)
If the baseline
argument is omitted an error will be
raised, as there are no study baselines estimated in this network.
# Not run
predict(contr_fit_FE, type = "response")
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
<- posterior_ranks(contr_fit_FE))
(contr_ranks #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[4] 3.51 0.73 2 3 3 4 5 2462 NA 1
#> rank[1] 4.64 0.77 2 5 5 5 5 2330 NA 1
#> rank[2] 1.05 0.28 1 1 1 1 2 2564 2522 1
#> rank[3] 3.51 0.92 2 3 4 4 5 3676 NA 1
#> rank[5] 2.29 0.67 1 2 2 2 4 2598 2496 1
plot(contr_ranks)
<- posterior_rank_probs(contr_fit_FE))
(contr_rankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5]
#> d[4] 0.00 0.05 0.48 0.38 0.09
#> d[1] 0.00 0.04 0.06 0.12 0.78
#> d[2] 0.96 0.03 0.01 0.00 0.00
#> d[3] 0.00 0.16 0.27 0.44 0.12
#> d[5] 0.04 0.72 0.18 0.06 0.01
plot(contr_rankprobs)
<- posterior_rank_probs(contr_fit_FE, cumulative = TRUE))
(contr_cumrankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5]
#> d[4] 0.00 0.05 0.53 0.91 1
#> d[1] 0.00 0.04 0.10 0.22 1
#> d[2] 0.96 0.99 1.00 1.00 1
#> d[3] 0.00 0.17 0.44 0.88 1
#> d[5] 0.04 0.76 0.93 0.99 1
plot(contr_cumrankprobs)
We now perform an analysis where some studies contribute arm-based data, and other contribute contrast-based data. Replicating Dias et al. (2011), we consider arm-based data from studies 1-3, and contrast-based data from studies 4-7.
<- parkinsons$studyn
studies <- parkinsons[studies %in% 1:3, ])
(parkinsons_arm #> studyn trtn y se n diff se_diff
#> 1 1 1 -1.22 0.504 54 NA 0.504
#> 2 1 3 -1.53 0.439 95 -0.31 0.668
#> 3 2 1 -0.70 0.282 172 NA 0.282
#> 4 2 2 -2.40 0.258 173 -1.70 0.382
#> 5 3 1 -0.30 0.505 76 NA 0.505
#> 6 3 2 -2.60 0.510 71 -2.30 0.718
#> 7 3 4 -1.20 0.478 81 -0.90 0.695
<- parkinsons[studies %in% 4:7, ])
(parkinsons_contr #> studyn trtn y se n diff se_diff
#> 8 4 3 -0.24 0.265 128 NA 0.265
#> 9 4 4 -0.59 0.354 72 -0.35 0.442
#> 10 5 3 -0.73 0.335 80 NA 0.335
#> 11 5 4 -0.18 0.442 46 0.55 0.555
#> 12 6 4 -2.20 0.197 137 NA 0.197
#> 13 6 5 -2.50 0.190 131 -0.30 0.274
#> 14 7 4 -1.80 0.200 154 NA 0.200
#> 15 7 5 -2.10 0.250 143 -0.30 0.320
We use the functions set_agd_arm()
and
set_agd_contrast()
to set up the respective data sources
within the network, and then combine together with
combine_network()
.
<- set_agd_arm(parkinsons_arm,
mix_arm_net study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n)
<- set_agd_contrast(parkinsons_contr,
mix_contr_net study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)
<- combine_network(mix_arm_net, mix_contr_net)
mix_net
mix_net#> A network with 3 AgD studies (arm-based), and 4 AgD studies (contrast-based).
#>
#> ------------------------------------------------------- AgD studies (arm-based) ----
#> Study Treatment arms
#> 1 2: 1 | 3
#> 2 2: 1 | 2
#> 3 3: 4 | 1 | 2
#>
#> Outcome type: continuous
#> -------------------------------------------------- AgD studies (contrast-based) ----
#> Study Treatment arms
#> 4 2: 4 | 3
#> 5 2: 4 | 3
#> 6 2: 4 | 5
#> 7 2: 4 | 5
#>
#> Outcome type: continuous
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 5
#> Total number of studies: 7
#> Reference treatment is: 4
#> Network is connected
The sample_size
argument is optional, but enables the
nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(mix_net, weight_edges = TRUE, weight_nodes = TRUE)
We fit both fixed effect (FE) and random effects (RE) models.
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.
<- nma(mix_net,
mix_fit_FE trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
#> Note: Setting "4" as the network reference treatment.
Basic parameter summaries are given by the print()
method:
mix_fit_FE#> A fixed effects NMA with a normal likelihood (identity link).
#> Inference for Stan model: normal.
#> 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[1] 0.54 0.01 0.48 -0.43 0.21 0.56 0.87 1.46 1400 1
#> d[2] -1.28 0.01 0.52 -2.29 -1.64 -1.27 -0.93 -0.31 1453 1
#> d[3] 0.05 0.01 0.33 -0.58 -0.17 0.05 0.27 0.72 2193 1
#> d[5] -0.31 0.00 0.21 -0.71 -0.44 -0.31 -0.17 0.11 3102 1
#> lp__ -43.33 0.05 1.88 -47.85 -44.38 -42.97 -41.97 -40.62 1734 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:58:10 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(mix_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using
the plot_prior_posterior()
function:
plot_prior_posterior(mix_fit_FE)
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
<- nma(mix_net,
mix_fit_RE seed = 437219664,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99)
#> Note: Setting "4" as the network reference treatment.
#> Warning: There were 2 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
We do see a small number of divergent transition errors, which cannot
simply be removed by increasing the value of the
adapt_delta
argument (by default set to 0.95
for RE models). To diagnose, we use the pairs()
method to
investigate where in the posterior distribution these divergences are
happening (indicated by red crosses):
pairs(mix_fit_RE, pars = c("d[3]", "delta[4: 4 vs. 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation. In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed. We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print()
method:
mix_fit_RE#> A random effects NMA with a normal likelihood (identity link).
#> Inference for Stan model: normal.
#> 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[1] 0.52 0.02 0.62 -0.68 0.14 0.52 0.90 1.73 1630 1.00
#> d[2] -1.32 0.02 0.69 -2.70 -1.74 -1.31 -0.90 0.01 1869 1.00
#> d[3] 0.01 0.01 0.46 -0.92 -0.26 0.02 0.28 0.88 2335 1.00
#> d[5] -0.28 0.01 0.41 -1.07 -0.49 -0.29 -0.08 0.55 2213 1.00
#> lp__ -51.85 0.11 3.31 -59.22 -53.81 -51.51 -49.50 -46.37 938 1.00
#> tau 0.39 0.01 0.37 0.01 0.13 0.28 0.52 1.35 601 1.01
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:58:21 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(mix_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(mix_fit_RE)
Model fit can be checked using the dic()
function:
<- dic(mix_fit_FE))
(mix_dic_FE #> Residual deviance: 9.3 (on 11 data points)
#> pD: 7
#> DIC: 16.4
<- dic(mix_fit_RE))
(mix_dic_RE #> Residual deviance: 9.5 (on 11 data points)
#> pD: 8.3
#> DIC: 17.8
Both models fit the data well, having posterior mean residual deviance close to the number of data points. The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the
corresponding plot()
method.
plot(mix_dic_FE)
plot(mix_dic_RE)
For comparison with Dias et al. (2011), we can produce
relative effects against placebo using the
relative_effects()
function with
trt_ref = 1
:
<- relative_effects(mix_fit_FE, trt_ref = 1))
(mix_releff_FE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[4] -0.54 0.48 -1.46 -0.87 -0.56 -0.21 0.43 1439 2275 1
#> d[2] -1.82 0.33 -2.46 -2.04 -1.82 -1.59 -1.14 5602 3215 1
#> d[3] -0.49 0.48 -1.43 -0.82 -0.49 -0.16 0.46 2255 2644 1
#> d[5] -0.84 0.52 -1.86 -1.21 -0.85 -0.49 0.17 1559 2448 1
plot(mix_releff_FE, ref_line = 0)
<- relative_effects(mix_fit_RE, trt_ref = 1))
(mix_releff_RE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[4] -0.52 0.62 -1.73 -0.90 -0.52 -0.14 0.68 1690 1682 1
#> d[2] -1.84 0.50 -2.86 -2.12 -1.83 -1.55 -0.86 4089 2312 1
#> d[3] -0.51 0.62 -1.72 -0.89 -0.50 -0.12 0.70 2941 2065 1
#> d[5] -0.80 0.76 -2.31 -1.24 -0.79 -0.38 0.70 1791 1643 1
plot(mix_releff_RE, ref_line = 0)
Following Dias et al. (2011), we produce absolute predictions
of the mean off-time reduction on each treatment assuming a Normal
distribution for the outcomes on treatment 1 (placebo) with mean \(-0.73\) and precision \(21\). We use the predict()
method, where the baseline
argument takes a
distr()
distribution object with which we specify the
corresponding Normal distribution, and we specify
trt_ref = 1
to indicate that the baseline distribution
corresponds to treatment 1. (Strictly speaking,
type = "response"
is unnecessary here, since the identity
link function was used.)
<- predict(mix_fit_FE,
mix_pred_FE baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
trt_ref = 1)
mix_pred_FE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[4] -1.26 0.54 -2.28 -1.63 -1.27 -0.89 -0.21 1564 2770 1
#> pred[1] -0.72 0.22 -1.14 -0.87 -0.72 -0.58 -0.29 3556 3638 1
#> pred[2] -2.54 0.40 -3.33 -2.81 -2.53 -2.27 -1.77 5100 3406 1
#> pred[3] -1.21 0.54 -2.28 -1.58 -1.21 -0.85 -0.13 2366 3038 1
#> pred[5] -1.57 0.57 -2.69 -1.95 -1.58 -1.17 -0.47 1656 2573 1
plot(mix_pred_FE)
<- predict(mix_fit_RE,
mix_pred_RE baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
trt_ref = 1)
mix_pred_RE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[4] -1.25 0.65 -2.51 -1.67 -1.24 -0.84 0.05 1799 1935 1
#> pred[1] -0.73 0.22 -1.17 -0.88 -0.73 -0.58 -0.31 4051 3817 1
#> pred[2] -2.57 0.54 -3.62 -2.89 -2.57 -2.24 -1.53 4082 2741 1
#> pred[3] -1.24 0.66 -2.52 -1.64 -1.24 -0.82 0.05 3016 2476 1
#> pred[5] -1.53 0.79 -3.07 -2.01 -1.53 -1.06 0.04 1830 1872 1
plot(mix_pred_RE)
If the baseline
argument is omitted, predictions of mean
off-time reduction will be produced for every arm-based study
in the network based on their estimated baseline response \(\mu_j\):
<- predict(mix_fit_FE, type = "response")
mix_pred_FE_studies
mix_pred_FE_studies#> ---------------------------------------------------------------------- Study: 1 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[1: 4] -1.66 0.46 -2.57 -1.97 -1.66 -1.35 -0.73 1773 2308 1
#> pred[1: 1] -1.12 0.43 -1.96 -1.42 -1.11 -0.82 -0.31 3431 2745 1
#> pred[1: 2] -2.94 0.51 -3.94 -3.28 -2.92 -2.60 -1.95 3163 2927 1
#> pred[1: 3] -1.61 0.39 -2.33 -1.87 -1.62 -1.34 -0.85 3321 2745 1
#> pred[1: 5] -1.96 0.50 -2.94 -2.31 -1.97 -1.63 -1.00 1872 2319 1
#>
#> ---------------------------------------------------------------------- Study: 2 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[2: 4] -1.17 0.51 -2.16 -1.53 -1.18 -0.83 -0.17 1429 2350 1
#> pred[2: 1] -0.64 0.26 -1.15 -0.81 -0.64 -0.46 -0.11 5027 3218 1
#> pred[2: 2] -2.45 0.24 -2.92 -2.62 -2.45 -2.28 -1.99 4600 3784 1
#> pred[2: 3] -1.12 0.52 -2.12 -1.48 -1.13 -0.77 -0.10 2135 2645 1
#> pred[2: 5] -1.48 0.55 -2.51 -1.85 -1.49 -1.11 -0.41 1540 2193 1
#>
#> ---------------------------------------------------------------------- Study: 3 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[3: 4] -1.12 0.42 -1.95 -1.40 -1.11 -0.84 -0.28 1776 2446 1
#> pred[3: 1] -0.58 0.36 -1.29 -0.82 -0.58 -0.34 0.12 4552 3294 1
#> pred[3: 2] -2.40 0.38 -3.17 -2.65 -2.40 -2.14 -1.65 4155 3226 1
#> pred[3: 3] -1.07 0.47 -2.00 -1.38 -1.07 -0.75 -0.13 2816 2808 1
#> pred[3: 5] -1.42 0.47 -2.36 -1.73 -1.42 -1.12 -0.52 2004 2545 1
plot(mix_pred_FE_studies)
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
<- posterior_ranks(mix_fit_FE))
(mix_ranks #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[4] 3.48 0.72 2 3 3 4 5 2118 NA 1
#> rank[1] 4.65 0.76 2 5 5 5 5 1863 NA 1
#> rank[2] 1.05 0.27 1 1 1 1 2 2313 2316 1
#> rank[3] 3.55 0.91 2 3 4 4 5 3138 NA 1
#> rank[5] 2.27 0.64 1 2 2 2 4 2398 2801 1
plot(mix_ranks)
<- posterior_rank_probs(mix_fit_FE))
(mix_rankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5]
#> d[4] 0.00 0.05 0.50 0.37 0.08
#> d[1] 0.00 0.04 0.06 0.11 0.79
#> d[2] 0.96 0.03 0.01 0.00 0.00
#> d[3] 0.00 0.16 0.26 0.46 0.12
#> d[5] 0.04 0.73 0.17 0.05 0.01
plot(mix_rankprobs)
<- posterior_rank_probs(mix_fit_FE, cumulative = TRUE))
(mix_cumrankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5]
#> d[4] 0.00 0.05 0.55 0.92 1
#> d[1] 0.00 0.04 0.10 0.21 1
#> d[2] 0.96 0.99 1.00 1.00 1
#> d[3] 0.00 0.16 0.41 0.88 1
#> d[5] 0.04 0.77 0.94 0.99 1
plot(mix_cumrankprobs)