A dose–effect network meta-analysis model with application in antidepressants using restricted cubic splines

Abstract

Network meta-analysis has been used to answer a range of clinical questions about the preferred intervention for a given condition. Although the effectiveness and safety of pharmacological agents depend on the dose administered, network meta-analysis applications typically ignore the role that drugs dosage plays in the results. This leads to more heterogeneity in the network. In this paper, we present a suite of network meta-analysis models that incorporate the dose–effect relationship using restricted cubic splines. We extend existing models into a dose–effect network meta-regression to account for study-level covariates and for groups of agents in a class-effect dose–effect network meta-analysis model. We apply our models to a network of aggregate data about the efficacy of 21 antidepressants and placebo for depression. We find that all antidepressants are more efficacious than placebo after a certain dose. Also, we identify the dose level at which each antidepressant's effect exceeds that of placebo and estimate the dose beyond which the effect of antidepressants no longer increases. When covariates were introduced to the model, we find that studies with small sample size tend to exaggerate antidepressants efficacy for several of the drugs. Our dose–effect network meta-analysis model with restricted cubic splines provides a flexible approach to modelling the dose–effect relationship in multiple interventions. Decision-makers can use our model to inform treatment choice.

Keywords

Evidence synthesis multiple treatments splines dose-response meta-regression

1 Introduction

Network meta-analysis (NMA) is a technique commonly used to simultaneously compare multiple agents.^1–3 Although the comparison between pharmacological agents is important, in practice clinicians always prescribe drugs at a particular dose, informed by the market authorisation, licensing of the product, dose–effect studies, and their experience. It is therefore important to know not only which pharmacological agents are preferable but how their advantage depends on the dose.

Health technology assessment agencies make recommendations that should, and sometimes do, specify the recommended dose range for several competing pharmacological agents. However, without a unified methodological approach to infer the relative effects of agent–dose combinations, contradictory information might be made available. For example, the guidelines produced by the National Institute of Health and Care Excellence state that no dose-dependency has been established within the therapeutic range of selective serotonin-reuptake inhibitors (SSRI) when treating people diagnosed with major depression, whereas the American Psychiatric Association guideline recommends titration up to the maximum tolerated dose.

In NMA, the first and often most challenging step is the definition of the nodes in the network with respect to the combination of agents and dose. When pharmacological agents are compared, an important decision faced early on is whether the dose of each agent is of interest, and consequently, whether the definition of each node incvolves the dose of the agent or not.

There are three main options when it comes to dealing with the dose of pharmacological interventions in NMA. Frequently, information about the dose is ignored and focus is placed only on the relative effects between agents (e.g. Cipriani et al.⁴). This approach may result in a network with increased heterogeneity and inconsistency. At the other end of the spectrum, one can consider each agent–dose combination as a different treatment that defines a different node in the network.⁵ This detailed and larger network will inevitably be at best sparse or even disconnected. A compromise is to model the dose–effect relationship for each agent by extending the dose–effect meta-analysis models.^6–8

The dose–effect relationship expresses the change in effect over different doses. In pairwise meta-analysis, the dose–effect curves are synthesised across studies. Such analyses can be conducted using the two-stage or one-stage methods in a frequentist^9,10 or Bayesian setting.¹¹ In NMA, the linear dose–effect model has been implemented,⁶ which, however, poorly reflects the natural dose–effect dynamics.¹² Del Giovane et al.⁷ addressed that by either considering an exchangeable effect for the different doses of a certain agent or assuming the dose–effect relationship as a monotonic, linear, or random walk. More recently, Mawdsley et al.⁸ have extended NMA to incorporate the Emax dose–effect model which is commonly used in pharmacometrics when determining the optimal dose. In clinical practice and for decision making, more flexibility in the assumed dose–effect shapes is desirable to better reflect a range of possible biological mechanisms of the various pharmacological agents.

With this paper, we aim to contribute to the growing literature about dose–effect models by describing a generic and flexible dose–effect NMA (DE-NMA) model with restricted cubic splines (RCS). Recent simulations showed that the RCS successfully capture a large range of functional shapes.¹¹ As residual heterogeneity and inconsistency (beyond what can be explained by different dosages) can occur, we extend the model into a dose–effect network meta-regression (DE-NMR) by incorporating study-level covariates.

The article is structured as follows. First, we present our motivating example. In Section 3, we present the DE-NMA model and two extensions: the DE-NMR and a DE-NMA that includes class effects. Next, we apply the models to the antidepressants network, and we then present the results. In Section 4, we discuss the strengths and limitations of the models, and we discuss other methods to estimate the dose-response shape, such as fractional polynomials.

The analysis code is implemented using Just Another Gibbs Sampler (JAGS)¹³ and R,¹⁴ and it is made available at Zenodo.¹⁵

2 Example: comparing the efficacy of 21 antidepressants

