Highlights

2 Treatment hierarchies in standard NMA models

In a standard NMA with no TCIs, the parameters representing treatment performance are given by $\psi _{20}^{(1)}, \psi _{30}^{(1)}, \dots , \psi _{G0}^{(1)}$ , where $\psi _{g0}^{(1)}$ is the expected relative effect of treatment g compared to treatment 1, and the index $0$ indicates that the parameter represents the main effect of the treatment and not an interaction, which will become important in the following section. Additionally, G is the number of treatments in the network, and $\psi _{10}^{(1)} = 0$ . We assume that the expected relative effect of treatment g versus h, $\psi _{g0}^{(h)}$ , is given by $\psi _{g0}^{(1)}-\psi _{h0}^{(1)}$ for any $g,h$ in $1, \dots , G$ . This is referred to as the consistency assumption in NMA literature.Reference Ades, Welton, Dias, Phillippo and Caldwell⁸ In the remainder of the article, we suppress the superscript $(1)$ and emphasize that all parameters are relative to treatment 1.

SUCRA values are calculated from the relative effects to obtain a treatment hierarchy that is relevant to the population represented by the NMA. The SUCRA value for treatment g is calculated in a standard Bayesian NMA as

(2.1) $$ \begin{align} \text{SUCRA}(g) &= \frac{\sum_{s = 1}^{G-1} \sum_{r = 1}^s p_{gr}}{G-1} \notag \\ &= \frac{G - \text{E}\left[rank(g)\right]}{G-1}, \end{align} $$

where $p_{gr}$ is the probability that treatment g takes rank r, that is, it has a more favorable response than exactly r treatments, and $\text {E}\left [rank(g)\right ]$ is the average rank of treatment g.Reference Salanti, Ades and Ioannidis^2, Reference Rücker and Schwarzer⁹

SUCRA values can be computed in a standard Bayesian NMA using samples from the posterior distribution. Assume that we have a sample of size M from the posterior distribution of $\psi _{10} = 0, \psi _{20}, \dots , \psi _{G0}$ , and let $\psi _{g0}{}_{(m)}$ represent the value of $\psi _{g0}$ in the $m^{\text {th}}$ posterior sample. Note that $\psi _{10} = 0$ by definition for every posterior sample. Then we can calculate SUCRA values as follows:

1. 1. For each sample $m = 1, \dots , M$ and each treatment $g = 1, \dots , G$ , determine the treatment rank according to

   $$ \begin{align*} rank(g)_{(m)} = G - \sum_{h \neq g} I\left\{\psi_{g0}{}_{(m)}> \psi_{h0}{}_{(m)}\right\}, \end{align*} $$
   where $I\{\cdot \} = 1$ when $\cdot $ is true and $0$ otherwise.

2. 2. Calculate the average rank for each treatment,

   $$ \begin{align*} \text{E}\left[rank(g)\right] = \frac{1}{M}\sum_{m = 1}^M rank(g)_{(m)}. \end{align*} $$

3. 3. Calculate $\text {SUCRA}(g)$ using

   $$ \begin{align*} \text{SUCRA}(g) = \frac{G-\text{E}\left[rank(g)\right]}{G-1}. \end{align*} $$

The SUCRA of treatment g in a standard NMA is interpreted as the proportion of competing treatments that g is expected to beat,Reference Salanti, Ades and Ioannidis² or the inverse-scaled average rank of treatment g,Reference Rücker and Schwarzer⁹ when it is employed in the population represented by the NMA. The resulting treatment hierarchy is not personalized to any particular patient and depends on the covariate distribution in the population.Reference Phillippo, Dias and Ades¹⁰

3 Treatment hierarchies in NMA models with TCIs

TCIs can be estimated using aggregate data NMA models with meta-regression and Individual Participant Data (IPD) NMA models, including recently proposed Individualized Treatment Rule (ITR) NMA models.Reference Riley, Lambert and Staessen^11, Reference Shen, Moodie and Golchi¹² In Supplementary Section S1, we describe the suitability of each model type for creating personalized treatment hierarchies. Briefly, we do not recommend using aggregate data NMA meta-regression to create personalized treatment hierarchies due to its risk of ecological bias and confounding.Reference Riley, Lambert and Staessen^11, Reference Donegan, Williamson, D’Alessandro and Smith^13– Reference Schmid, Stark, Berlin, Landais and Lau¹⁵ IPD NMA can be used to produce personalized treatment hierarchies when care is taken in specifying the model.Reference Riley, Dias and Donegan^6, Reference Riley, Debray and Fisher^7, Reference Hua, Burke, Crowther, Ensor, Tudur Smith and Riley¹⁶

In an NMA model that includes TCIs, the interpretation of the parameter $\psi _{g0}$ does not correspond to the expected relative effect on average in the population. Instead, $\psi _{g0}$ corresponds to the expected relative effect of treatment g compared to treatment 1 in the reference group. To illustrate, suppose there are Q effect-modifying covariates included as TCIs in an NMA model. The expected relative effect of treatment g versus treatment 1 for an individual with covariate values $X_1 = x_1, X_2 = x_2, \dots , X_Q = x_Q$ is given by

(3.1) $$ \begin{align} \psi_{g0} + \psi_{g1} x_1 + \psi_{g2} x_2 + \dots + \psi_{gQ} x_Q, \end{align} $$

where $\psi _{g1}, \dots , \psi _{gQ}$ are the coefficients representing the interaction between each covariate and treatment g compared to treatment 1, that is, $\psi _{gq}, q> 0$ is the change in the relative effect of treatment g compared to treatment 1 for a one unit increase in covariate q, and $\psi _{g0}$ is the expected relative effect when all $x_1, \dots , x_Q$ equal 0.

