library(multinma)
options(mc.cores = parallel::detectCores())
This vignette describes the analysis of treatments for
moderate-to-severe plaque psoriasis from an HTA report (Woolacott et al. 2006), replicating the
analysis in NICE Technical Support Document 2 (Dias et al. 2011). The
data are available in this package as hta_psoriasis
:
head(hta_psoriasis)
#> studyn studyc year trtn trtc sample_size PASI50 PASI75 PASI90
#> 1 1 Elewski 2004 1 Supportive care 193 12 5 1
#> 2 1 Elewski 2004 2 Etanercept 25 mg 196 59 46 21
#> 3 1 Elewski 2004 3 Etanercept 50 mg 194 54 56 40
#> 4 2 Gottlieb 2003 1 Supportive care 55 5 1 0
#> 5 2 Gottlieb 2003 2 Etanercept 25 mg 57 23 11 6
#> 6 3 Lebwohl 2003 1 Supportive care 122 13 5 1
Outcomes are ordered multinomial success/failure to achieve 50%, 75%, or 90% reduction in symptoms on the Psoriasis Area and Severity Index (PASI) scale. Some studies report ordered outcomes at all three cutpoints, others only one or two:
::filter(hta_psoriasis, studyc %in% c("Elewski", "Gordon", "ACD2058g", "Altmeyer"))
dplyr#> studyn studyc year trtn trtc sample_size PASI50 PASI75 PASI90
#> 1 1 Elewski 2004 1 Supportive care 193 12 5 1
#> 2 1 Elewski 2004 2 Etanercept 25 mg 196 59 46 21
#> 3 1 Elewski 2004 3 Etanercept 50 mg 194 54 56 40
#> 4 5 Gordon 2003 1 Supportive care 187 18 8 NA
#> 5 5 Gordon 2003 4 Efalizumab 369 118 98 NA
#> 6 6 ACD2058g 2004 1 Supportive care 170 25 NA NA
#> 7 6 ACD2058g 2004 4 Efalizumab 162 99 NA NA
#> 8 10 Altmeyer 1994 1 Supportive care 51 NA 1 NA
#> 9 10 Altmeyer 1994 6 Fumaderm 49 NA 12 NA
Here, the outcome counts are given as “exclusive” counts. That is, for a study reporting all outcomes (e.g. Elewski), the counts represent the categories 50 < PASI < 75, 75 < PASI < 90, and 90 < PASI < 100, and the corresponding columns are named by the lower end of the interval. Missing values are used where studies only report a subset of the outcomes. For a study reporting only two outcomes, say PASI50 and PASI75 as in Gordon, the counts represent the categories 50 < PASI < 75 and 75 < PASI < 100. For a study reporting only one outcome, say PASI70 as in Altmeyer, the count represents 70 < PASI < 100. We also need the count for the lowest category (i.e. no higher outcomes achieved), which is equal to the sample size minus the counts in the other observed categories.
We begin by setting up the network. We have arm-level ordered
multinomial count data, so we use the function
set_agd_arm()
. The function multi()
helps us
to specify the ordered outcomes correctly.
<- set_agd_arm(hta_psoriasis,
pso_net study = paste(studyc, year),
trt = trtc,
r = multi(r0 = sample_size - rowSums(cbind(PASI50, PASI75, PASI90), na.rm = TRUE),
PASI50, PASI75, PASI90,inclusive = FALSE,
type = "ordered"))
pso_net#> A network with 16 AgD studies (arm-based).
#>
#> ------------------------------------------------------- AgD studies (arm-based) ----
#> Study Treatment arms
#> ACD2058g 2004 2: Supportive care | Efalizumab
#> ACD2600g 2004 2: Supportive care | Efalizumab
#> Altmeyer 1994 2: Supportive care | Fumaderm
#> Chaudari 2001 2: Supportive care | Infliximab
#> Elewski 2004 3: Supportive care | Etanercept 25 mg | Etanercept 50 mg
#> Ellis 1991 3: Supportive care | Ciclosporin | Ciclosporin
#> Gordon 2003 2: Supportive care | Efalizumab
#> Gottlieb 2003 2: Supportive care | Etanercept 25 mg
#> Gottlieb 2004 3: Supportive care | Infliximab | Infliximab
#> Guenther 1991 2: Supportive care | Ciclosporin
#> ... plus 6 more studies
#>
#> Outcome type: ordered (4 categories)
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 8
#> Total number of studies: 16
#> Reference treatment is: Supportive care
#> Network is connected
Plot the network structure.
plot(pso_net, weight_edges = TRUE, weight_nodes = TRUE) +
# Nudge the legend over
::theme(legend.box.spacing = ggplot2::unit(0.75, "in"),
ggplot2plot.margin = ggplot2::margin(0.1, 0, 0.1, 0.75, "in"))
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"
, using a probit link
function link = "probit"
. We use \(\mathrm{N}(0, 10^2)\) prior distributions
for the treatment effects \(d_k\), and
\(\mathrm{N}(0, 100^2)\) prior
distributions for the 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 = 10))
#> A Normal prior distribution: location = 0, scale = 10.
#> 50% of the prior density lies between -6.74 and 6.74.
#> 95% of the prior density lies between -19.6 and 19.6.
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.
We also need to specify prior distributions for the latent cutpoints
\(c_\textrm{PASI75}\) and \(c_\textrm{PASI90}\) on the underlying scale
- here the PASI standardised mean difference due to the probit link (the
cutpoint \(c_\textrm{PASI50}=0\)). To
make these easier to reason about, we actually specify priors on the
differences between adjacent cutpoints, e.g. \(c_\textrm{PASI90} - c_\textrm{PASI75}\) and
\(c_\textrm{PASI75} -
c_\textrm{PASI50}\). These can be given any positive-valued prior
distribution, and Stan will automatically impose the necessary ordering
constraints behind the scenes. We choose to give these implicit flat
priors flat()
.
The model is fitted using the nma()
function.
<- nma(pso_net,
pso_fit_FE trt_effects = "fixed",
link = "probit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_aux = flat())
#> Note: Setting "Supportive care" as the network reference treatment.
Basic parameter summaries are given by the print()
method:
pso_fit_FE#> A fixed effects NMA with a ordered likelihood (probit link).
#> Inference for Stan model: ordered_multinomial.
#> 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[Ciclosporin] 1.93 0.01 0.34 1.31 1.68 1.91 2.15 2.63 1320 1
#> d[Efalizumab] 1.19 0.00 0.06 1.08 1.15 1.19 1.23 1.30 1911 1
#> d[Etanercept 25 mg] 1.51 0.00 0.09 1.33 1.44 1.51 1.57 1.69 2287 1
#> d[Etanercept 50 mg] 1.92 0.00 0.10 1.73 1.85 1.92 1.98 2.11 2514 1
#> d[Fumaderm] 1.48 0.01 0.49 0.61 1.14 1.45 1.79 2.51 2348 1
#> d[Infliximab] 2.34 0.01 0.28 1.83 2.16 2.33 2.52 2.91 2842 1
#> d[Methotrexate] 1.62 0.01 0.44 0.77 1.32 1.62 1.92 2.51 1540 1
#> lp__ -3516.17 0.09 3.58 -3524.40 -3518.36 -3515.86 -3513.56 -3510.21 1454 1
#> cc[PASI50] 0.00 NaN 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN
#> cc[PASI75] 0.76 0.00 0.03 0.70 0.73 0.75 0.78 0.82 5385 1
#> cc[PASI90] 1.56 0.00 0.05 1.46 1.53 1.56 1.60 1.67 6851 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 18:07:36 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).
Note: the treatment effects are the opposite sign to those in TSD 2 (Dias et al. 2011). This is because we parameterise the linear predictor as \(\mu_j + d_k + c_m\), rather than \(\mu_j + d_k - c_m\). The interpretation here thus follows that of a standard binomial probit (or logit) regression; SMDs (or log ORs) greater than zero mean that the treatment increases the probability of an event compared to the comparator (and less than zero mean a reduction in probability). Here higher outcomes are positive, and all of the active treatments are estimated to increase the response (i.e. a greater reduction) on the PASI scale compared to the network reference (supportive care).
By default, summaries of the study-specific intercepts \(\mu_j\) are hidden, but could be examined
by changing the pars
argument:
# Not run
print(pso_fit_FE, pars = c("d", "mu", "cc"))
The prior and posterior distributions can be compared visually using
the plot_prior_posterior()
function:
plot_prior_posterior(pso_fit_FE)
Focusing specifically on the cutpoints we see that these are highly identified by the data, which is why the implicit flat priors work for these parameters.
plot_prior_posterior(pso_fit_FE, prior = "aux")
We now fit a random effects model using the nma()
function with trt_effects = "random"
. Again, we use \(\mathrm{N}(0, 10^2)\) prior distributions
for the treatment effects \(d_k\),
\(\mathrm{N}(0, 100^2)\) prior
distributions for the study-specific intercepts \(\mu_j\), implicit flat prior distributions
for the latent cutpoints, and we additionally use a \(\textrm{half-N}(2.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 = 10))
#> A Normal prior distribution: location = 0, scale = 10.
#> 50% of the prior density lies between -6.74 and 6.74.
#> 95% of the prior density lies between -19.6 and 19.6.
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 = 2.5))
#> A half-Normal prior distribution: location = 0, scale = 2.5.
#> 50% of the prior density lies between 0 and 1.69.
#> 95% of the prior density lies between 0 and 4.9.
Fitting the RE model
<- nma(pso_net,
pso_fit_RE trt_effects = "random",
link = "probit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_aux = flat(),
prior_het = half_normal(scale = 2.5),
adapt_delta = 0.99)
#> Note: Setting "Supportive care" as the network reference treatment.
Basic parameter summaries are given by the print()
method:
pso_fit_RE#> A random effects NMA with a ordered likelihood (probit link).
#> Inference for Stan model: ordered_multinomial.
#> 4 chains, each with iter=5000; warmup=2500; thin=1;
#> post-warmup draws per chain=2500, total post-warmup draws=10000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> d[Ciclosporin] 2.02 0.01 0.42 1.29 1.74 1.99 2.27 2.95 2437 1
#> d[Efalizumab] 1.19 0.00 0.18 0.80 1.10 1.19 1.27 1.56 3907 1
#> d[Etanercept 25 mg] 1.53 0.00 0.24 1.04 1.40 1.52 1.65 2.03 5361 1
#> d[Etanercept 50 mg] 1.92 0.00 0.28 1.35 1.78 1.92 2.06 2.49 4927 1
#> d[Fumaderm] 1.48 0.01 0.61 0.34 1.09 1.46 1.85 2.78 8575 1
#> d[Infliximab] 2.32 0.00 0.38 1.55 2.08 2.31 2.55 3.08 6738 1
#> d[Methotrexate] 1.71 0.01 0.63 0.57 1.31 1.68 2.06 3.09 3771 1
#> lp__ -3523.36 0.20 6.74 -3537.03 -3527.97 -3523.18 -3518.57 -3510.96 1170 1
#> tau 0.31 0.01 0.22 0.02 0.15 0.27 0.43 0.84 860 1
#> cc[PASI50] 0.00 NaN 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN
#> cc[PASI75] 0.76 0.00 0.03 0.70 0.73 0.76 0.78 0.82 18439 1
#> cc[PASI90] 1.56 0.00 0.05 1.46 1.53 1.56 1.60 1.67 21740 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jan 9 18:08:43 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(pso_fit_RE, pars = c("d", "cc", "mu", "delta"))
The prior and posterior distributions can be compared visually using
the plot_prior_posterior()
function:
plot_prior_posterior(pso_fit_RE, prior = c("trt", "aux", "het"))
Model fit can be checked using the dic()
function:
<- dic(pso_fit_FE))
(dic_FE #> Residual deviance: 75.1 (on 58 data points)
#> pD: 25.6
#> DIC: 100.6
<- dic(pso_fit_RE))
(dic_RE #> Residual deviance: 62.7 (on 58 data points)
#> pD: 33.4
#> DIC: 96.1
The random effects model has a lower DIC and the residual deviance is closer to the number of data points, so is preferred in this case.
We can also examine the residual deviance contributions with the
corresponding plot()
method.
plot(dic_FE)
plot(dic_RE)
Most data points are fit well, with posterior mean residual deviances close to the degrees of freedom. The Meffert 1997 study has a substantially higher residual deviance contribution, which could be investigated further to see why this study appears to be an outlier.
Dias et al. (2011) produce absolute predictions of
probability of achieving responses at each PASI cutoff, assuming a
Normal distribution for the baseline probit probability of PASI50
response on supportive care with mean \(-1.097\) and precision \(123\). We can replicate these results using
the predict()
method. The baseline
argument
takes a distr()
distribution object, with which we specify
the corresponding Normal distribution. We set
type = "response"
to produce predicted probabilities
(type = "link"
would produce predicted probit
probabilities).
<- predict(pso_fit_FE,
pred_FE baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5),
type = "response")
pred_FE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[Supportive care, PASI50] 0.14 0.02 0.10 0.12 0.14 0.15 0.18 3631 3432 1
#> pred[Supportive care, PASI75] 0.03 0.01 0.02 0.03 0.03 0.04 0.05 3798 3795 1
#> pred[Supportive care, PASI90] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 4176 3581 1
#> pred[Ciclosporin, PASI50] 0.78 0.10 0.57 0.72 0.79 0.85 0.94 1465 1547 1
#> pred[Ciclosporin, PASI75] 0.53 0.13 0.28 0.43 0.52 0.62 0.79 1462 1773 1
#> pred[Ciclosporin, PASI90] 0.24 0.11 0.08 0.16 0.23 0.31 0.50 1481 1647 1
#> pred[Efalizumab, PASI50] 0.54 0.04 0.45 0.51 0.53 0.57 0.62 3061 3193 1
#> pred[Efalizumab, PASI75] 0.25 0.04 0.19 0.23 0.25 0.28 0.33 3153 3293 1
#> pred[Efalizumab, PASI90] 0.07 0.02 0.04 0.06 0.07 0.08 0.11 3409 3327 1
#> pred[Etanercept 25 mg, PASI50] 0.66 0.05 0.56 0.63 0.66 0.69 0.75 2707 3322 1
#> pred[Etanercept 25 mg, PASI75] 0.37 0.05 0.27 0.33 0.37 0.40 0.46 2738 3788 1
#> pred[Etanercept 25 mg, PASI90] 0.13 0.03 0.08 0.11 0.12 0.14 0.19 2906 3155 1
#> pred[Etanercept 50 mg, PASI50] 0.79 0.04 0.71 0.77 0.79 0.82 0.86 2779 3495 1
#> pred[Etanercept 50 mg, PASI75] 0.53 0.05 0.42 0.49 0.53 0.56 0.63 2769 3678 1
#> pred[Etanercept 50 mg, PASI90] 0.23 0.04 0.16 0.20 0.23 0.26 0.32 2892 3481 1
#> pred[Fumaderm, PASI50] 0.63 0.17 0.31 0.52 0.64 0.76 0.92 2793 2222 1
#> pred[Fumaderm, PASI75] 0.37 0.17 0.10 0.24 0.35 0.48 0.76 2781 2275 1
#> pred[Fumaderm, PASI90] 0.15 0.12 0.02 0.06 0.11 0.20 0.45 2763 2260 1
#> pred[Infliximab, PASI50] 0.88 0.06 0.75 0.85 0.89 0.93 0.97 2927 2764 1
#> pred[Infliximab, PASI75] 0.68 0.10 0.47 0.62 0.68 0.75 0.86 2890 2912 1
#> pred[Infliximab, PASI90] 0.38 0.11 0.19 0.30 0.37 0.45 0.61 2904 2750 1
#> pred[Methotrexate, PASI50] 0.68 0.15 0.37 0.59 0.70 0.80 0.93 1654 1948 1
#> pred[Methotrexate, PASI75] 0.42 0.16 0.13 0.30 0.41 0.53 0.75 1642 1918 1
#> pred[Methotrexate, PASI90] 0.17 0.11 0.03 0.09 0.15 0.23 0.45 1658 1987 1
plot(pred_FE)
<- predict(pso_fit_RE,
pred_RE baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5),
type = "response")
pred_RE#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[Supportive care, PASI50] 0.14 0.02 0.10 0.12 0.14 0.15 0.18 9886 9993 1
#> pred[Supportive care, PASI75] 0.03 0.01 0.02 0.03 0.03 0.04 0.05 10471 9755 1
#> pred[Supportive care, PASI90] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 11189 9998 1
#> pred[Ciclosporin, PASI50] 0.80 0.10 0.57 0.74 0.81 0.88 0.97 2891 2681 1
#> pred[Ciclosporin, PASI75] 0.56 0.15 0.28 0.45 0.55 0.66 0.87 2857 2603 1
#> pred[Ciclosporin, PASI90] 0.28 0.14 0.08 0.18 0.25 0.35 0.62 2896 2773 1
#> pred[Efalizumab, PASI50] 0.53 0.08 0.37 0.49 0.54 0.58 0.69 5018 3190 1
#> pred[Efalizumab, PASI75] 0.26 0.06 0.14 0.22 0.25 0.29 0.40 5080 3220 1
#> pred[Efalizumab, PASI90] 0.07 0.03 0.03 0.06 0.07 0.09 0.14 5195 3201 1
#> pred[Etanercept 25 mg, PASI50] 0.66 0.09 0.46 0.61 0.67 0.72 0.83 5965 4297 1
#> pred[Etanercept 25 mg, PASI75] 0.38 0.09 0.20 0.32 0.37 0.43 0.58 5967 4145 1
#> pred[Etanercept 25 mg, PASI90] 0.14 0.06 0.05 0.10 0.13 0.16 0.27 6030 4166 1
#> pred[Etanercept 50 mg, PASI50] 0.79 0.08 0.59 0.75 0.80 0.84 0.92 5795 3901 1
#> pred[Etanercept 50 mg, PASI75] 0.53 0.11 0.30 0.47 0.53 0.59 0.74 5828 4022 1
#> pred[Etanercept 50 mg, PASI90] 0.24 0.09 0.09 0.19 0.23 0.28 0.44 5865 4122 1
#> pred[Fumaderm, PASI50] 0.63 0.19 0.22 0.50 0.64 0.78 0.95 8804 6212 1
#> pred[Fumaderm, PASI75] 0.37 0.20 0.06 0.22 0.35 0.51 0.82 8801 6264 1
#> pred[Fumaderm, PASI90] 0.16 0.14 0.01 0.06 0.11 0.21 0.55 8856 6381 1
#> pred[Infliximab, PASI50] 0.87 0.08 0.67 0.83 0.89 0.93 0.98 7199 5104 1
#> pred[Infliximab, PASI75] 0.67 0.13 0.38 0.59 0.68 0.76 0.89 7201 5220 1
#> pred[Infliximab, PASI90] 0.37 0.14 0.13 0.28 0.37 0.46 0.67 7243 5377 1
#> pred[Methotrexate, PASI50] 0.70 0.18 0.29 0.58 0.72 0.83 0.98 4456 2702 1
#> pred[Methotrexate, PASI75] 0.45 0.21 0.10 0.29 0.43 0.59 0.89 4447 2722 1
#> pred[Methotrexate, PASI90] 0.21 0.17 0.02 0.09 0.16 0.28 0.67 4419 2758 1
plot(pred_RE)
If instead of information on the baseline PASI 50 response probit
probability we have PASI 50 event counts, we can use these to construct
a Beta distribution for the baseline probability of PASI 50 response.
For example, if 56 out of 408 individuals achieved PASI 50 response on
supportive care in the target population of interest, the appropriate
Beta distribution for the response probability would be \(\textrm{Beta}(56, 408-56)\). We can specify
this Beta distribution for the baseline response using the
baseline_type = "reponse"
argument (the default is
"link"
, used above for the baseline probit
probability).
<- predict(pso_fit_FE,
pred_FE_beta baseline = distr(qbeta, 56, 408-56),
baseline_type = "response",
type = "response")
pred_FE_beta#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[Supportive care, PASI50] 0.14 0.02 0.10 0.12 0.14 0.15 0.17 3940 4002 1
#> pred[Supportive care, PASI75] 0.03 0.01 0.02 0.03 0.03 0.04 0.05 4150 3816 1
#> pred[Supportive care, PASI90] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 4410 3897 1
#> pred[Ciclosporin, PASI50] 0.78 0.10 0.58 0.72 0.79 0.85 0.94 1410 1921 1
#> pred[Ciclosporin, PASI75] 0.53 0.13 0.29 0.43 0.52 0.62 0.79 1412 1768 1
#> pred[Ciclosporin, PASI90] 0.24 0.11 0.08 0.16 0.23 0.30 0.50 1427 1631 1
#> pred[Efalizumab, PASI50] 0.54 0.04 0.46 0.51 0.54 0.56 0.61 3267 3490 1
#> pred[Efalizumab, PASI75] 0.25 0.03 0.19 0.23 0.25 0.27 0.32 3330 3666 1
#> pred[Efalizumab, PASI90] 0.07 0.01 0.04 0.06 0.07 0.08 0.10 3498 3813 1
#> pred[Etanercept 25 mg, PASI50] 0.66 0.04 0.57 0.63 0.66 0.69 0.74 2871 2833 1
#> pred[Etanercept 25 mg, PASI75] 0.37 0.05 0.28 0.33 0.37 0.40 0.46 2982 3380 1
#> pred[Etanercept 25 mg, PASI90] 0.13 0.03 0.08 0.11 0.12 0.14 0.18 3225 3574 1
#> pred[Etanercept 50 mg, PASI50] 0.79 0.04 0.72 0.77 0.79 0.82 0.86 2867 3036 1
#> pred[Etanercept 50 mg, PASI75] 0.53 0.05 0.43 0.49 0.52 0.56 0.62 3050 3340 1
#> pred[Etanercept 50 mg, PASI90] 0.23 0.04 0.16 0.20 0.23 0.25 0.31 3249 3249 1
#> pred[Fumaderm, PASI50] 0.63 0.16 0.32 0.52 0.64 0.76 0.92 2796 2015 1
#> pred[Fumaderm, PASI75] 0.37 0.17 0.11 0.24 0.34 0.48 0.75 2834 2037 1
#> pred[Fumaderm, PASI90] 0.15 0.12 0.02 0.06 0.11 0.20 0.45 2801 1995 1
#> pred[Infliximab, PASI50] 0.88 0.05 0.76 0.85 0.89 0.92 0.97 2960 2277 1
#> pred[Infliximab, PASI75] 0.68 0.10 0.48 0.62 0.69 0.75 0.87 2921 2330 1
#> pred[Infliximab, PASI90] 0.38 0.11 0.19 0.30 0.37 0.45 0.61 2932 2409 1
#> pred[Methotrexate, PASI50] 0.68 0.15 0.37 0.59 0.70 0.79 0.92 1582 1997 1
#> pred[Methotrexate, PASI75] 0.42 0.16 0.14 0.29 0.41 0.53 0.74 1577 2074 1
#> pred[Methotrexate, PASI90] 0.17 0.11 0.03 0.09 0.15 0.23 0.44 1590 2137 1
plot(pred_FE_beta)
<- predict(pso_fit_RE,
pred_RE_beta baseline = distr(qbeta, 56, 408-56),
baseline_type = "response",
type = "response")
pred_RE_beta#> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[Supportive care, PASI50] 0.14 0.02 0.11 0.13 0.14 0.15 0.17 9843 9606 1
#> pred[Supportive care, PASI75] 0.03 0.01 0.02 0.03 0.03 0.04 0.05 10410 9613 1
#> pred[Supportive care, PASI90] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 12274 9519 1
#> pred[Ciclosporin, PASI50] 0.80 0.10 0.57 0.74 0.81 0.88 0.97 3035 2838 1
#> pred[Ciclosporin, PASI75] 0.56 0.15 0.28 0.45 0.55 0.66 0.87 3003 2786 1
#> pred[Ciclosporin, PASI90] 0.28 0.14 0.08 0.18 0.25 0.35 0.62 3060 2825 1
#> pred[Efalizumab, PASI50] 0.53 0.07 0.38 0.49 0.54 0.58 0.68 4688 3388 1
#> pred[Efalizumab, PASI75] 0.26 0.06 0.14 0.22 0.25 0.29 0.39 4765 3360 1
#> pred[Efalizumab, PASI90] 0.07 0.03 0.03 0.06 0.07 0.09 0.14 4943 3392 1
#> pred[Etanercept 25 mg, PASI50] 0.66 0.09 0.47 0.61 0.67 0.71 0.83 5847 4057 1
#> pred[Etanercept 25 mg, PASI75] 0.38 0.09 0.20 0.32 0.37 0.43 0.58 5850 4035 1
#> pred[Etanercept 25 mg, PASI90] 0.14 0.06 0.05 0.10 0.13 0.16 0.27 5991 4071 1
#> pred[Etanercept 50 mg, PASI50] 0.79 0.08 0.59 0.75 0.80 0.84 0.92 5774 3766 1
#> pred[Etanercept 50 mg, PASI75] 0.53 0.11 0.30 0.47 0.53 0.59 0.74 5801 3891 1
#> pred[Etanercept 50 mg, PASI90] 0.24 0.09 0.09 0.19 0.23 0.28 0.44 5909 4035 1
#> pred[Fumaderm, PASI50] 0.63 0.19 0.22 0.50 0.64 0.78 0.95 8902 6259 1
#> pred[Fumaderm, PASI75] 0.37 0.20 0.06 0.22 0.35 0.50 0.82 8891 6286 1
#> pred[Fumaderm, PASI90] 0.16 0.14 0.01 0.06 0.11 0.21 0.55 8960 6369 1
#> pred[Infliximab, PASI50] 0.87 0.08 0.67 0.84 0.89 0.93 0.98 7029 5043 1
#> pred[Infliximab, PASI75] 0.67 0.13 0.38 0.59 0.68 0.76 0.89 7031 5042 1
#> pred[Infliximab, PASI90] 0.37 0.14 0.13 0.28 0.36 0.46 0.67 7095 5071 1
#> pred[Methotrexate, PASI50] 0.70 0.18 0.30 0.58 0.72 0.84 0.98 4471 2855 1
#> pred[Methotrexate, PASI75] 0.45 0.21 0.10 0.29 0.43 0.59 0.89 4472 2878 1
#> pred[Methotrexate, PASI90] 0.21 0.17 0.02 0.09 0.16 0.28 0.67 4484 2720 1
plot(pred_RE_beta)
(Notice that these results are equivalent to those calculated above using the Normal distribution for the baseline probit probability, since these event counts correspond to the same probit probability.)
We can modify the plots using standard ggplot2
functions. For example, to plot the cutpoints together with a colour
coding (instead of split into facets):
library(ggplot2)
plot(pred_RE, position = position_dodge(width = 0.75)) +
facet_null() +
aes(colour = Category) +
scale_colour_brewer(palette = "Blues")
If the baseline
argument is omitted, predicted
probabilities will be produced for every study in the network based on
their estimated baseline probit probability \(\mu_j\).
Treatment rankings, rank probabilities, and cumulative rank
probabilities can also be produced. We set
lower_better = FALSE
since higher outcome categories are
better (the outcomes are positive).
<- posterior_ranks(pso_fit_RE, lower_better = FALSE))
(pso_ranks #> mean sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[Supportive care] 7.99 0.11 8 8 8 8 8 6794 NA 1
#> rank[Ciclosporin] 2.77 1.27 1 2 3 4 5 6291 6805 1
#> rank[Efalizumab] 6.35 0.81 4 6 7 7 7 5929 NA 1
#> rank[Etanercept 25 mg] 4.92 1.07 3 4 5 6 7 6353 5633 1
#> rank[Etanercept 50 mg] 3.05 1.22 1 2 3 4 6 4814 5452 1
#> rank[Fumaderm] 4.91 1.93 1 3 5 7 7 8323 6866 1
#> rank[Infliximab] 1.77 1.17 1 1 1 2 5 3843 4744 1
#> rank[Methotrexate] 4.24 1.86 1 3 4 6 7 5385 5667 1
plot(pso_ranks)
<- posterior_rank_probs(pso_fit_RE, lower_better = FALSE))
(pso_rankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6] p_rank[7] p_rank[8]
#> d[Supportive care] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.99
#> d[Ciclosporin] 0.17 0.29 0.27 0.17 0.07 0.02 0.00 0.00
#> d[Efalizumab] 0.00 0.00 0.01 0.03 0.10 0.36 0.52 0.00
#> d[Etanercept 25 mg] 0.00 0.01 0.08 0.21 0.39 0.26 0.05 0.00
#> d[Etanercept 50 mg] 0.08 0.30 0.27 0.24 0.09 0.02 0.01 0.00
#> d[Fumaderm] 0.07 0.09 0.10 0.11 0.16 0.19 0.27 0.01
#> d[Infliximab] 0.60 0.18 0.12 0.06 0.02 0.01 0.00 0.00
#> d[Methotrexate] 0.09 0.12 0.15 0.18 0.17 0.15 0.14 0.00
plot(pso_rankprobs)
<- posterior_rank_probs(pso_fit_RE, lower_better = FALSE, cumulative = TRUE))
(pso_cumrankprobs #> p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6] p_rank[7] p_rank[8]
#> d[Supportive care] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 1
#> d[Ciclosporin] 0.17 0.46 0.73 0.90 0.98 1.00 1.00 1
#> d[Efalizumab] 0.00 0.00 0.01 0.03 0.13 0.48 1.00 1
#> d[Etanercept 25 mg] 0.00 0.02 0.10 0.31 0.70 0.95 1.00 1
#> d[Etanercept 50 mg] 0.08 0.38 0.65 0.89 0.97 0.99 1.00 1
#> d[Fumaderm] 0.07 0.15 0.25 0.37 0.53 0.72 0.99 1
#> d[Infliximab] 0.60 0.78 0.90 0.96 0.99 1.00 1.00 1
#> d[Methotrexate] 0.09 0.21 0.36 0.54 0.71 0.85 1.00 1
plot(pso_cumrankprobs)