We illustrate our different models using a network of double-blind fixed-dose randomised controlled trials (RCT) that compare antidepressants for depression (see Figure 1(a) and Appendix Figure 1). The primary outcome is efficacy measured as the total number of patients who had more than 50% reduction in symptoms (response rate).¹⁶ The participants of the included studies were adults diagnosed with unipolar major depressive disorder. The dataset is a superset of one used to compare 21 antidepressants and placebo according to their efficacy, acceptability, and safety.⁴ In that NMA, Cipriani et al.⁴ synthesised only arms with agents administered at approved doses (as fixed or flexible schedule), while we included all trial arms regardless of the dosage. More details about inclusion criteria, the search strategy, data extraction, and risk of bias in these studies can be found in Cipriani et al.⁴

Figure 1.

(a) Network meta-analysis of 21 antidepressants and placebo. The width of the lines is proportional to the number of trials comparing each pair of agents. This plot was produced using the plot() function from the R package MBNMAdose. (b) The dose distribution for the 21 antidepressants.

Our dataset includes 170 RCTs comparing 21 antidepressants with placebo or another active treatment. Trazodone is excluded from the primary analysis (because only one dose level is examined in the included studies), yet included in the class-effect model (it belongs with nefazodone in the same class, serotonin antagonist and reuptake inhibitor (SARIs)). The trials report 457 different fixed dose-per-drug treatments and include 54,048 participants. In Appendix Table 1, we summarise the number of events, the sample size, the number of studies, the number of different doses, and the class for each drug. We present the distribution of observed doses per drug in Figure 1(b).

A subset of our data (only SSRIs except for fluvoxamine) has been previously analysed using dose–effect pairwise meta-analysis, thus ignoring the differences between the individual drugs (using frequentist¹⁷ or Bayesian dose–effect meta-analysis¹¹).

3 Methods

We first present the DE-NMA model, and then we extend it by incorporating covariates or by assuming class effects between the exposures. As most studies in any NMA are RCTs, we assume the case where each arm of a trial has been randomised to an agent at a different dose. We also present the model assuming a dichotomous outcome using RCS for the association between dose and effects. However, the model can be easily adapted for any assumed shape (e.g. linear, quadratic, etc.) and any type of outcome.

3.1 Notation

Table 1 summarises the notation we used. Suppose we compared K agents ( $k = 1, \dots, K$ ) in ns studies ( $i = 1, \dots, n s$ ) that report the dose level j. For each dose, $x_{i j k}$ , we observe the number of events $r_{i j k}$ and the sample size $n_{i j k}$ (dichotomous outcome). Additionally, we have information on a study-level covariate $Z_{i}$ . In the class-effect model, the c index refers to the class of the agent. Note that we differentiated between agent and treatment when the latter refers to the given dose of a certain agent.

Table 1.
Notations for dose–effect network meta-analysis (DE-NMA).

$i = 1, \dots, n s$ Study id

$j$ Index for the dose levels in study $i$

$k = 1, \dots, K$ Agent

$c = 1, \dots, C$ Exposure clusters

$p = 1, \dots, P$ Number of dose transformations associated with the dose-effect shape. For a linear shape $p = 1$ and for quadratic and restricted cubic splines $p = 2$

$x_{i j k}$ The jth dose in study i for agent k

$x_{d k}$ Minimal dose d for agent k

$Z_{i}$ Covariate in study i

$r_{i j k}$ Number of events in dose j within study i for agent k

$n_{i j}$ Sample size in dose j within study i for agent k

3.2 DE-NMA model with a placebo arm

We defined the DE-NMA model as an extension of the standard NMA. We describe it as a hierarchical model with three layers; we first estimated the dose–effect association within each study, and then we synthesised the shapes across studies and agents. For simplicity, we present the model assuming that the network includes placebo.

3.2.1 Dose–effect model within each study

For ease of understanding and notation, we begin our description of a dose–effect model for a network of trials that have ‘placebo’ as a common comparator. We relax this initial assumption at the end of Section 3.2.1, when we describe a dose–effect model for a network of studies with different controls.

Let us assume that within each study i, the number of events follows a binomial distribution

r_{i j k} \sim Binomial (p_{i j k}, n_{i j k})

with

p_{i j k}

being the probability of an event to occur. We choose a transformation of these probabilities based on the measure of the relative effect that we are interested in. We set the transformation to the logit transformation for odds ratio (OR)