Now, consider creating a treatment hierarchy from a model that includes TCIs. If we were to use only $\psi _{g0}$ , $g = 1, \dots , G$ to create the hierarchy for such a model, the hierarchy would apply to patients whose covariate values are all equal to zero. The hierarchy will be meaningless for patients with a different covariate profile. Clearly, for a treatment hierarchy to be meaningful in an NMA with TCIs, it should be intentionally created for a particular set of covariate values.

Let $\text {E}\left [rank(g)\mid x_{1}, \dots , x_Q\right ]$ be the average rank of treatment g given covariates $x_{1}, \dots , x_Q$ . The average ranks are determined using the expected relative effects defined in equation (3.1). They can be calculated using samples from the posterior of $\psi _{gq}$ , $g = 1, \dots $ , G and $q = 0$ , $\dots $ , Q

1. 1. For each sample $m = 1, \dots , M$ and each treatment $g = 1, \dots , G$ , determine the expected relative treatment effect given $x_{1}, \dots , x_Q$ via

   $$ \begin{align*} RE_{1g}{}_{(m)}(x_{1}, \dots, x_Q) = \psi_{g0}{}_{(m)} + \psi_{g1}{}_{(m)} x_1 + \psi_{g2}{}_{(m)} x_2 + \dots + \psi_{gQ}{}_{(m)} x_Q \quad \text{for } g = 2, ..., G, \end{align*} $$
   and $RE_{11} {}_{(m)}(x_{1}, \dots , x_Q) = 0$ , where $RE_{1g} {}_{(m)}(x_{1}, \dots , x_Q)$ represents the expected relative effect of treatment g versus 1 given covariates $x_{1}, \dots , x_Q$ in sample m.

2. 2. For each sample $m = 1, \dots , M$ and each treatment $g = 1, \dots , G$ , determine the treatment rank according to

   $$ \begin{align*} rank(g)_{(m)} = G - \sum_{h \neq g} I\left\{RE_{1g} {}_{(m)}(x_{1}, \dots, x_Q)> RE_{1h} {}_{(m)}(x_{1}, \dots, x_Q)\right\}. \end{align*} $$

3. 3. Calculate the average rank for each treatment given the covariates,

   $$ \begin{align*} \text{E}\left[rank(g)\mid x_{1}, \dots, x_Q\right] = \frac{1}{M}\sum_{m = 1}^M rank(g)_{(m)}. \end{align*} $$

4. 4. Calculate $\text {SUCRA}(g)$ using

   $$ \begin{align*} \text{SUCRA}(g) = \frac{G-\text{E}\left[rank(g)\mid x_{1}, \dots, x_Q\right]}{G-1}. \end{align*} $$

The treatment hierarchy produced in this way is personalized to a patient with the covariate profile defined by $x_{1}, \dots , x_Q$ .

4 Example: Pharmacological treatments for MDD

We illustrate the use and impact of personalized treatment hierarchies using data from three studies investigating treatment for MDD. A systematic review process was not followed to select the studies; rather, the studies were selected based on the availability of IPD. Details on the study designs are given in Supplementary Section S2. The data provide information on six pharmacological treatments: Bupropion, Citalopram + Bupropion, Citalopram + Buspirone, Escitalopram, Sertraline, and Venlafaxine. The response is defined as the negative of the 17-question Hamilton Rating Scale for Depression (HRSD-17), so that larger values are preferred. The following covariates, considered to be potential effect modifiers for depression, were included as TCIs in the analysis: age, sex, marital status, years of education, employment status, household size, age of depression onset, number of depressive episodes, chronicity of current MDD episode, and HRSD-17 before receiving any treatment in the study.Reference Kessler, van and Wardenaar^17, Reference Perlman, Benrimoh and Israel¹⁸ We fit a two-stage common effect Bayesian ITR NMA model with IPD to estimate expected relative treatment effects in the MDD network,Reference Shen, Moodie and Golchi¹² described in detail in Supplementary Section S3.

We computed personalized treatment hierarchies using SUCRA for two hypothetical patients, Patient A and Patient B, each with unique covariate profiles: Patient A is employed while Patient B is unemployed, Patient A has experienced three or fewer depressive episodes while Patient B has experienced more than three, Patient A is male and Patient B is not, and Patient A lives in a household with two more members than Patient B. The specific covariate values for each patient are described in Supplementary Section S3.3.

4.1 Results

A plot showing point estimates and 95% credible intervals for treatment effects and TCIs is shown in Supplementary Section S3.4. All credible intervals contain zero. This is not too surprising, since only three studies were included in the analysis. In Figure 1, we show the posterior distribution of the expected relative effects as defined in equation (3.1) of each treatment compared to Sertraline for Patients A and B. It is clear that there are differences in the posterior densities between the patients. The SUCRA values for Patients A and B reflect the differences in the posteriors, and are compared in Table 1. For example, Citalopram + Bupropion is ranked first for Patient A, while it is ranked only fourth out of the six treatments for Patient B. Table 1 also compares the SUCRA and treatment rankings that would result from using a naive method that only considered the main treatment effects, $\psi _{g0}$ . The naive method gives the same hierarchy for both patients, and this differs from both personalized hierarchies, particularly for Patient B, who deviates more from the reference profile.

Table 1

SUCRA values and ranking of each treatment using a naive nonpersonalized approach, and using a personalized approach for Patient A and Patient B

Figure 1

Posterior distributions of the relative effects of each treatment for Patients A and B.