Example: Plaque psoriasis HTA report

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:

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:

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.

Setting up the network

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.

pso_net <- set_agd_arm(hta_psoriasis, 
                       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
  ggplot2::theme(legend.box.spacing = ggplot2::unit(0.75, "in"),
                 plot.margin = ggplot2::margin(0.1, 0, 0.1, 0.75, "in"))

Meta-analysis models

We fit both fixed effect (FE) and random effects (RE) models.

Fixed effect meta-analysis

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.

pso_fit_FE <- nma(pso_net, 
                  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.35     1.30     1.69     1.91     2.15     2.64  1523    1
#> d[Efalizumab]           1.19    0.00 0.06     1.07     1.15     1.19     1.23     1.30  2250    1
#> d[Etanercept 25 mg]     1.51    0.00 0.10     1.32     1.44     1.51     1.57     1.70  1794    1
#> d[Etanercept 50 mg]     1.92    0.00 0.10     1.72     1.85     1.92     1.99     2.12  2056    1
#> d[Fumaderm]             1.47    0.01 0.49     0.61     1.14     1.44     1.77     2.50  2529    1
#> d[Infliximab]           2.33    0.01 0.27     1.83     2.15     2.33     2.52     2.87  2829    1
#> d[Methotrexate]         1.63    0.01 0.45     0.77     1.33     1.62     1.92     2.53  1719    1
#> lp__                -3405.15    0.10 3.60 -3412.90 -3407.43 -3404.79 -3402.53 -3399.25  1327    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.74     0.76     0.78     0.82  5385    1
#> cc[PASI90]              1.56    0.00 0.05     1.46     1.53     1.56     1.60     1.67  5914    1
#> 
#> Samples were drawn using NUTS(diag_e) at Tue May 23 12:01:09 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).

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")

Random effects meta-analysis

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

pso_fit_RE <- nma(pso_net, 
                  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.04    0.01 0.44     1.30     1.74     2.00     2.29     2.98  3155    1
#> d[Efalizumab]           1.19    0.00 0.18     0.80     1.10     1.19     1.28     1.58  4356    1
#> d[Etanercept 25 mg]     1.53    0.00 0.25     1.00     1.40     1.52     1.65     2.06  4941    1
#> d[Etanercept 50 mg]     1.93    0.00 0.28     1.33     1.79     1.93     2.07     2.52  4523    1
#> d[Fumaderm]             1.50    0.01 0.64     0.29     1.09     1.47     1.88     2.82  7767    1
#> d[Infliximab]           2.31    0.00 0.39     1.54     2.08     2.31     2.55     3.09  7015    1
#> d[Methotrexate]         1.73    0.01 0.66     0.57     1.30     1.69     2.09     3.15  4093    1
#> lp__                -3410.42    0.19 6.73 -3424.09 -3414.94 -3410.10 -3405.61 -3398.14  1245    1
#> tau                     0.32    0.01 0.23     0.02     0.16     0.27     0.44     0.88   920    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.74     0.76     0.78     0.82 15406    1
#> cc[PASI90]              1.56    0.00 0.05     1.46     1.53     1.56     1.60     1.67 18542    1
#> 
#> Samples were drawn using NUTS(diag_e) at Tue May 23 12:03:22 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(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 comparison

Model fit can be checked using the dic() function:

(dic_FE <- dic(pso_fit_FE))
#> Residual deviance: 74.9 (on 58 data points)
#>                pD: 25.5
#>               DIC: 100.4

(dic_RE <- dic(pso_fit_RE))
#> Residual deviance: 62.5 (on 58 data points)
#>                pD: 33.5
#>               DIC: 96

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.

Further results

Predicted probabilities of response

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).

pred_FE <- predict(pso_fit_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.00
#> pred[Supportive care, PASI75]  0.03 0.01 0.02 0.03 0.03 0.04  0.05     3807     3516 1.00
#> pred[Supportive care, PASI90]  0.00 0.00 0.00 0.00 0.00 0.00  0.01     4068     3394 1.00
#> pred[Ciclosporin, PASI50]      0.78 0.10 0.57 0.72 0.79 0.86  0.94     1601     1559 1.01
#> pred[Ciclosporin, PASI75]      0.53 0.13 0.28 0.43 0.52 0.62  0.80     1597     1637 1.01
#> pred[Ciclosporin, PASI90]      0.25 0.11 0.08 0.16 0.23 0.31  0.51     1620     1581 1.00
#> pred[Efalizumab, PASI50]       0.54 0.04 0.45 0.51 0.54 0.56  0.62     2923     3442 1.00
#> pred[Efalizumab, PASI75]       0.25 0.04 0.19 0.23 0.25 0.28  0.33     3036     3325 1.00
#> pred[Efalizumab, PASI90]       0.07 0.02 0.04 0.06 0.07 0.08  0.11     3343     3688 1.00
#> pred[Etanercept 25 mg, PASI50] 0.66 0.05 0.56 0.63 0.66 0.69  0.75     2119     3171 1.00
#> pred[Etanercept 25 mg, PASI75] 0.37 0.05 0.27 0.33 0.36 0.40  0.47     2154     3381 1.00
#> pred[Etanercept 25 mg, PASI90] 0.13 0.03 0.08 0.11 0.12 0.15  0.19     2243     3551 1.00
#> pred[Etanercept 50 mg, PASI50] 0.79 0.04 0.71 0.77 0.79 0.82  0.86     2311     3362 1.00
#> pred[Etanercept 50 mg, PASI75] 0.53 0.05 0.42 0.49 0.53 0.56  0.63     2320     3169 1.00
#> pred[Etanercept 50 mg, PASI90] 0.23 0.04 0.16 0.20 0.23 0.26  0.32     2443     3047 1.00
#> pred[Fumaderm, PASI50]         0.63 0.16 0.31 0.51 0.64 0.75  0.92     2891     2000 1.00
#> pred[Fumaderm, PASI75]         0.36 0.17 0.10 0.23 0.34 0.47  0.75     2895     2022 1.00
#> pred[Fumaderm, PASI90]         0.14 0.11 0.02 0.06 0.11 0.19  0.45     2900     2149 1.00
#> pred[Infliximab, PASI50]       0.88 0.05 0.76 0.85 0.89 0.92  0.97     2936     2860 1.00
#> pred[Infliximab, PASI75]       0.68 0.10 0.48 0.61 0.68 0.75  0.86     2991     2924 1.00
#> pred[Infliximab, PASI90]       0.38 0.11 0.20 0.30 0.37 0.45  0.60     2994     2911 1.00
#> pred[Methotrexate, PASI50]     0.69 0.15 0.36 0.59 0.70 0.80  0.93     1769     2043 1.00
#> pred[Methotrexate, PASI75]     0.42 0.16 0.13 0.30 0.41 0.53  0.76     1760     2038 1.00
#> pred[Methotrexate, PASI90]     0.17 0.11 0.03 0.09 0.15 0.23  0.46     1769     1966 1.00
plot(pred_FE)

pred_RE <- predict(pso_fit_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    10253    10100    1
#> pred[Supportive care, PASI90]  0.00 0.00 0.00 0.00 0.00 0.00  0.01    11065     9464    1
#> pred[Ciclosporin, PASI50]      0.80 0.11 0.57 0.74 0.82 0.89  0.97     3820     3219    1
#> pred[Ciclosporin, PASI75]      0.56 0.16 0.28 0.45 0.56 0.67  0.87     3838     3154    1
#> pred[Ciclosporin, PASI90]      0.28 0.15 0.08 0.18 0.26 0.36  0.64     3897     3202    1
#> pred[Efalizumab, PASI50]       0.54 0.08 0.37 0.49 0.54 0.58  0.70     5661     3984    1
#> pred[Efalizumab, PASI75]       0.26 0.07 0.14 0.22 0.25 0.29  0.40     5716     3974    1
#> pred[Efalizumab, PASI90]       0.08 0.03 0.03 0.06 0.07 0.09  0.15     5834     3972    1
#> pred[Etanercept 25 mg, PASI50] 0.66 0.09 0.46 0.61 0.66 0.72  0.84     5913     4067    1
#> pred[Etanercept 25 mg, PASI75] 0.38 0.10 0.19 0.32 0.37 0.43  0.59     5965     4186    1
#> pred[Etanercept 25 mg, PASI90] 0.14 0.06 0.05 0.10 0.13 0.16  0.29     6104     4392    1
#> pred[Etanercept 50 mg, PASI50] 0.79 0.08 0.59 0.75 0.80 0.84  0.93     5405     3886    1
#> pred[Etanercept 50 mg, PASI75] 0.53 0.11 0.29 0.47 0.53 0.59  0.76     5412     4016    1
#> pred[Etanercept 50 mg, PASI90] 0.24 0.09 0.09 0.19 0.23 0.28  0.46     5540     4115    1
#> pred[Fumaderm, PASI50]         0.63 0.20 0.21 0.49 0.65 0.78  0.96     8267     5197    1
#> pred[Fumaderm, PASI75]         0.38 0.21 0.06 0.22 0.35 0.51  0.84     8338     5146    1
#> pred[Fumaderm, PASI90]         0.16 0.15 0.01 0.06 0.12 0.22  0.57     8273     5428    1
#> pred[Infliximab, PASI50]       0.87 0.08 0.66 0.83 0.89 0.93  0.98     7628     4995    1
#> pred[Infliximab, PASI75]       0.66 0.13 0.37 0.58 0.68 0.76  0.89     7639     5157    1
#> pred[Infliximab, PASI90]       0.37 0.14 0.13 0.28 0.36 0.46  0.67     7644     5164    1
#> pred[Methotrexate, PASI50]     0.70 0.18 0.30 0.58 0.72 0.84  0.98     4827     3070    1
#> pred[Methotrexate, PASI75]     0.45 0.21 0.10 0.29 0.43 0.60  0.91     4850     3085    1
#> pred[Methotrexate, PASI90]     0.21 0.18 0.02 0.09 0.16 0.29  0.70     4889     3149    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).

pred_FE_beta <- predict(pso_fit_FE, 
                        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     4082     4057    1
#> pred[Supportive care, PASI90]  0.00 0.00 0.00 0.00 0.00 0.00  0.01     4371     3724    1
#> pred[Ciclosporin, PASI50]      0.78 0.10 0.58 0.72 0.79 0.85  0.94     1658     1987    1
#> pred[Ciclosporin, PASI75]      0.53 0.13 0.29 0.43 0.52 0.62  0.79     1658     1965    1
#> pred[Ciclosporin, PASI90]      0.24 0.11 0.08 0.16 0.23 0.31  0.50     1684     1926    1
#> pred[Efalizumab, PASI50]       0.54 0.04 0.46 0.51 0.54 0.56  0.61     3290     3569    1
#> pred[Efalizumab, PASI75]       0.25 0.03 0.19 0.23 0.25 0.27  0.32     3323     3713    1
#> pred[Efalizumab, PASI90]       0.07 0.01 0.05 0.06 0.07 0.08  0.10     3473     3569    1
#> pred[Etanercept 25 mg, PASI50] 0.66 0.05 0.57 0.63 0.66 0.69  0.74     2500     3212    1
#> pred[Etanercept 25 mg, PASI75] 0.37 0.05 0.28 0.33 0.37 0.40  0.46     2508     3037    1
#> pred[Etanercept 25 mg, PASI90] 0.13 0.03 0.08 0.11 0.12 0.14  0.18     2557     3604    1
#> pred[Etanercept 50 mg, PASI50] 0.79 0.04 0.72 0.77 0.79 0.82  0.86     2768     3281    1
#> pred[Etanercept 50 mg, PASI75] 0.53 0.05 0.43 0.49 0.53 0.56  0.62     2737     3331    1
#> pred[Etanercept 50 mg, PASI90] 0.23 0.04 0.16 0.20 0.23 0.26  0.31     2758     3251    1
#> pred[Fumaderm, PASI50]         0.63 0.16 0.31 0.52 0.63 0.75  0.93     2808     2212    1
#> pred[Fumaderm, PASI75]         0.36 0.17 0.11 0.24 0.34 0.47  0.75     2808     2332    1
#> pred[Fumaderm, PASI90]         0.14 0.11 0.02 0.06 0.11 0.19  0.45     2809     2318    1
#> pred[Infliximab, PASI50]       0.88 0.05 0.76 0.85 0.89 0.92  0.96     2954     2774    1
#> pred[Infliximab, PASI75]       0.68 0.10 0.48 0.61 0.68 0.75  0.85     3005     2828    1
#> pred[Infliximab, PASI90]       0.38 0.10 0.19 0.30 0.37 0.44  0.59     3001     2788    1
#> pred[Methotrexate, PASI50]     0.69 0.15 0.37 0.59 0.70 0.80  0.93     1840     2180    1
#> pred[Methotrexate, PASI75]     0.42 0.16 0.14 0.30 0.41 0.53  0.75     1831     2092    1
#> pred[Methotrexate, PASI90]     0.17 0.11 0.03 0.09 0.15 0.23  0.45     1839     2194    1
plot(pred_FE_beta)

pred_RE_beta <- predict(pso_fit_RE, 
                        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    10425     9634    1
#> pred[Supportive care, PASI90]  0.00 0.00 0.00 0.00 0.00 0.00  0.01    11490     9632    1
#> pred[Ciclosporin, PASI50]      0.80 0.11 0.57 0.74 0.82 0.88  0.97     3902     3056    1
#> pred[Ciclosporin, PASI75]      0.56 0.15 0.28 0.45 0.56 0.67  0.87     3924     2995    1
#> pred[Ciclosporin, PASI90]      0.28 0.14 0.08 0.18 0.26 0.36  0.64     3978     3263    1
#> pred[Efalizumab, PASI50]       0.54 0.08 0.38 0.49 0.54 0.58  0.70     5446     4339    1
#> pred[Efalizumab, PASI75]       0.26 0.06 0.14 0.22 0.25 0.29  0.40     5511     4216    1
#> pred[Efalizumab, PASI90]       0.07 0.03 0.03 0.06 0.07 0.09  0.15     5635     4455    1
#> pred[Etanercept 25 mg, PASI50] 0.66 0.09 0.46 0.61 0.66 0.72  0.83     6179     4147    1
#> pred[Etanercept 25 mg, PASI75] 0.38 0.10 0.20 0.32 0.37 0.43  0.59     6233     4095    1
#> pred[Etanercept 25 mg, PASI90] 0.14 0.06 0.05 0.10 0.13 0.16  0.28     6367     4196    1
#> pred[Etanercept 50 mg, PASI50] 0.79 0.08 0.58 0.75 0.80 0.84  0.93     5479     3544    1
#> pred[Etanercept 50 mg, PASI75] 0.53 0.11 0.29 0.47 0.53 0.59  0.76     5476     3674    1
#> pred[Etanercept 50 mg, PASI90] 0.24 0.09 0.09 0.19 0.23 0.28  0.45     5580     3691    1
#> pred[Fumaderm, PASI50]         0.63 0.20 0.21 0.49 0.65 0.79  0.96     8326     5222    1
#> pred[Fumaderm, PASI75]         0.38 0.21 0.06 0.22 0.35 0.51  0.83     8392     5294    1
#> pred[Fumaderm, PASI90]         0.16 0.15 0.01 0.06 0.12 0.22  0.56     8347     5384    1
#> pred[Infliximab, PASI50]       0.87 0.08 0.67 0.84 0.89 0.93  0.98     7664     5059    1
#> pred[Infliximab, PASI75]       0.66 0.13 0.37 0.59 0.67 0.76  0.89     7669     5301    1
#> pred[Infliximab, PASI90]       0.37 0.14 0.13 0.28 0.36 0.46  0.67     7677     5326    1
#> pred[Methotrexate, PASI50]     0.70 0.18 0.30 0.58 0.72 0.84  0.98     4891     3247    1
#> pred[Methotrexate, PASI75]     0.45 0.21 0.10 0.29 0.43 0.60  0.91     4913     3171    1
#> pred[Methotrexate, PASI90]     0.21 0.18 0.02 0.09 0.16 0.29  0.70     4951     3248    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\).

Ranks and rank probabilities

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).

(pso_ranks <- posterior_ranks(pso_fit_RE, lower_better = FALSE))
#>                        mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
#> rank[Supportive care]  7.99 0.12    8   8   8   8     8     5976       NA    1
#> rank[Ciclosporin]      2.76 1.28    1   2   3   4     5     7445     7599    1
#> rank[Efalizumab]       6.33 0.82    4   6   7   7     7     5703     5550    1
#> rank[Etanercept 25 mg] 4.94 1.07    3   4   5   6     7     7249     5668    1
#> rank[Etanercept 50 mg] 3.05 1.24    1   2   3   4     6     5641     5104    1
#> rank[Fumaderm]         4.86 1.97    1   3   5   7     7     7970     5590    1
#> rank[Infliximab]       1.84 1.21    1   1   1   2     5     3683     4661    1
#> rank[Methotrexate]     4.23 1.88    1   3   4   6     7     6051     5623    1
plot(pso_ranks)

(pso_rankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE))
#>                     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.16      0.08      0.02      0.00      0.00
#> d[Efalizumab]            0.00      0.00      0.01      0.02      0.10      0.36      0.51      0.00
#> d[Etanercept 25 mg]      0.00      0.01      0.08      0.21      0.38      0.26      0.05      0.00
#> d[Etanercept 50 mg]      0.08      0.30      0.26      0.24      0.09      0.02      0.01      0.00
#> d[Fumaderm]              0.08      0.09      0.10      0.12      0.16      0.18      0.27      0.01
#> d[Infliximab]            0.57      0.20      0.13      0.06      0.03      0.01      0.00      0.00
#> d[Methotrexate]          0.10      0.12      0.15      0.18      0.16      0.14      0.15      0.00
plot(pso_rankprobs)

(pso_cumrankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE, cumulative = TRUE))
#>                     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.49      1.00         1
#> d[Etanercept 25 mg]      0.00      0.02      0.10      0.31      0.69      0.95      1.00         1
#> d[Etanercept 50 mg]      0.08      0.38      0.64      0.88      0.97      0.99      1.00         1
#> d[Fumaderm]              0.08      0.17      0.26      0.38      0.54      0.72      0.99         1
#> d[Infliximab]            0.57      0.76      0.89      0.95      0.98      1.00      1.00         1
#> d[Methotrexate]          0.10      0.21      0.37      0.54      0.71      0.85      1.00         1
plot(pso_cumrankprobs)