library(multinma)
options(mc.cores = parallel::detectCores())
This vignette describes the analysis of data on the number of new
cases of diabetes in 22 trials of 6 antihypertensive drugs (Elliott and Meyer 2007; Dias et al. 2011). The data are
available in this package as diabetes
:
head(diabetes)
#> studyn studyc trtn trtc r n time
#> 1 1 MRC-E 1 Diuretic 43 1081 5.8
#> 2 1 MRC-E 2 Placebo 34 2213 5.8
#> 3 1 MRC-E 3 Beta Blocker 37 1102 5.8
#> 4 2 EWPH 1 Diuretic 29 416 4.7
#> 5 2 EWPH 2 Placebo 20 424 4.7
#> 6 3 SHEP 1 Diuretic 140 1631 3.0
We begin by setting up the network. We have arm-level count data
giving the number of new cases of diabetes (r
) out of the
total (n
) in each arm, so we use the function
set_agd_arm()
. For computational efficiency, we let “Beta
Blocker” be set as the network reference treatment by default. Elliott and Meyer (2007) and Dias et
al. (2011) use
“Diuretic” as the reference, but it is a simple matter to transform the
results after fitting the NMA model.1
<- set_agd_arm(diabetes,
db_net study = studyc,
trt = trtc,
r = r,
n = n)
db_net#> A network with 22 AgD studies (arm-based).
#>
#> ------------------------------------------------------- AgD studies (arm-based) ----
#> Study Treatment arms
#> AASK 3: Beta Blocker | ACE Inhibitor | CCB
#> ALLHAT 3: ACE Inhibitor | CCB | Diuretic
#> ALPINE 2: ARB | Diuretic
#> ANBP-2 2: ACE Inhibitor | Diuretic
#> ASCOT 2: Beta Blocker | CCB
#> CAPPP 2: Beta Blocker | ACE Inhibitor
#> CHARM 2: ARB | Placebo
#> DREAM 2: ACE Inhibitor | Placebo
#> EWPH 2: Diuretic | Placebo
#> FEVER 2: CCB | Placebo
#> ... plus 12 more studies
#>
#> Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 6
#> Total number of studies: 22
#> Reference treatment is: Beta Blocker
#> Network is connected
We also have details of length of follow-up in years in each trial
(time
), which we will use as an offset with a cloglog link
function to model the data as rates. We do not have to specify this in
the function set_agd_arm()
: any additional columns in the
data (e.g. offsets or covariates, here the column time
)
will automatically be made available in the network.
Plot the network structure.
plot(db_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. We specify
that a cloglog link will be used with link = "cloglog"
(the
Binomial likelihood is the default for these data), and specify the log
follow-up time offset using the regression formula
regression = ~offset(log(time))
.
<- nma(db_net,
db_fit_FE trt_effects = "fixed",
link = "cloglog",
regression = ~offset(log(time)),
prior_intercept = normal(scale = 100),
prior_trt = 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.
#> Note: Setting "Beta Blocker" as the network reference treatment.
Basic parameter summaries are given by the print()
method:
db_fit_FE#> A fixed effects NMA with a binomial likelihood (cloglog link).
#> Regression model: ~offset(log(time)).
#> 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
#> d[ACE Inhibitor] -0.30 0.00 0.05 -0.39 -0.33 -0.30 -0.27 -0.21 1227
#> d[ARB] -0.39 0.00 0.05 -0.49 -0.43 -0.39 -0.36 -0.30 2010
#> d[CCB] -0.20 0.00 0.03 -0.26 -0.22 -0.20 -0.17 -0.13 1671
#> d[Diuretic] 0.06 0.00 0.06 -0.05 0.02 0.06 0.09 0.17 1523
#> d[Placebo] -0.19 0.00 0.05 -0.29 -0.22 -0.19 -0.16 -0.09 1179
#> lp__ -38119.60 0.09 3.70 -38127.60 -38121.85 -38119.25 -38117.00 -38113.27 1748
#> Rhat
#> d[ACE Inhibitor] 1
#> d[ARB] 1
#> d[CCB] 1
#> d[Diuretic] 1
#> d[Placebo] 1
#> lp__ 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:56:28 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(db_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using
the plot_prior_posterior()
function:
plot_prior_posterior(db_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(db_net,
db_fit_RE trt_effects = "random",
link = "cloglog",
regression = ~offset(log(time)),
prior_intercept = normal(scale = 10),
prior_trt = normal(scale = 10),
prior_het = half_normal(scale = 5),
init_r = 0.5)
#> 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.
#> Note: Setting "Beta Blocker" as the network reference treatment.
Basic parameter summaries are given by the print()
method:
db_fit_RE#> A random effects NMA with a binomial likelihood (cloglog link).
#> Regression model: ~offset(log(time)).
#> 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
#> d[ACE Inhibitor] -0.33 0.00 0.08 -0.49 -0.38 -0.33 -0.28 -0.18 1438
#> d[ARB] -0.40 0.00 0.10 -0.59 -0.46 -0.40 -0.34 -0.21 1606
#> d[CCB] -0.17 0.00 0.07 -0.30 -0.21 -0.17 -0.13 -0.03 1802
#> d[Diuretic] 0.07 0.00 0.09 -0.10 0.01 0.07 0.13 0.25 1737
#> d[Placebo] -0.21 0.00 0.09 -0.40 -0.27 -0.21 -0.16 -0.04 1237
#> lp__ -38070.21 0.22 6.87 -38084.72 -38074.58 -38069.88 -38065.62 -38057.55 967
#> tau 0.13 0.00 0.05 0.06 0.10 0.13 0.16 0.24 845
#> Rhat
#> d[ACE Inhibitor] 1.00
#> d[ARB] 1.00
#> d[CCB] 1.00
#> d[Diuretic] 1.00
#> d[Placebo] 1.00
#> lp__ 1.00
#> tau 1.01
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 17:56:51 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(db_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(db_fit_RE, prior = c("trt", "het"))
Model fit can be checked using the dic()
function:
<- dic(db_fit_FE))
(dic_FE #> Residual deviance: 78.5 (on 48 data points)
#> pD: 27.3
#> DIC: 105.8
<- dic(db_fit_RE))
(dic_RE #> Residual deviance: 53.4 (on 48 data points)
#> pD: 38.1
#> DIC: 91.5
The FE model is a very poor fit to the data, with a residual deviance much higher than the number of data points. The RE model fits the data better, and has a much lower DIC; we prefer the RE model.
We can also examine the residual deviance contributions with the
corresponding plot()
method.
plot(dic_FE)
plot(dic_RE)
For comparison with Elliott and Meyer (2007) and Dias et al. (2011), we can produce relative effects
against “Diuretic” using the relative_effects()
function
with trt_ref = "Diuretic"
:
<- relative_effects(db_fit_FE, trt_ref = "Diuretic"))
(db_releff_FE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Beta Blocker] -0.06 0.06 -0.17 -0.09 -0.06 -0.02 0.05 1534 2566 1
#> d[ACE Inhibitor] -0.36 0.05 -0.47 -0.39 -0.36 -0.32 -0.25 4582 3450 1
#> d[ARB] -0.45 0.06 -0.58 -0.49 -0.45 -0.41 -0.33 3181 3052 1
#> d[CCB] -0.25 0.05 -0.36 -0.29 -0.25 -0.22 -0.15 2595 3234 1
#> d[Placebo] -0.25 0.06 -0.36 -0.28 -0.25 -0.21 -0.14 3992 2883 1
plot(db_releff_FE, ref_line = 0)
<- relative_effects(db_fit_RE, trt_ref = "Diuretic"))
(db_releff_RE #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Beta Blocker] -0.07 0.09 -0.25 -0.13 -0.07 -0.01 0.10 1776 2284 1
#> d[ACE Inhibitor] -0.40 0.09 -0.58 -0.46 -0.40 -0.34 -0.24 4083 3200 1
#> d[ARB] -0.47 0.11 -0.69 -0.54 -0.47 -0.40 -0.26 3293 2854 1
#> d[CCB] -0.24 0.09 -0.41 -0.30 -0.24 -0.18 -0.07 3577 3316 1
#> d[Placebo] -0.29 0.09 -0.47 -0.34 -0.28 -0.23 -0.12 3835 2964 1
plot(db_releff_RE, ref_line = 0)
Dias et al. (2011) produce absolute predictions of
the probability of developing diabetes after three years, assuming a
Normal distribution on the baseline cloglog probability of developing
diabetes on diuretic treatment with mean \(-4.2\) and precision \(1.11\). We can replicate these results
using the predict()
method. We specify a data frame of
newdata
, containing the time
offset(s) at
which to produce predictions (here only 3 years). The
baseline
argument takes a distr()
distribution
object with which we specify the corresponding Normal distribution on
the baseline cloglog probability, and we set
trt_ref = "Diuretic"
to indicate that the baseline
distribution corresponds to “Diuretic” rather than the network reference
“Beta Blocker”. We set type = "response"
to produce
predicted event probabilities (type = "link"
would produce
predicted cloglog probabilities).
<- predict(db_fit_FE,
db_pred_FE newdata = data.frame(time = 3),
baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
trt_ref = "Diuretic",
type = "response")
db_pred_FE#> ------------------------------------------------------------------ Study: New 1 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[New 1: Beta Blocker] 0.06 0.06 0.01 0.02 0.04 0.07 0.23 3871 3891 1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.00 0.02 0.03 0.06 0.18 3895 3852 1
#> pred[New 1: ARB] 0.04 0.04 0.00 0.01 0.03 0.05 0.16 3899 3889 1
#> pred[New 1: CCB] 0.05 0.05 0.01 0.02 0.03 0.06 0.20 3887 3852 1
#> pred[New 1: Diuretic] 0.06 0.06 0.01 0.02 0.04 0.08 0.24 3905 3891 1
#> pred[New 1: Placebo] 0.05 0.05 0.01 0.02 0.03 0.06 0.20 3904 3813 1
plot(db_pred_FE)
<- predict(db_fit_RE,
db_pred_RE newdata = data.frame(time = 3),
baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
trt_ref = "Diuretic",
type = "response")
db_pred_RE#> ------------------------------------------------------------------ Study: New 1 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[New 1: Beta Blocker] 0.06 0.06 0.01 0.02 0.04 0.08 0.25 3784 3932 1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.00 0.02 0.03 0.06 0.18 3789 3853 1
#> pred[New 1: ARB] 0.04 0.05 0.00 0.01 0.03 0.05 0.17 3784 3892 1
#> pred[New 1: CCB] 0.05 0.06 0.01 0.02 0.03 0.06 0.21 3779 3930 1
#> pred[New 1: Diuretic] 0.07 0.07 0.01 0.02 0.04 0.08 0.26 3776 3852 1
#> pred[New 1: Placebo] 0.05 0.05 0.01 0.02 0.03 0.06 0.21 3772 3852 1
plot(db_pred_RE)
If the baseline
and newdata
arguments are
omitted, predicted probabilities will be produced for every study in the
network based on their follow-up times and estimated baseline cloglog
probabilities \(\mu_j\):
<- predict(db_fit_RE, type = "response")
db_pred_RE_studies
db_pred_RE_studies#> ------------------------------------------------------------------- Study: AASK ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[AASK: Beta Blocker] 0.17 0.02 0.14 0.16 0.17 0.18 0.20 4723 3246 1
#> pred[AASK: ACE Inhibitor] 0.13 0.01 0.10 0.12 0.12 0.13 0.15 3772 3061 1
#> pred[AASK: ARB] 0.12 0.01 0.09 0.11 0.12 0.13 0.15 3487 2967 1
#> pred[AASK: CCB] 0.15 0.02 0.12 0.14 0.14 0.16 0.18 4204 3078 1
#> pred[AASK: Diuretic] 0.18 0.02 0.14 0.17 0.18 0.19 0.22 3476 3153 1
#> pred[AASK: Placebo] 0.14 0.02 0.11 0.13 0.14 0.15 0.17 2844 3267 1
#>
#> ----------------------------------------------------------------- Study: ALLHAT ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ALLHAT: Beta Blocker] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 2494 2348 1
#> pred[ALLHAT: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03 0.04 3888 2440 1
#> pred[ALLHAT: ARB] 0.03 0.00 0.02 0.03 0.03 0.03 0.04 3598 2730 1
#> pred[ALLHAT: CCB] 0.04 0.00 0.03 0.03 0.04 0.04 0.05 3558 2108 1
#> pred[ALLHAT: Diuretic] 0.05 0.01 0.04 0.04 0.05 0.05 0.06 3860 2665 1
#> pred[ALLHAT: Placebo] 0.03 0.00 0.03 0.03 0.03 0.04 0.05 3687 2481 1
#>
#> ----------------------------------------------------------------- Study: ALPINE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ALPINE: Beta Blocker] 0.03 0.01 0.01 0.02 0.02 0.03 0.05 5764 3098 1
#> pred[ALPINE: ACE Inhibitor] 0.02 0.01 0.01 0.01 0.02 0.02 0.03 6749 3056 1
#> pred[ALPINE: ARB] 0.02 0.01 0.01 0.01 0.02 0.02 0.03 7128 2936 1
#> pred[ALPINE: CCB] 0.02 0.01 0.01 0.02 0.02 0.03 0.04 6932 3156 1
#> pred[ALPINE: Diuretic] 0.03 0.01 0.01 0.02 0.03 0.03 0.05 7259 2950 1
#> pred[ALPINE: Placebo] 0.02 0.01 0.01 0.02 0.02 0.03 0.04 6941 2811 1
#>
#> ----------------------------------------------------------------- Study: ANBP-2 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ANBP-2: Beta Blocker] 0.07 0.01 0.05 0.06 0.07 0.07 0.09 2612 2178 1
#> pred[ANBP-2: ACE Inhibitor] 0.05 0.01 0.04 0.04 0.05 0.05 0.06 4351 2710 1
#> pred[ANBP-2: ARB] 0.05 0.01 0.03 0.04 0.05 0.05 0.06 3664 2424 1
#> pred[ANBP-2: CCB] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 3639 2552 1
#> pred[ANBP-2: Diuretic] 0.07 0.01 0.06 0.07 0.07 0.08 0.09 4509 2419 1
#> pred[ANBP-2: Placebo] 0.05 0.01 0.04 0.05 0.05 0.06 0.07 4125 2600 1
#>
#> ------------------------------------------------------------------ Study: ASCOT ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[ASCOT: Beta Blocker] 0.11 0.00 0.10 0.11 0.11 0.11 0.12 5041 2398 1
#> pred[ASCOT: ACE Inhibitor] 0.08 0.01 0.07 0.08 0.08 0.09 0.10 2056 2404 1
#> pred[ASCOT: ARB] 0.08 0.01 0.06 0.07 0.08 0.08 0.09 1963 2425 1
#> pred[ASCOT: CCB] 0.10 0.01 0.08 0.09 0.10 0.10 0.11 2089 2170 1
#> pred[ASCOT: Diuretic] 0.12 0.01 0.10 0.11 0.12 0.13 0.14 2014 2428 1
#> pred[ASCOT: Placebo] 0.09 0.01 0.08 0.09 0.09 0.10 0.11 1512 2075 1
#>
#> ------------------------------------------------------------------ Study: CAPPP ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[CAPPP: Beta Blocker] 0.07 0.00 0.07 0.07 0.07 0.08 0.08 4890 3103 1
#> pred[CAPPP: ACE Inhibitor] 0.05 0.00 0.05 0.05 0.05 0.06 0.06 1680 2220 1
#> pred[CAPPP: ARB] 0.05 0.01 0.04 0.05 0.05 0.05 0.06 2092 2140 1
#> pred[CAPPP: CCB] 0.06 0.00 0.05 0.06 0.06 0.07 0.07 2585 2575 1
#> pred[CAPPP: Diuretic] 0.08 0.01 0.07 0.08 0.08 0.09 0.10 2135 2501 1
#> pred[CAPPP: Placebo] 0.06 0.01 0.05 0.06 0.06 0.06 0.07 1449 1855 1
#>
#> ------------------------------------------------------------------ Study: CHARM ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[CHARM: Beta Blocker] 0.09 0.01 0.07 0.08 0.09 0.10 0.12 2361 2250 1
#> pred[CHARM: ACE Inhibitor] 0.07 0.01 0.05 0.06 0.07 0.07 0.09 4092 2726 1
#> pred[CHARM: ARB] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 4581 3002 1
#> pred[CHARM: CCB] 0.08 0.01 0.06 0.07 0.08 0.08 0.10 3558 2717 1
#> pred[CHARM: Diuretic] 0.10 0.01 0.07 0.09 0.10 0.10 0.13 4112 2808 1
#> pred[CHARM: Placebo] 0.07 0.01 0.06 0.07 0.07 0.08 0.10 4625 2875 1
#>
#> ------------------------------------------------------------------ Study: DREAM ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[DREAM: Beta Blocker] 0.23 0.03 0.18 0.21 0.23 0.24 0.29 2492 2621 1
#> pred[DREAM: ACE Inhibitor] 0.17 0.02 0.13 0.16 0.17 0.18 0.21 4212 3075 1
#> pred[DREAM: ARB] 0.16 0.02 0.12 0.14 0.16 0.17 0.21 3712 2893 1
#> pred[DREAM: CCB] 0.20 0.03 0.15 0.18 0.19 0.21 0.25 3651 2834 1
#> pred[DREAM: Diuretic] 0.24 0.03 0.19 0.22 0.24 0.26 0.31 4067 2569 1
#> pred[DREAM: Placebo] 0.19 0.02 0.15 0.17 0.19 0.20 0.24 4509 2798 1
#>
#> ------------------------------------------------------------------- Study: EWPH ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[EWPH: Beta Blocker] 0.06 0.01 0.04 0.05 0.06 0.07 0.09 3506 2967 1
#> pred[EWPH: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.04 0.05 0.07 4771 3198 1
#> pred[EWPH: ARB] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 4052 3280 1
#> pred[EWPH: CCB] 0.05 0.01 0.04 0.05 0.05 0.06 0.08 4108 3243 1
#> pred[EWPH: Diuretic] 0.07 0.01 0.05 0.06 0.07 0.07 0.10 4852 3138 1
#> pred[EWPH: Placebo] 0.05 0.01 0.03 0.04 0.05 0.06 0.07 4749 2947 1
#>
#> ------------------------------------------------------------------ Study: FEVER ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[FEVER: Beta Blocker] 0.04 0.01 0.03 0.04 0.04 0.04 0.05 2565 2642 1
#> pred[FEVER: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03 0.04 3706 2837 1
#> pred[FEVER: ARB] 0.03 0.00 0.02 0.02 0.03 0.03 0.04 3737 2459 1
#> pred[FEVER: CCB] 0.03 0.00 0.03 0.03 0.03 0.04 0.05 3697 2988 1
#> pred[FEVER: Diuretic] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 3942 2928 1
#> pred[FEVER: Placebo] 0.03 0.00 0.02 0.03 0.03 0.04 0.04 3795 2544 1
#>
#> ----------------------------------------------------------------- Study: HAPPHY ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[HAPPHY: Beta Blocker] 0.02 0 0.02 0.02 0.02 0.03 0.03 4877 3398 1
#> pred[HAPPHY: ACE Inhibitor] 0.02 0 0.01 0.02 0.02 0.02 0.02 3655 2613 1
#> pred[HAPPHY: ARB] 0.02 0 0.01 0.02 0.02 0.02 0.02 3736 3039 1
#> pred[HAPPHY: CCB] 0.02 0 0.02 0.02 0.02 0.02 0.03 3830 2783 1
#> pred[HAPPHY: Diuretic] 0.03 0 0.02 0.02 0.03 0.03 0.03 3133 2826 1
#> pred[HAPPHY: Placebo] 0.02 0 0.02 0.02 0.02 0.02 0.02 2870 3062 1
#>
#> ------------------------------------------------------------------- Study: HOPE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[HOPE: Beta Blocker] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 2246 2122 1
#> pred[HOPE: ACE Inhibitor] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 3893 3079 1
#> pred[HOPE: ARB] 0.04 0.01 0.03 0.04 0.04 0.04 0.05 3666 2986 1
#> pred[HOPE: CCB] 0.05 0.01 0.04 0.04 0.05 0.05 0.07 3231 2836 1
#> pred[HOPE: Diuretic] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 4129 2995 1
#> pred[HOPE: Placebo] 0.05 0.01 0.04 0.04 0.05 0.05 0.06 4106 2752 1
#>
#> ---------------------------------------------------------------- Study: INSIGHT ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[INSIGHT: Beta Blocker] 0.07 0.01 0.05 0.06 0.06 0.07 0.08 2646 2628 1
#> pred[INSIGHT: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05 0.06 3590 2713 1
#> pred[INSIGHT: ARB] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 3423 2317 1
#> pred[INSIGHT: CCB] 0.06 0.01 0.04 0.05 0.05 0.06 0.07 3740 3088 1
#> pred[INSIGHT: Diuretic] 0.07 0.01 0.05 0.06 0.07 0.07 0.09 4217 2346 1
#> pred[INSIGHT: Placebo] 0.05 0.01 0.04 0.05 0.05 0.06 0.07 3560 2604 1
#>
#> ----------------------------------------------------------------- Study: INVEST ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[INVEST: Beta Blocker] 0.08 0.00 0.08 0.08 0.08 0.08 0.09 6040 2876 1
#> pred[INVEST: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06 0.07 2019 2434 1
#> pred[INVEST: ARB] 0.06 0.01 0.05 0.05 0.06 0.06 0.07 1978 2354 1
#> pred[INVEST: CCB] 0.07 0.00 0.06 0.07 0.07 0.07 0.08 2229 2585 1
#> pred[INVEST: Diuretic] 0.09 0.01 0.07 0.08 0.09 0.09 0.11 2069 2556 1
#> pred[INVEST: Placebo] 0.07 0.01 0.06 0.06 0.07 0.07 0.08 1500 1955 1
#>
#> ------------------------------------------------------------------- Study: LIFE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[LIFE: Beta Blocker] 0.08 0.00 0.07 0.08 0.08 0.08 0.09 7078 2746 1
#> pred[LIFE: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06 0.07 2268 2786 1
#> pred[LIFE: ARB] 0.06 0.01 0.05 0.05 0.06 0.06 0.07 1975 2442 1
#> pred[LIFE: CCB] 0.07 0.01 0.06 0.07 0.07 0.07 0.08 2941 2434 1
#> pred[LIFE: Diuretic] 0.09 0.01 0.07 0.08 0.09 0.09 0.11 2351 2509 1
#> pred[LIFE: Placebo] 0.07 0.01 0.05 0.06 0.07 0.07 0.08 1728 2165 1
#>
#> ------------------------------------------------------------------ Study: MRC-E ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[MRC-E: Beta Blocker] 0.03 0 0.02 0.03 0.03 0.03 0.04 3514 3056 1
#> pred[MRC-E: ACE Inhibitor] 0.02 0 0.02 0.02 0.02 0.02 0.03 4767 3028 1
#> pred[MRC-E: ARB] 0.02 0 0.01 0.02 0.02 0.02 0.03 4670 3162 1
#> pred[MRC-E: CCB] 0.03 0 0.02 0.02 0.02 0.03 0.03 4412 2808 1
#> pred[MRC-E: Diuretic] 0.03 0 0.02 0.03 0.03 0.03 0.04 4325 2982 1
#> pred[MRC-E: Placebo] 0.02 0 0.02 0.02 0.02 0.03 0.03 4216 3534 1
#>
#> ----------------------------------------------------------------- Study: NORDIL ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[NORDIL: Beta Blocker] 0.05 0.00 0.04 0.05 0.05 0.05 0.06 6403 3275 1
#> pred[NORDIL: ACE Inhibitor] 0.04 0.00 0.03 0.03 0.04 0.04 0.04 2326 2387 1
#> pred[NORDIL: ARB] 0.03 0.00 0.03 0.03 0.03 0.04 0.04 2037 2491 1
#> pred[NORDIL: CCB] 0.04 0.00 0.04 0.04 0.04 0.04 0.05 2673 2639 1
#> pred[NORDIL: Diuretic] 0.05 0.01 0.04 0.05 0.05 0.06 0.06 2339 2475 1
#> pred[NORDIL: Placebo] 0.04 0.00 0.03 0.04 0.04 0.04 0.05 1833 2091 1
#>
#> ------------------------------------------------------------------ Study: PEACE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[PEACE: Beta Blocker] 0.14 0.02 0.10 0.13 0.14 0.15 0.18 2153 2171 1
#> pred[PEACE: ACE Inhibitor] 0.10 0.01 0.08 0.09 0.10 0.11 0.13 3960 2836 1
#> pred[PEACE: ARB] 0.09 0.01 0.07 0.09 0.09 0.10 0.13 3461 2694 1
#> pred[PEACE: CCB] 0.12 0.02 0.09 0.11 0.12 0.13 0.15 3122 2631 1
#> pred[PEACE: Diuretic] 0.15 0.02 0.11 0.13 0.15 0.16 0.19 3703 2711 1
#> pred[PEACE: Placebo] 0.11 0.01 0.09 0.10 0.11 0.12 0.14 4168 2706 1
#>
#> ------------------------------------------------------------------ Study: SCOPE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[SCOPE: Beta Blocker] 0.06 0.01 0.05 0.06 0.06 0.07 0.09 2151 2008 1
#> pred[SCOPE: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05 0.06 3607 2638 1
#> pred[SCOPE: ARB] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 3985 2702 1
#> pred[SCOPE: CCB] 0.06 0.01 0.04 0.05 0.05 0.06 0.07 2959 2261 1
#> pred[SCOPE: Diuretic] 0.07 0.01 0.05 0.06 0.07 0.08 0.09 3938 2644 1
#> pred[SCOPE: Placebo] 0.05 0.01 0.04 0.05 0.05 0.06 0.07 4557 2312 1
#>
#> ------------------------------------------------------------------- Study: SHEP ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[SHEP: Beta Blocker] 0.09 0.01 0.06 0.08 0.09 0.09 0.12 2162 2221 1
#> pred[SHEP: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 4209 2968 1
#> pred[SHEP: ARB] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 3542 2638 1
#> pred[SHEP: CCB] 0.07 0.01 0.05 0.07 0.07 0.08 0.10 3452 2651 1
#> pred[SHEP: Diuretic] 0.09 0.01 0.07 0.08 0.09 0.10 0.12 4340 2571 1
#> pred[SHEP: Placebo] 0.07 0.01 0.05 0.06 0.07 0.08 0.09 4156 2845 1
#>
#> ----------------------------------------------------------------- Study: STOP-2 ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[STOP-2: Beta Blocker] 0.05 0.00 0.04 0.05 0.05 0.06 0.06 3564 2894 1
#> pred[STOP-2: ACE Inhibitor] 0.04 0.00 0.03 0.04 0.04 0.04 0.05 2323 2639 1
#> pred[STOP-2: ARB] 0.04 0.00 0.03 0.03 0.04 0.04 0.04 2316 2854 1
#> pred[STOP-2: CCB] 0.05 0.00 0.04 0.04 0.05 0.05 0.05 3358 3095 1
#> pred[STOP-2: Diuretic] 0.06 0.01 0.05 0.05 0.06 0.06 0.07 2722 2768 1
#> pred[STOP-2: Placebo] 0.04 0.00 0.03 0.04 0.04 0.05 0.05 1929 1935 1
#>
#> ------------------------------------------------------------------ Study: VALUE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[VALUE: Beta Blocker] 0.20 0.03 0.15 0.18 0.19 0.21 0.25 2623 2133 1
#> pred[VALUE: ACE Inhibitor] 0.15 0.02 0.11 0.13 0.14 0.16 0.19 3370 2479 1
#> pred[VALUE: ARB] 0.14 0.02 0.10 0.13 0.14 0.15 0.17 3700 2326 1
#> pred[VALUE: CCB] 0.17 0.02 0.13 0.16 0.17 0.18 0.21 3514 2400 1
#> pred[VALUE: Diuretic] 0.21 0.03 0.16 0.19 0.21 0.22 0.27 3360 2698 1
#> pred[VALUE: Placebo] 0.16 0.02 0.12 0.15 0.16 0.17 0.21 3284 2180 1
plot(db_pred_RE_studies)
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
<- posterior_ranks(db_fit_RE))
(db_ranks #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[Beta Blocker] 5.18 0.43 5 5 5 5 6 1971 NA 1
#> rank[ACE Inhibitor] 1.85 0.54 1 2 2 2 3 3345 2913 1
#> rank[ARB] 1.27 0.52 1 1 1 1 3 3185 3079 1
#> rank[CCB] 3.69 0.54 3 3 4 4 4 2977 2623 1
#> rank[Diuretic] 5.80 0.41 5 6 6 6 6 2280 NA 1
#> rank[Placebo] 3.21 0.60 2 3 3 4 4 2826 2278 1
plot(db_ranks)
<- posterior_rank_probs(db_fit_RE))
(db_rankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
#> d[Beta Blocker] 0.00 0.00 0.00 0.02 0.78 0.2
#> d[ACE Inhibitor] 0.23 0.71 0.06 0.01 0.00 0.0
#> d[ARB] 0.77 0.20 0.03 0.00 0.00 0.0
#> d[CCB] 0.00 0.02 0.27 0.69 0.01 0.0
#> d[Diuretic] 0.00 0.00 0.00 0.00 0.20 0.8
#> d[Placebo] 0.01 0.07 0.64 0.28 0.01 0.0
plot(db_rankprobs)
<- posterior_rank_probs(db_fit_RE, cumulative = TRUE))
(db_cumrankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
#> d[Beta Blocker] 0.00 0.00 0.00 0.02 0.8 1
#> d[ACE Inhibitor] 0.23 0.93 0.99 1.00 1.0 1
#> d[ARB] 0.77 0.97 1.00 1.00 1.0 1
#> d[CCB] 0.00 0.02 0.30 0.99 1.0 1
#> d[Diuretic] 0.00 0.00 0.00 0.00 0.2 1
#> d[Placebo] 0.01 0.08 0.71 0.99 1.0 1
plot(db_cumrankprobs)
The gain in efficiency here from using “Beta Blocker” as the network reference treatment instead of “Diuretic” is considerable - around 4-8 times, in terms of effective samples per second. The functions in this package will always attempt to choose a default network reference treatment that maximises computational efficiency and stability. If you have chosen an alternative network reference treatment and the model runs very slowly or has low effective sample size, this is a likely cause.↩︎