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 1641
#> d[ARB] -0.40 0.00 0.05 -0.48 -0.43 -0.40 -0.37 -0.31 2179
#> d[CCB] -0.20 0.00 0.03 -0.26 -0.22 -0.20 -0.18 -0.14 2246
#> d[Diuretic] 0.06 0.00 0.06 -0.06 0.02 0.06 0.10 0.17 1834
#> d[Placebo] -0.19 0.00 0.05 -0.29 -0.22 -0.19 -0.16 -0.09 1526
#> lp__ -37970.33 0.09 3.71 -37978.57 -37972.63 -37969.98 -37967.64 -37964.11 1883
#> 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 May 23 11:35:23 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(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.50 -0.38 -0.33 -0.28 -0.17 2159
#> d[ARB] -0.40 0.00 0.10 -0.61 -0.46 -0.40 -0.34 -0.22 2202
#> d[CCB] -0.17 0.00 0.06 -0.30 -0.21 -0.17 -0.13 -0.04 2071
#> d[Diuretic] 0.07 0.00 0.09 -0.11 0.02 0.07 0.13 0.24 2105
#> d[Placebo] -0.22 0.00 0.09 -0.41 -0.27 -0.21 -0.16 -0.05 1595
#> lp__ -37980.58 0.21 6.84 -37995.04 -37984.88 -37980.23 -37975.76 -37968.11 1035
#> tau 0.13 0.00 0.04 0.06 0.10 0.13 0.16 0.24 973
#> 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 May 23 11:35:57 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(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.3 (on 48 data points)
#> pD: 27.1
#> DIC: 105.3
<- dic(db_fit_RE))
(dic_RE #> Residual deviance: 53.4 (on 48 data points)
#> pD: 38.3
#> DIC: 91.7
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.10 -0.06 -0.02 0.06 1848 2424 1
#> d[ACE Inhibitor] -0.36 0.05 -0.47 -0.40 -0.36 -0.32 -0.25 4545 3730 1
#> d[ARB] -0.45 0.06 -0.58 -0.50 -0.45 -0.41 -0.33 3685 3094 1
#> d[CCB] -0.25 0.05 -0.36 -0.29 -0.25 -0.22 -0.15 2860 3091 1
#> d[Placebo] -0.25 0.06 -0.36 -0.28 -0.25 -0.21 -0.14 4309 3367 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.24 -0.13 -0.07 -0.02 0.11 2151 2445 1
#> d[ACE Inhibitor] -0.40 0.09 -0.58 -0.46 -0.40 -0.34 -0.23 3870 2782 1
#> d[ARB] -0.47 0.11 -0.70 -0.54 -0.47 -0.40 -0.26 3847 2884 1
#> d[CCB] -0.24 0.08 -0.41 -0.29 -0.24 -0.18 -0.08 4088 2657 1
#> d[Placebo] -0.29 0.09 -0.48 -0.34 -0.29 -0.23 -0.12 3605 2815 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.08 0.24 3991 3703 1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.00 0.02 0.03 0.06 0.18 3971 3777 1
#> pred[New 1: ARB] 0.04 0.04 0.00 0.01 0.03 0.05 0.17 3995 3739 1
#> pred[New 1: CCB] 0.05 0.05 0.00 0.02 0.03 0.06 0.20 3998 3776 1
#> pred[New 1: Diuretic] 0.06 0.07 0.01 0.02 0.04 0.08 0.25 3971 3703 1
#> pred[New 1: Placebo] 0.05 0.05 0.01 0.02 0.03 0.06 0.20 3977 3814 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.07 0.01 0.02 0.04 0.08 0.25 3904 3932 1
#> pred[New 1: ACE Inhibitor] 0.05 0.05 0.00 0.02 0.03 0.06 0.18 3907 3790 1
#> pred[New 1: ARB] 0.04 0.05 0.00 0.01 0.03 0.05 0.18 3911 3873 1
#> pred[New 1: CCB] 0.05 0.06 0.01 0.02 0.03 0.06 0.22 3912 3914 1
#> pred[New 1: Diuretic] 0.07 0.07 0.01 0.02 0.04 0.08 0.26 3929 3830 1
#> pred[New 1: Placebo] 0.05 0.06 0.01 0.02 0.03 0.06 0.21 3915 3792 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 6113 3302 1
#> pred[AASK: ACE Inhibitor] 0.12 0.01 0.10 0.12 0.12 0.13 0.15 4371 3018 1
#> pred[AASK: ARB] 0.12 0.01 0.09 0.11 0.12 0.13 0.15 4151 3047 1
#> pred[AASK: CCB] 0.14 0.01 0.12 0.13 0.14 0.15 0.18 5305 3073 1
#> pred[AASK: Diuretic] 0.18 0.02 0.14 0.17 0.18 0.19 0.22 3844 2608 1
#> pred[AASK: Placebo] 0.14 0.02 0.11 0.13 0.14 0.15 0.17 3432 3069 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.05 3079 2347 1
#> pred[ALLHAT: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03 0.04 4528 2696 1
#> pred[ALLHAT: ARB] 0.03 0.00 0.02 0.03 0.03 0.03 0.04 3789 2644 1
#> pred[ALLHAT: CCB] 0.04 0.00 0.03 0.03 0.04 0.04 0.05 4169 2782 1
#> pred[ALLHAT: Diuretic] 0.05 0.01 0.04 0.04 0.05 0.05 0.06 4586 2697 1
#> pred[ALLHAT: Placebo] 0.03 0.00 0.03 0.03 0.03 0.04 0.04 4279 2783 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.03 0.03 0.05 5867 3087 1
#> pred[ALPINE: ACE Inhibitor] 0.02 0.01 0.01 0.01 0.02 0.02 0.03 6436 3000 1
#> pred[ALPINE: ARB] 0.02 0.01 0.01 0.01 0.02 0.02 0.03 6262 2998 1
#> pred[ALPINE: CCB] 0.02 0.01 0.01 0.02 0.02 0.03 0.04 6415 3109 1
#> pred[ALPINE: Diuretic] 0.03 0.01 0.01 0.02 0.03 0.03 0.05 6335 3048 1
#> pred[ALPINE: Placebo] 0.02 0.01 0.01 0.02 0.02 0.03 0.04 6311 3190 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 2754 2018 1
#> pred[ANBP-2: ACE Inhibitor] 0.05 0.01 0.04 0.04 0.05 0.05 0.06 4032 2236 1
#> pred[ANBP-2: ARB] 0.05 0.01 0.03 0.04 0.05 0.05 0.06 3825 2586 1
#> pred[ANBP-2: CCB] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 3686 2217 1
#> pred[ANBP-2: Diuretic] 0.07 0.01 0.05 0.07 0.07 0.08 0.10 4172 2395 1
#> pred[ANBP-2: Placebo] 0.05 0.01 0.04 0.05 0.05 0.06 0.07 4062 2488 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 5897 2964 1
#> pred[ASCOT: ACE Inhibitor] 0.08 0.01 0.07 0.08 0.08 0.09 0.10 2612 2602 1
#> pred[ASCOT: ARB] 0.08 0.01 0.06 0.07 0.08 0.08 0.09 2549 2660 1
#> pred[ASCOT: CCB] 0.10 0.01 0.08 0.09 0.10 0.10 0.11 2393 2514 1
#> pred[ASCOT: Diuretic] 0.12 0.01 0.10 0.11 0.12 0.13 0.14 2392 2476 1
#> pred[ASCOT: Placebo] 0.09 0.01 0.08 0.09 0.09 0.10 0.11 1909 2065 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 5623 2793 1
#> pred[CAPPP: ACE Inhibitor] 0.05 0.00 0.05 0.05 0.05 0.06 0.06 2519 2399 1
#> pred[CAPPP: ARB] 0.05 0.01 0.04 0.05 0.05 0.05 0.06 2612 2319 1
#> pred[CAPPP: CCB] 0.06 0.00 0.05 0.06 0.06 0.07 0.07 2956 2975 1
#> pred[CAPPP: Diuretic] 0.08 0.01 0.07 0.08 0.08 0.08 0.10 2650 2376 1
#> pred[CAPPP: Placebo] 0.06 0.01 0.05 0.06 0.06 0.06 0.07 1959 2036 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 2917 2500 1
#> pred[CHARM: ACE Inhibitor] 0.07 0.01 0.05 0.06 0.07 0.07 0.09 4579 2806 1
#> pred[CHARM: ARB] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 4843 2573 1
#> pred[CHARM: CCB] 0.08 0.01 0.06 0.07 0.08 0.08 0.10 4191 2907 1
#> pred[CHARM: Diuretic] 0.10 0.02 0.07 0.09 0.10 0.11 0.13 4258 2656 1
#> pred[CHARM: Placebo] 0.07 0.01 0.06 0.07 0.07 0.08 0.09 5066 2653 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.17 0.21 0.23 0.25 0.30 2510 1934 1
#> pred[DREAM: ACE Inhibitor] 0.17 0.02 0.13 0.16 0.17 0.18 0.22 3595 1943 1
#> pred[DREAM: ARB] 0.16 0.02 0.12 0.14 0.16 0.17 0.21 3733 2400 1
#> pred[DREAM: CCB] 0.20 0.03 0.15 0.18 0.20 0.21 0.26 3538 2316 1
#> pred[DREAM: Diuretic] 0.24 0.03 0.19 0.22 0.24 0.26 0.32 3700 2533 1
#> pred[DREAM: Placebo] 0.19 0.02 0.15 0.17 0.19 0.20 0.24 3850 2212 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 3634 2857 1
#> pred[EWPH: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.04 0.05 0.07 5026 2897 1
#> pred[EWPH: ARB] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 4730 2810 1
#> pred[EWPH: CCB] 0.05 0.01 0.04 0.05 0.05 0.06 0.08 4889 2818 1
#> pred[EWPH: Diuretic] 0.07 0.01 0.05 0.06 0.07 0.07 0.09 5531 2995 1
#> pred[EWPH: Placebo] 0.05 0.01 0.03 0.04 0.05 0.06 0.07 5620 2892 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.05 0.06 3066 2216 1
#> pred[FEVER: ACE Inhibitor] 0.03 0.00 0.02 0.03 0.03 0.03 0.04 4512 2410 1
#> pred[FEVER: ARB] 0.03 0.00 0.02 0.03 0.03 0.03 0.04 4270 2400 1
#> pred[FEVER: CCB] 0.04 0.01 0.03 0.03 0.03 0.04 0.05 4424 2271 1
#> pred[FEVER: Diuretic] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 4500 2649 1
#> pred[FEVER: Placebo] 0.03 0.00 0.03 0.03 0.03 0.04 0.04 4406 2696 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 6405 3123 1
#> pred[HAPPHY: ACE Inhibitor] 0.02 0 0.01 0.02 0.02 0.02 0.02 4832 2670 1
#> pred[HAPPHY: ARB] 0.02 0 0.01 0.02 0.02 0.02 0.02 4342 3265 1
#> pred[HAPPHY: CCB] 0.02 0 0.02 0.02 0.02 0.02 0.03 5247 2991 1
#> pred[HAPPHY: Diuretic] 0.03 0 0.02 0.02 0.03 0.03 0.03 3651 2485 1
#> pred[HAPPHY: Placebo] 0.02 0 0.02 0.02 0.02 0.02 0.02 3754 2882 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 3050 2368 1
#> pred[HOPE: ACE Inhibitor] 0.04 0.01 0.03 0.04 0.04 0.05 0.05 4373 2509 1
#> pred[HOPE: ARB] 0.04 0.01 0.03 0.04 0.04 0.04 0.05 4522 2309 1
#> pred[HOPE: CCB] 0.05 0.01 0.04 0.04 0.05 0.05 0.07 4191 2234 1
#> pred[HOPE: Diuretic] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 4315 2580 1
#> pred[HOPE: Placebo] 0.05 0.01 0.03 0.04 0.05 0.05 0.06 4625 2448 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.09 3356 2492 1
#> pred[INSIGHT: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05 0.06 4483 2395 1
#> pred[INSIGHT: ARB] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 4337 2426 1
#> pred[INSIGHT: CCB] 0.06 0.01 0.04 0.05 0.05 0.06 0.07 4688 2562 1
#> pred[INSIGHT: Diuretic] 0.07 0.01 0.05 0.06 0.07 0.08 0.09 4647 2362 1
#> pred[INSIGHT: Placebo] 0.05 0.01 0.04 0.05 0.05 0.06 0.07 4378 2355 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 7583 2680 1
#> pred[INVEST: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06 0.07 2601 2766 1
#> pred[INVEST: ARB] 0.06 0.01 0.05 0.05 0.06 0.06 0.07 2402 2499 1
#> pred[INVEST: CCB] 0.07 0.00 0.06 0.07 0.07 0.07 0.08 2562 2716 1
#> pred[INVEST: Diuretic] 0.09 0.01 0.07 0.08 0.09 0.09 0.10 2482 2441 1
#> pred[INVEST: Placebo] 0.07 0.01 0.06 0.06 0.07 0.07 0.08 1933 1987 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 7801 3242 1
#> pred[LIFE: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.06 0.07 2903 2613 1
#> pred[LIFE: ARB] 0.06 0.01 0.04 0.05 0.06 0.06 0.07 2546 2831 1
#> pred[LIFE: CCB] 0.07 0.01 0.06 0.07 0.07 0.07 0.08 3133 2323 1
#> pred[LIFE: Diuretic] 0.09 0.01 0.07 0.08 0.09 0.09 0.11 2854 2524 1
#> pred[LIFE: Placebo] 0.07 0.01 0.05 0.06 0.07 0.07 0.08 2132 2160 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 3941 2968 1
#> pred[MRC-E: ACE Inhibitor] 0.02 0 0.02 0.02 0.02 0.02 0.03 5669 3190 1
#> pred[MRC-E: ARB] 0.02 0 0.01 0.02 0.02 0.02 0.03 5191 3061 1
#> pred[MRC-E: CCB] 0.03 0 0.02 0.02 0.02 0.03 0.03 5142 2939 1
#> pred[MRC-E: Diuretic] 0.03 0 0.02 0.03 0.03 0.03 0.04 4241 2917 1
#> pred[MRC-E: Placebo] 0.02 0 0.02 0.02 0.02 0.03 0.03 5258 3355 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 7564 3108 1
#> pred[NORDIL: ACE Inhibitor] 0.04 0.00 0.03 0.03 0.04 0.04 0.04 3202 2973 1
#> pred[NORDIL: ARB] 0.03 0.00 0.03 0.03 0.03 0.04 0.04 2827 2657 1
#> pred[NORDIL: CCB] 0.04 0.00 0.04 0.04 0.04 0.04 0.05 3294 3229 1
#> pred[NORDIL: Diuretic] 0.05 0.01 0.04 0.05 0.05 0.06 0.06 2918 2796 1
#> pred[NORDIL: Placebo] 0.04 0.00 0.03 0.04 0.04 0.04 0.05 2412 2293 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 2879 2096 1
#> pred[PEACE: ACE Inhibitor] 0.10 0.01 0.08 0.09 0.10 0.11 0.13 4361 2270 1
#> pred[PEACE: ARB] 0.09 0.01 0.07 0.09 0.09 0.10 0.13 4242 2696 1
#> pred[PEACE: CCB] 0.12 0.02 0.09 0.11 0.12 0.13 0.15 4018 2452 1
#> pred[PEACE: Diuretic] 0.15 0.02 0.11 0.13 0.15 0.16 0.19 4068 2226 1
#> pred[PEACE: Placebo] 0.11 0.01 0.09 0.10 0.11 0.12 0.14 4432 2528 1
#>
#> ------------------------------------------------------------------ Study: SCOPE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[SCOPE: Beta Blocker] 0.07 0.01 0.05 0.06 0.06 0.07 0.09 3094 2672 1
#> pred[SCOPE: ACE Inhibitor] 0.05 0.01 0.03 0.04 0.05 0.05 0.06 4688 2704 1
#> pred[SCOPE: ARB] 0.04 0.01 0.03 0.04 0.04 0.05 0.06 5458 3159 1
#> pred[SCOPE: CCB] 0.06 0.01 0.04 0.05 0.05 0.06 0.07 4418 2414 1
#> pred[SCOPE: Diuretic] 0.07 0.01 0.05 0.06 0.07 0.08 0.10 4049 2585 1
#> pred[SCOPE: Placebo] 0.05 0.01 0.04 0.05 0.05 0.06 0.07 4947 2972 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 2719 2199 1
#> pred[SHEP: ACE Inhibitor] 0.06 0.01 0.05 0.06 0.06 0.07 0.08 3945 2443 1
#> pred[SHEP: ARB] 0.06 0.01 0.04 0.05 0.06 0.06 0.08 4039 2522 1
#> pred[SHEP: CCB] 0.07 0.01 0.05 0.07 0.07 0.08 0.10 3705 2674 1
#> pred[SHEP: Diuretic] 0.09 0.01 0.07 0.08 0.09 0.10 0.12 4530 2687 1
#> pred[SHEP: Placebo] 0.07 0.01 0.05 0.06 0.07 0.08 0.09 4668 2858 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 4582 2742 1
#> pred[STOP-2: ACE Inhibitor] 0.04 0.00 0.03 0.04 0.04 0.04 0.05 3084 2615 1
#> pred[STOP-2: ARB] 0.04 0.00 0.03 0.03 0.04 0.04 0.05 2982 2499 1
#> pred[STOP-2: CCB] 0.05 0.00 0.04 0.04 0.05 0.05 0.05 3992 2687 1
#> pred[STOP-2: Diuretic] 0.06 0.01 0.05 0.05 0.06 0.06 0.07 3148 2326 1
#> pred[STOP-2: Placebo] 0.04 0.00 0.03 0.04 0.04 0.05 0.05 2375 2326 1
#>
#> ------------------------------------------------------------------ Study: VALUE ----
#>
#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[VALUE: Beta Blocker] 0.20 0.02 0.15 0.18 0.20 0.21 0.25 2916 2288 1
#> pred[VALUE: ACE Inhibitor] 0.15 0.02 0.11 0.13 0.14 0.16 0.19 4341 2696 1
#> pred[VALUE: ARB] 0.14 0.02 0.10 0.13 0.14 0.15 0.17 4244 2623 1
#> pred[VALUE: CCB] 0.17 0.02 0.13 0.16 0.17 0.18 0.21 4060 2306 1
#> pred[VALUE: Diuretic] 0.21 0.03 0.16 0.19 0.21 0.23 0.27 3999 2690 1
#> pred[VALUE: Placebo] 0.16 0.02 0.12 0.15 0.16 0.17 0.21 4373 2649 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.19 0.43 5 5 5 5 6 2514 NA 1
#> rank[ACE Inhibitor] 1.86 0.56 1 2 2 2 3 3920 3432 1
#> rank[ARB] 1.27 0.52 1 1 1 1 3 3350 2907 1
#> rank[CCB] 3.70 0.52 3 3 4 4 4 3990 2423 1
#> rank[Diuretic] 5.79 0.42 5 6 6 6 6 2427 NA 1
#> rank[Placebo] 3.19 0.61 2 3 3 4 4 3402 2603 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.01 0.79 0.2
#> d[ACE Inhibitor] 0.23 0.69 0.07 0.01 0.00 0.0
#> d[ARB] 0.76 0.21 0.02 0.00 0.00 0.0
#> d[CCB] 0.00 0.02 0.27 0.70 0.01 0.0
#> d[Diuretic] 0.00 0.00 0.00 0.00 0.20 0.8
#> d[Placebo] 0.01 0.08 0.63 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.01 0.8 1
#> d[ACE Inhibitor] 0.23 0.92 0.99 1.00 1.0 1
#> d[ARB] 0.76 0.97 1.00 1.00 1.0 1
#> d[CCB] 0.00 0.02 0.29 0.99 1.0 1
#> d[Diuretic] 0.00 0.00 0.00 0.00 0.2 1
#> d[Placebo] 0.01 0.09 0.72 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.↩︎