logit (p_{i j k}) = {\begin{matrix} u_{i}, & placebo \\ u_{i} + δ_{i j k}, & active agent \end{matrix}

u_{i}

is the log-odds of the event on the placebo arm in study i. The term

δ_{i j k}

denotes the underlying parameter for the effect of agent k in study i at dose

x_{i j k}

(dose level j). It is the effect of agent k in study i at dose

x_{i j k}

relative to placebo (or the minimum dose in the study i; see the end of the section). If the log function instead of the logit is used to transform the probabilities, the model will estimate risk ratios instead of OR.

The parameter, $δ_{i j k}$ , can be modelled then assuming a common- or exchangeable-effect model, see Table 2. For the common-effect model, the underlying true effect is assumed to be equal in all studies, so we set

δ_{i j k} = Δ_{i j k}

Table 2.
List of potential assumptions for the parameters in the dose–effect network meta-analysis (DE-NMA) model.

Assumptions about the effect parameter $δ_{i j k}$

Assumption 1.1 – exchangeable

$δ_{i j k} \sim N (Δ_{i j k}, τ^{2})$

Assumption 1.2 – common

$δ_{i j k} = Δ_{i j k}$

Assumptions about the pth within-study shape parameter $β_{p, i k}$

Assumption 2.1 – exchangeable

$β_{p, i k} \sim N (B_{p, k}, σ_{β, p}^{2})$

Assumption 2.2 – common

$β_{p, i k} = B_{p, k}$

Assumptions about the pth between-agents shape parameter: $B_{p, k}$

Assumption 3.1 – independent

$B_{p, k} = b_{p, k}$

Assumption 3.2 – exchangeable

$B_{p, k} \sim N (b_{p}, σ_{B, p}^{2})$

Assumption 3.3 – common

$B_{p, k} = b_{p}$

Assumptions for the class-effect model

Assumption 3.4 – exchangeable class-effect across agents $k_{c}$ belonging to class c

$B_{p, k_{c}} \sim N (b_{p, c}, σ_{B, p}^{2})$

Assumption 3.5 – common class-effect across agents $k_{c}$ belonging to class c

$B_{p, k_{c}} = b_{p, c}$

For exchangeable-effect model, $δ_{i j k}$ are assumed to come from a common normal distribution with mean $Δ_{i j k}$ and variance $τ^{2}$ ,

δ_{i j k} \sim N (Δ_{i j k}, τ^{2})

The heterogeneity

τ^{2}

reflects between-studies variability, and it is assumed to be independent of the dose and agent. For multi-arm trials with more than one active agent examined, there are more than one

δ_{i j k}

per study, and as they are calculated using the same reference arm, they shall be jointly modelled using a multivariate normal distribution as in standard NMA.²

Το incorporate the dose–effect relationship in the model, we linked the parameters $Δ_{i j k}$ to the transformed doses under an assigned function F, which we will call the dose–effect function:

Δ_{i j k} = F (x_{i j k}; β_{1, i k}, β_{2, i k}, \dots, β_{P, i k})

(1)

The function F can take various forms and the shape is defined by a set of P parameters

β_{1, i k}, \dots, β_{P, i k}

. In addition to that, F can be set differently for each agent k;

F_{k}

. Here we will set F to be an RCS – the same for all agents.

The general form of the RCS with $κ$ knots $t_{1}, \dots, t_{κ}$ is defined as follows:

\begin{aligned} F (x_{i j k}; β_{1, i k}, β_{2, i k}, \dots, β_{P, i k}) \\ = β_{1, i k} x_{i j k} + β_{2, i k} f_{2} (x_{i j k}) + \dots + β_{(κ - 1), i k} f_{(κ - 1)} (x_{i j k}) \end{aligned}

(2)

where for

m = 1, \dots, (κ - 2)

f_{(m + 1)} (x_{i j k}) = (x_{i j k} - t_{m})_{+}^{3} - \frac{t_{κ} - t_{m}}{t_{κ} - t_{k - 1}} (x_{i j k} - t_{κ - 1})_{+}^{3} + \frac{t_{κ - 1} - t_{m}}{t_{κ} - t_{κ - 1}} (x_{i j k} - t_{κ})_{+}^{3} .

with

(x)_{+} = x if x > 0

and 0 otherwise. For more details, see Section 2.4.5 in Harrell.¹⁸

Setting three knots ( $κ = 3$ ) will reduce F in equation (2) into a function with two coefficients. Then the dose–effect relationship becomes expressed by the linear and the spline terms; $Δ_{i j k} = β_{1, i k} x_{i D k} + β_{2, i k} f (x_{i j k})$ . We use three knots for the remainder of the paper. A discussion about selecting the number of knots and their location can be found elsewhere.^11,18

When the study i does not have a placebo arm, we can choose an agent R at the minimum dose level r as the study-specific reference treatment. Then, the relative treatment effect $Δ_{i (j r) (k R)}$ refers to the effect of agent k at dose level j versus agent R at dose level r; it is modelled as

Δ_{i (j r) (k R)} = F (x_{i j k}; β_{1, i k}, β_{2, i k}, \dots, β_{P, i k}) - F (x_{i r R}; β_{1, i R}, β_{2, i R}, \dots, β_{P, i R})

3.2.2 Dose–effect model across studies and agents

To synthesise the dose–effect parameters $β_{p, i k}$ across studies, we employed the following assumptions (see Table 2). We can assume each agent-specific pth shape parameter $β_{p, i k}$ to be either exchangeable $β_{p, i k} \sim N (B_{p, k}, σ_{β, p}^{2})$ (Assumption 2.1) or equal $β_{p, i k} = B_{p, k}$ (Assumption 2.2) across studies. We can simplify Assumption 2.1 by setting a common shape heterogeneity $σ_{β, p} = σ_{β}$ .

Across agents, we can relate the shape parameters based on three possible assumptions (see Table 2). For agent k, we have a set of P shape parameters $B_{p, k}$ ; it can be either independent $B_{p, k}$ (Assumption 3.1), it can have a common normal distribution $B_{p, k} \sim N (b_{p}, σ_{B, p}^{2})$ (Assumption 3.2), or it can be fixed to a single value as $B_{p, k} = b_{p}$ (Assumption 3.3). The latter assumption requires harmonisation of doses, so all agents’ doses are measured on the same scale.

Let us define $Δ_{. (a c) (A C)}$ as the expectation for the log OR between treatment A at dose $x_{a}$ versus treatment C at dose $x_{c}$ . Now, to estimate the dose–effect curve between the two non-referent agents A and C, we can use consistency equations

\begin{aligned} Δ_{. (a c) (A C)} = & Δ_{. a A} - Δ_{. c C} \\ = & B_{1, A} x_{a} + B_{2, A} f (x_{a}) - [B_{1, C} x_{c} + B_{2, C} f (x_{c})] \end{aligned}

(3)

where

Δ_{. a A}

refers to the study-specific treatment effect of agent A at dose

x_{a}

versus placebo.

3.3 DE-NMR model

We can extend the DE-NMA to DE-NMR by adding dose-covariate interaction terms. Assuming an RCS interaction between the covariate and the dose, the DE-NMA model (in Section 3.2) can be updated to

logit (p_{i j k}) = {\begin{matrix} u_{i} & placebo \\ u_{i} + δ_{i j k} + Z_{i} F (x_{i j k}; γ_{1, i k}, γ_{2, i k}, . . ., γ_{P, i k}) & active agent \end{matrix}

The term

Z_{i}

represents a study-level covariate. The parameters

γ_{1, i k}, γ_{2, i k}, . . ., γ_{P, i k}

are expressing the impact of the dose-covariate interaction effect on the relative treatment effect. In most cases however, a linear interaction term should suffice (and would be estimable) so that

F (x_{i j k}; γ_{1, i k}, γ_{2, i k}, . . ., γ_{P, i k}) = γ_{1, i k} x_{i j k} .

Across studies, we can assume either exchangeable-effect; $γ_{m, i k} \sim N (Γ_{m, k}, τ_{γ}^{2}$ ) or a common-effect model; $γ_{m, i k} = Γ_{m, k}$ . Across agents, we can model $Γ_{m, k}$ under one of the following three alternatives; estimate each one independently; $Γ_{m, k} = g_{m, k}$ , assume exchangeable dose-covariate interaction effect $Γ_{m, k} \sim N (g_{m}, τ_{Γ}^{2})$ , or presume a common dose-covariate interaction effect $Γ_{m, k} = g_{m}$ .

We assume consistency for dose-covariate interaction effects per treatment comparison, that is, for the impact of dose-variable m interaction on the parameter effect between two active agents $k_{1}, k_{2}$ ; it is $Γ_{m, k_{1} v s k_{2}} = Γ_{m, k_{1}} - Γ_{m, k_{2}}$ . This means that when we assume $Γ_{m, k} = g_{m}$ , $Γ_{m, k_{1} v s k_{2}} = 0$ .

3.4 DE-NMA model accounting for clusters

Often it might be desirable to group agents in classes and then estimate the class effect alongside agent effects. The assumptions for the shape parameters behind such a model are added as Assumptions 3.4 and 3.5 in Table 2. Such parameters for agents $k_{c}$ belonging to class c can be assumed either exchangeable $B_{p, k_{c}} \sim N (b_{p, c}, σ_{B, p}^{2})$ or common $B_{p, k_{c}} = b_{p, c}$ . Then the parameters $b_{p, c}$ are estimated independently for each class c.

When classes are considered, the doses of the agents within a given class need to be measured on the same (or equivalent scales) to calculate meaningful class effects. For example, to estimate a dose–effect of all SSRIs, we will need first to transform the dose of each different SSRI into the same fluoxetine-equivalent scale.

3.5 Estimating an absolute mean effect for each agent at each dose level and calculating a treatment hierarchy

With many treatments and doses, results are more easily presented and understood using absolute estimands, such as the response probability $P_{j k}$ for a specific dose j of a certain agent k. In a Bayesian setting, this can be done by combining the estimated dose–effect parameters with the response probability for placebo $P_{0}$ . The latter can be computed outside the DE-NMA model by placing a binomial distribution for the corresponding events $r_{i 0}$ with sample size $n_{i 0}$ and probability of the event to occur in placebo arm $p_{i 0}$

\begin{aligned} r_{i 0} \sim & Binom (p_{i 0}, n_{i 0}) \\ \log it (p_{i 0}) \sim & N (logit (P_{0}), σ_{0}^{2}) \end{aligned}

Next, the predicted probability of the event to occur at dose j and agent k is

P_{j k} * = expit (logit ({\tilde{P}}_{0}) + F ((x_{j k}; {\tilde{B}}_{1, k}), (x_{j k}; {\tilde{B}}_{2, k}), \dots, (x_{j k}; {\tilde{B}}_{p, k})) + \tilde{g} \times co v_{pred} \times x_{j k}),

where the tilde (∼) refers to the posterior of parameters.

These probabilities may be used then to rank the agents according to their efficacy at each dose level. However, to make comparison easy, one might need to transform the doses into a single scale using equivalence formulae, if available.

4 Application in dose–effect of antidepressants

4.1 Implementation of the models and diagnostics

We conducted a DE-NMA under five different model specifications. M1 is the primary dose-effect NMA model, and then we added three DE-NMR models (M2–M4) for covariates, risk of bias (low vs high), study publication year (centred at 2010), and the variance of log OR (to evaluate small study effects). In M5, we accounted for class effects instead of the agent effects as listed in Appendix Table 1. All models employ Assumptions 1.1 and 2.2 (Table 2). M1 to M4 additionally employ Assumption 3.1. We set a common dose-covariate interaction effect across studies and agents in M2 to M4 ( $Γ_{m, k} = g_{m})$ . In M5, class effects are modelled using Assumption 3.5 where all doses are transformed to fluoxetine-equivalent dose using previously established transformation.

We modelled the dose–effect relationship with RCS with three knots. Because agents have different dose ranges, knots are placed for each agent at 25%, 50% and 75% percentiles of the corresponding observed dose range. We investigated the sensitivity of the estimated curve to knots position, only for M1 by placing knots at 10%, 20% and 30% percentiles.

All parameters were estimated using the JAGS program which is implemented via R. We assessed the overall performance of the model using the deviance information criterion (DIC) statistic and leverage plots. The values of DIC can be used to compare between different models but they need to have the same likelihood and data. The model provides the best balance between model fit and complexity when it has the lowest DIC.

We estimated the parameters with Markov Chain Monte Carlo (MCMC) using three chains with $1 \times 10^{4}$ iterations, $4 \times 10^{3}$ burn-in, and a thinning of one. We set a minimally informative prior for the placebo effect $u_{i} \sim N (0, 10^{3})$ and the shape parameters $B_{p, k} \sim N (0, 10^{3})$ , $b_{p} \sim N (0, 10^{3})$ . The two heterogeneity parameters are given a uniform prior $τ, σ_{p, B} \sim Unif (0, 5)$ . For the covariate effect in DE-NMR (M2–M4), we set $g_{m} \sim N (0, 10^{3})$ . For the placebo response model, we placed $logit (P_{0}) \sim N (0, 10^{3})$ and $σ_{0} \sim Unif (0, 5)$ .

We used the rcs function from the rms package to compute the RCS transformations.¹⁹ The codes are available at Zenodo library.¹⁵ We used different numerical and graphical methods (using the coda package²⁰) to investigate the convergence of the MCMC. The results are provided as a posterior median with the 95% credible interval (CrI).

4.2 Results

In Figure 2, we depict the absolute dose–effect relationship for each antidepressant along with the overall placebo effect for M1. The resulting estimates for the linear and spline coefficients are shown in Appendix Table 3–4. The response to placebo is estimated at a mean of 36.2% (95% CrI 34.4%–38.0%) (blue line). All antidepressants are more effective than placebo after some dose level which differed by an agent. However, for some agents, there is a lot of uncertainty particularly at high dose levels (except clomipramine at a low dose level where we have no data). The efficacy initially increased up to a specific dose only to flatten out after a given dose for most agents. For example, the efficacy of duloxetine increased until 75 mg, then it leveled out after that. We identified moderate to small differences in the estimated curves from M1 when we changed knot positions (see Appendix Figure 6). However, the overall conclusions do not change with the change in knot locations. We investigated the convergence of M1 model parameters in Appendix Table 2.

Figure 2.

Dose–effect network meta-analysis summary curve for each antidepressant. The blue line depicts the effect estimated from all placebo arms in the network (36.2%) with its 95% credible region. The red line represents the absolute response to each antidepressant (estimated from model M1) and the shaded area is its 95% credible region.

In M2 to M4 models, we estimate the dose-effect curves assuming three different covariates. In Appendix Figure 3, M2 suggests that studies with high risk of bias (RoB) tend to overestimate on average the efficacy compared to low RoB for some drugs such as bupropion. The average of bupropion efficacy is also more exaggerated in older studies (M3); this is additionally observed for some other antidepressants, see Appendix Figure 2. In M4, the efficacy of most antidepressants is on average higher in small studies (or studies with large variance of logOR) compared to studies with large sample size (Figure 3). In Appendix Table 5, we summarize the findings and the performance for all models.

In Appendix Figure 5, we present the contribution of each observation to $p_{D}$ in y-axis and to ${\bar{D}}_{r e s}$ in x-axis along with the overall model fit measures DIC, $p_{D}$ and ${\bar{D}}_{r e s}$ . The DIC is 790 for the M1, M2, M3 models and it is slightly declined to 789 for M4 model with the variance of logOR as a covariate (Appendix Table 5).

In Figure 4, we show the absolute probabilities under the class-effect model M5. As expected, for classes with many drugs such as SSRIs and serotonin-norepinephrine reuptake inhibitors we gained precision compared to the agent-level models M1.

Figure 3.

Dose-effect network meta-regression summary curve of each one of the 20 antidepressants using the study variance of log odds ratio as a covariate (estimated from model M4). The dose-effect curves are depicted for studies with low variance at 0.027 (green) and with large variance at 0.95 (red). The dotted lines represent the 95% credible interval.

Figure 4.

Dose-effect network meta-analysis summary curve for each of the 9 drug classes (see Appendix Table 1). The blue line depicts the effect estimated from all placebo arms in the network (36.2%) and its 95% credible region. The red line represents the absolute response to each drug class (estimated from model M5) and the shaded area is its 95% credible region.

5 Discussion

We present a dose-effect NMA model to synthesise evidence from trials that compare multiple agents at different dosages. To model the dose–effect relationship, we choose an RCS model to take advantage of their flexibility. We added two extensions to the model: a DE-NMR to account for study-level covariates and groups of agents in a class-effect model. We implemented various DE-NMA models in a network of antidepressants and placebo, and the resulting dose–effect shapes are in line with clinical expectations and previous findings.^4,11,17 Introducing covariates allows us investigating how the dose-effect curve changes at different values of the covariate. These changes are substantial for antidepressants when we added the logOR variance as a covariate. Modelling class effects resulted in more precise estimates of the dose–effect association. We, additionally, identified the specific dose range in which the antidepressant effect exceeds the placebo effect and beyond which dose the effect no longer increases.

Some limitations of the DE-NMA models need to be acknowledged. First, the findings from such analyses can be sensitive to the assumptions about the dose–effect shapes (whether it is an assumed polynomial or splines). Besides a sensitivity analysis, researchers can a priori narrow down the set of possible shapes to the ones that best reflect the known biological behaviour of agents. If needed, the goodness of fit statistics can guide the final choice when enough data is available. When several models provide an equally good fit, Bayesian model averaging can be used. The location of knots in RCS requires particular attention. The estimation of the model could be sensitive to the location of knots, and a sensitivity analysis is recommended to explore any impact on the results.¹¹ Although some researchers argue that the location of knots is not problematic in general,^18,21 we have previously found that positioning the knots at places where shift changes in the effect are expected might be a good strategy.¹¹

Second, there are often very few observations for the same agent to estimate the dose–effect relationship with precision high enough to inform clinical practice. In such cases, the analysis might require considering other sources of information, such as informative priors for the shape, and coefficients of the association based on an external source, such as observational studies, or clinical expertise. Alternatively, we can impose stronger assumptions by borrowing information internally, such as assuming class effects or even exchangeable dose–effect coefficients across all agents. This assumption will improve the parameters’ identifiability, and it also enables us to analyse a disconnected network. This approach, however, requires the doses to be harmonised across the agents and assumes exchangeable dose–effect shapes across agents, which might be difficult to justify in practice.

At present, we only synthesised fixed-dose studies. Studies with a flexible-dose schedule, where the dose is increased up to a maximum targeted level, depending on the patients’ response and acceptability, require special attention. The analysis of post-randomisation dose adjustments requires causal modelling and individual participant data. Synthesising fixed and flexible-dose studies is challenging and results will require careful interpretation. Finally, we did not examine potential inconsistencies in the data; this can be done using newly introduced methods.²²

Technically, dose-response meta-analysis with RCS with three knots requires two studies with at least one of them having three different dose levels. However, issues of precision, model fit, and heterogeneity question the utility of results from such analyses. Depending on the sparseness of the outcome and the complexity of the underlying dose-response shape, a substantial amount of data might be required to obtain useful results from dose-response meta-analyses.

In the present study, we only synthesised fixed-dose randomised studies where all patients in a study arm were prescribed and took the same dose of the same antidepressant. That means that dose is not a ‘patient-level’ characteristic aggregated over the study arm, but an arm-level characteristic. Consequently, aggregation bias is unlikely in the dose-response association with fixed-dose randomised trials. Including studies with a flexible-dose schedule, where the dose is increased up to a maximum targeted level according to the patients’ response and acceptability, warrants special attention. The analysis of post-randomisation adjustments of the dose requires causal modelling and individual participant data. Synthesising fixed and flexible-dose studies is challenging and results need cautious interpretation.

There is a variety of functional forms to model the dose–effect relationship in NMA, such as the Emax model.⁸ The Emax model is widely used in the drug development context where the focus is on studying drug safety and finding optimal doses (e.g. finding the dose at which half of the maximum effect is achieved). In clinical practice, however, the interest is on estimating the dose–effect relationship for the whole dose range. In this context, the parameters of the Emax model are of less interest and the dynamics of the function makes it less likely to portray the underlying true dose–effect association. In contrast, the RCS offer sufficient flexibility to capture the biological behaviour of agents with only a few parameters (only two parameters when we set three knots). This is particularly important in larger dose levels where the efficacy of many pharmacological agents is expected to level out.

A fractional polynomial is another alternative to model the dose–effect relationship. They have been shown to perform well when modelling longitudinal data in NMA²³ but have not been implemented in the DE-NMA context yet. However, they can be less appealing when modelling dose–effect associations. Fractional polynomials are non-local functions which means they can be less efficient in detecting the multiple changes in drug dynamics.²⁴ Although fractional polynomials might be useful in the dose-findings studies where the focus is on the safest dose,²⁵ the RCS might be preferable in (network) meta-analysis contexts where we can benefit from the locality of the RCS to place knots at the expected changing points based on clinical or biological knowledge.

Little work has been done systematically comparing the performance of various functions in the dose-response context. Zhang et al.²⁶ conducted a dose–effect meta-analysis to model the sleep duration and the risk of all-cause mortality, assuming different dose–effect shapes. They found that RCS performed well, while fractional polynomials yielded unreasonable results at 5 and 6 hours of sleep. Additionally, fractional polynomials need intensive computations to find the optimal powers, which is cumbersome to implement in a Bayesian setting. Further work is needed in this direction to study and compare different dose–effect shapes and pinpoint the advantages and limitations of the fractional polynomials in this context.

Our study's model is an extension of our previous work in pairwise meta-analysis.¹¹ Dose–effect pairwise meta-analysis models require transforming the doses into a common scale across agents, which is not always straightforward or even possible. DE-NMA allows us to compare multiple agents simultaneously, using their original doses. It can also answer key questions about what treatments are preferable and what dose can maximise the relative effects. These results from the DE-NMA model are important for drug guideline developers, health technology assessment agencies, and of course patients and their treating clinicians.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802211070256 - Supplemental material for A dose–effect network meta-analysis model with application in antidepressants using restricted cubic splines

Supplemental material, sj-pdf-1-smm-10.1177_09622802211070256 for A dose–effect network meta-analysis model with application in antidepressants using restricted cubic splines by Tasnim Hamza, Toshi A Furukawa, Nicola Orsini, Andrea Cipriani, Cynthia P Iglesias and Georgia Salanti in Statistical Methods in Medical Research

Footnotes

Acknowledgements

This work has been done within the HTx Horizon 2020 project. HTx is supported by the European Union, lasting for 5 years from January 2019. The main aim of HTx is to create a framework for the Next Generation Health Technology Assessment (HTA) to support patient-centered, societally oriented, real-time decision-making on access to and reimbursement for health technologies throughout Europe. The views expressed are those of Prof. Andrea Cipriani and not necessarily those of the UK National Health Service, the NIHR, or the UK Department of Health.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship and/or publication of this article: TH, CI and GS are funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 825162. AC is supported by the National Institute for Health Research (NIHR) Oxford Cognitive Health Clinical Research Facility, by an NIHR Research Professorship (grant RP-2017-08-ST2-006), by the NIHR Oxford and Thames Valley Applied Research Collaboration and by the NIHR Oxford Health Biomedical Research Centre (grant BRC-1215-20005). TAF reports grants and personal fees from Mitsubishi-Tanabe, personal fees from MSD, personal fees from SONY, grants and personal fees from Shionogi, outside the submitted work; In addition, TAF has a patent 2020-548587 concerning smartphone CBT apps pending, and intellectual properties for Kokoro-app licensed to Mitsubishi-Tanabe.

ORCID iDs

Tasnim Hamza

Toshi A Furukawa

Cynthia P Iglesias

Georgia Salanti

Supplemental material

Supplemental material for this article is available online.

Correction (June 2024):

The article has been updated since original publication. For further details, please see

References

Caldwell

Ades

Higgins

JPT

. Simultaneous comparison of multiple treatments: Combining direct and indirect evidence. Br Med J 2005; 331: 897–900.

Ades

. Combination of direct and indirect evidence in mixed treatment comparisons. Stat Med 2004; 23: 3105–3124.

Efthimiou

Debray

TPA

van Valkenhoef

, et al. Get real in network meta-analysis: A review of the methodology. Res Synth Methods 2016; 7: 236–263.

Cipriani

Furukawa

Salanti

, et al. Comparative efficacy and acceptability of 21 antidepressant drugs for the acute treatment of adults with major depressive disorder: A systematic review and network meta-analysis. Lancet 2018; 391: 1357–1366.

Manzoli

Salanti

De Vito

, et al. Immunogenicity and adverse events of avian influenza A H5N1 vaccine in healthy adults: Multiple-treatments meta-analysis. Lancet Infect Dis 2009; 9: 482–492.

Thorlund

Mills

, et al. Comparative efficacy of triptans for the abortive treatment of migraine: A multiple treatment comparison meta-analysis. Cephalalgia 2014; 34: 258–267.

Del Giovane

Vacchi

Mavridis

, et al. Network meta-analysis models to account for variability in treatment definitions: Application to dose effects. Stat Med 2013; 32: 25–39.

Mawdsley

Bennetts

Dias

, et al. Model-based network meta-analysis: A framework for evidence synthesis of clinical trial data. CPT Pharmacometrics Syst Pharmacol 2016; 5: 393–401.

Orsini

Wolk

, et al. Meta-analysis for linear and nonlinear dose-response relations: Examples, an evaluation of approximations, and software. Am J Epidemiol 2012; 175: 66–73.

10.

Crippa

Discacciati

Bottai

, et al. One-stage dose-response meta-analysis for aggregated data. Stat Methods Med Res 2019; 28: 1579–1596.

11.

Hamza

Cipriani

Furukawa

, et al. A Bayesian dose-response meta-analysis model: A simulations study and application. Stat Methods Med Res 2021.

12.

Bagnardi

Zambon

Quatto

, et al. Flexible meta-regression functions for modeling aggregate dose-response data, with an application to alcohol and mortality. Am J Epidemiol 2004; 159: 1077–1086.

13.

Plummer

. JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. 2003.

14.

RStudio Team. RStudio: Integrated development environment for R. RStudio, Inc., 2019. http://www.rstudio.com/.

15.

Hamza

Furukawa

Orsini

, et al. A dose-effect network meta-analysis model: An application in antidepressants. Zenodo 2021. https://doi.org/10.5281/zenodo.4680410.

16.

Cipriani

Furukawa

Salanti

. Comparative efficacy and acceptability of 12 new-generation antidepressants: A multiple-treatments meta-analysis. 2009. https://www.ncbi.nlm.nih.gov/pubmed/19185342.

17.

Furukawa

Cipriani

Cowen

, et al. Optimal dose of selective serotonin reuptake inhibitors, venlafaxine, and mirtazapine in major depression: a systematic review and dose-response meta-analysis. Lancet Psychiatry 2019; 6: 601–609.

18.

Harrell

. Regression modelling strategies: with applications to linear models, logistic regression, and survival analysis. New York: Springer; 2015.

19.

Harrell

. rms: Regression modeling strategies. R package version 5.1-3.1. https://CRAN.R-project.org/package=rms.

20.

Plummer

Best

Cowles

, et al. CODA: Convergence diagnosis and output analysis for MCMC. R News 6: 7–11.

21.

Stone

. Comment: generalized additive models. Stat Sci 1986; 1: 312–314.

22.

Pedder

Dias

Boucher

, et al. Methods to assess evidence consistency in dose-response model based network meta-analysis. Stat Med 2022; 41: 625–44.

23.

Jansen

Vieira

Cope

. Network meta-analysis of longitudinal data using fractional polynomials. Stat Med 2015; 34: 2294–2311.

24.

Royston

Sauerbrei

. Multivariable model-building: A pragmatic approach to regression analysis based on fractional polynomials for modelling continuous variables. New York: John Wiley & Sons, 2008.

25.

Faes

Geys

Aerts

, et al. Use of fractional polynomials for dose-response modelling and quantitative risk assessment in developmental toxicity studies. Stat Modelling 2003; 3: 109–125.

26.

Zhang

Jia

, et al. Introduction to methodology of dose-response meta-analysis for binary outcome: with application on software. J Evidence Based Med 2018: 1–5.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

1.09 MB

$i = 1, \dots, n s$	Study id
$j$	Index for the dose levels in study $i$
$k = 1, \dots, K$	Agent
$c = 1, \dots, C$	Exposure clusters
$p = 1, \dots, P$	Number of dose transformations associated with the dose-effect shape. For a linear shape $p = 1$ and for quadratic and restricted cubic splines $p = 2$
$x_{i j k}$	The jth dose in study i for agent k
$x_{d k}$	Minimal dose d for agent k
$Z_{i}$	Covariate in study i
$r_{i j k}$	Number of events in dose j within study i for agent k
$n_{i j}$	Sample size in dose j within study i for agent k