2.1 Overview of modelling

The DT-MBNMA framework addresses how treatment effects evolve over time and vary across doses by combining three components:

1. 1. Response model (Section 2.2): The model uses aggregate data extracted from publications (means and standard errors [SEs] at each timepoint for each study arm). These summaries are assumed to follow a multivariate normal distribution, with expected values come from parametric functions describing time–course and dose–response (Table 1). SEs are treated as known, extracted from publications while correlations are either set to zero, estimated using simple time series models, or uses observed correlations.

2. 2. Parameter model (Section 2.3): The model estimates structural parameters that govern response functions rather than separate effects for each dose-timepoint combination. Response functions are defined based on subject-matter knowledge. For example, a time–course might be characterized by maximum effect and rate of onset; dose–response by maximum effect and the dose producing half-maximal response. These parameters vary across studies through study-specific reference levels (fixed effects) and treatment-specific deviations (random effects). Dose-–response functions map dose levels to relative effects, enabling information sharing across doses, timepoints, and trials.

3. 3. Random effects specification and diagnostics (Section 2.4): Between-study variability in treatment effects is handled through random effects with a covariance matrix that allows parameters to correlate (e.g., treatments differing in maximum effect may also differ in time-to-effect). Multi-arm trials receive special handling to account for correlation between arms sharing a common reference. This random effects structure provides diagnostics: large between-study variability suggests model misspecification or network inconsistency.

Table 1

Response functions for joint time–course and dose–response modelling

Table 1. Long description

The table has five rows and five columns. The columns are Model, Function f dot, Parameters, Time-course, and Dose-response. Row one is Log-linear with function b times log open parenthesis x plus one close parenthesis, parameter normal prior on b, and x in Dose-response. Row two is Linear with function b times x, parameter normal prior on b, and x in Dose-response. Row three is E max with function E max times x all over exp open parenthesis x sub 50 close parenthesis plus x, parameters normal priors on E max and log x sub 50, and x in both Time-course and Dose-response. Row four is Exponential with function E max times open parenthesis one minus exp open parenthesis minus exp open parenthesis k T close parenthesis times x close parenthesis close parenthesis, parameters E max and log k T, and x in Time-course. Row five is Bateman with function A times open parenthesis exp open parenthesis k T sub one close parenthesis all over exp open parenthesis k T sub one close parenthesis minus exp open parenthesis k T sub two close parenthesis close parenthesis times open parenthesis exp open parenthesis minus exp open parenthesis k T sub two close parenthesis times x minus exp open parenthesis k T sub one close parenthesis times x close parenthesis close parenthesis, parameters A, log k T sub two, and log k T sub one, and x in Time-course. The note states that x denotes function explored for modeling time-course, dose-response, or both.

Note: The mark ‘x’ denotes function that explored for modelling time–course, dose–response or both.

The framework nests simpler approaches: replacing parametric dose–response with discrete treatment effects yields T-MBNMA; using only end-of-treatment timepoints yields standard NMA. This nested conceptualization allows assessment of whether added complexity improves evidence synthesis.

2.2 Response model

We consider the typical setting where only aggregate-level trial data are available from publications or trial reports. Let y _ij denote a vector of summary statistics for study (i = 1 ,…, k) and arm (j = 1 ,…, n_i) across the time–course profile (t = [0 ,…, τ_ij]) with dimension (n_i x (τ_ij+1)). Here n_i denotes the number of arms in study i, and τ_ij denotes the end-of-treatment timepoint in study i and treatment arm j. Our framework applies to outcomes measured repeatedly at protocol-defined timepoints where time–course relationships are biologically plausible. Repeated observations on the individual are correlated over time and therefore it is reasonable to specify a multivariate normal likelihood for the aggregate data:

(1) $$\begin{align}{\mathbf{y}}_{ij}\sim N\left({\boldsymbol{\mu}}_{ij},{\boldsymbol{\Sigma}}_{ij}\;\right),\end{align}$$

where μ _ij is the vector of modelled outcomes generated from the response functions in Table 1, which are motivated in Section 2.2 of the article. The covariance matrix is denoted by ∑ _ij. If the clinical outcome is given in absolute values, then an intercept should be added to the equations in Table 1, addressing the outcome at time = 0.

The diagonal elements in ∑ _ij are the within-arm variances (se ² _ij), and in line with standard NMA, are assumed known. The within-arm SEs can be extracted or digitized from publications, whereas within-arm correlations seldom are reported. Therefore, we follow the approach of Pedder et al.Reference Pedder, Dias, Bennetts, Boucher and Welton ³ by either setting within-arm correlations to zero (i.e., a univariate normal likelihood) or estimating within-arm correlations based on the aggregate response summaries. Based on the work of Pedder et al.,Reference Pedder, Dias, Bennetts, Boucher and Welton ³ an autoregressive of first order structure (AR-1) appears reasonable because covariances are dependent on the time distance between clinical visits. In this article, data are simulated at the individual patient level (IPD). In the illustrative dataset, IPD is available to investigate the consequence of ignoring or approximating the within-arm correlation structure because correlations estimated from IPD will be at the individual level, which may be different to summary-level correlations.

2.3 Parameter model

In statistical terms, a random effects DT-MBNMA is a nonlinear mixed effect (NLME) model. NLME models are widely used in pharmacometrics, particularly for characterizing the time–course of drug concentrations and their effects on the body. The variation within and between studies is modelled through the structural parameters (e.g., time-constant (kT) or maximum effect (E _maxT) for the exponential function) of the response functions in Table 1, which are presented in detail in Section 2.4. To generalize, the model for the p’th structural parameter ( ${\theta}_{ij}^{(p)}$ ) can be written as

(2) $$\begin{align}g\left({\theta}_{ij}^{(p)}\right)=\left\{\begin{array}{cc}{\mu}_i^{(p)},\kern15em & \kern5em j=1\;(\textit{reference}),\\ {}{\mu}_i^{(p)}+f\left({\mathbf{X}}_{\boldsymbol{ij}}{\boldsymbol{\unicode{x3b2}}}^{(p)},{d}_{ij}\right)+{\gamma}_{ij}^{(p)},\kern4.2em & \kern0.24em j>1,\end{array}\right.\end{align}$$

where g(.) is a link function, e.g., the time-constants are log transformed, ensuring that they are positive real numbers on the time scale. The ${\mu}_i^{(p)}$ is the study-specific level for the reference treatment in the ith study, which is treated as a fixed effect and can be viewed as a nuisance parameter. The function f(.) is a parametric dose–response function and commonly used response functions are presented in Section 2.2 of the article. The response function f(.) maps the treatment doses to relative treatment effects and ${\mathbf{X}}_{\boldsymbol{ij}}$ is a design matrix, encoding agent specific effects, and ${\boldsymbol{\unicode{x3b2}}}^{(p)}$ is a vector of agent specific dose–response parameters (e.g., median effective dose (ED50) and maximum effect (E _maxD)). The proposed framework has the flexibility to model time–course profiles that are dose-dependent.

2.4 Random effects specification and diagnostics

The ${\gamma}_{ij}^{(p)}$ are the study-specific treatment effects, which are assumed to be exchangeable and handled as random effects. If ${\gamma}_{ij}^{(p)}$ are set to zero or dropped, then equation (2) reduces to a fixed effects model. Let ${\boldsymbol{\unicode{x3b3}}}_i$ be of vector of study-specific random effects of length n _i–1 x n _p, where n _i is the number of arms of study i as before, and n _p is the number of structural parameters (e.g., 3 for the Bateman equation, Table 1). Then, it is assumed:

(3) $$\begin{align}{\boldsymbol{\unicode{x3b3}}}_i\sim N\left(\boldsymbol{0},\boldsymbol{\Omega} \otimes {\mathbf{T}}_i\right),\end{align}$$

where $\boldsymbol{\Omega}$ is an unstructured covariance matrix (n_p x n_p) for parameters, ${\mathbf{T}}_i$ is the known covariance matrix (n _i–1 by n _i–1) with 1 on the diagonal elements and 1/2 on the off-diagonal elements, accounting for correlated random effects in multi-arm studies.Reference Higgins and Whitehead ⁷ The Kronecker product is denoted by $\otimes$ .

The framework outlined in equation (2) accommodates T-MBNMA,Reference Pedder, Dias, Bennetts, Boucher and Welton ³ where the agent-by-dose combinations are viewed as separate treatments. If the parametric dose–response function f(.) in equation (2) is replaced by a vector of relative treatment effects for the agent-by-dose combinations, then the model corresponds to a T-MBNMA.Reference Pedder, Dias, Bennetts, Boucher and Welton ³ This relaxes the assumptions around the parametric form of the dose–response, i.e., monotonicity, but does assume consistency. The random effects specification provides diagnostics for assessing dose–response model fit and the consistency assumption. Misspecification of the dose–response relationship increases between-study standard deviations (SDs) for the structural parameters defining treatment effects, providing a useful diagnostic for selection of an appropriate dose–response function.Reference Pedder, Dias, Boucher, Bennetts, Mawdsley and Welton ⁸ Between-study heterogeneity, quantified through Ω, reflects variability in structural parameters not explained by the specified dose–response and time–course models.

Consistency can be assessed by comparing model fit and between-study SDs across nested model specifications. Fitting all direct connections in the network extends the framework to the unrelated mean effects model (UME), relaxing both dose–response assumptions and the consistency assumption.Reference Pedder, Dias, Boucher, Bennetts, Mawdsley and Welton ⁸ Substantially better model fit or lower between-study SDs in the UME is suggestive of inconsistency. Deviance–deviance plots comparing the parametric models to the UME provide a graphical diagnostic for detecting inconsistency. Stable between-study SDs when moving from UME to T-MBNMA to DT-MBNMA, with points clustering near the line of equality in deviance–deviance plots, indicates that parametric restrictions are compatible with the data and the consistency assumption is plausibly satisfied.

Residual heterogeneity following dose–response and time–course modelling may indicate effect modification by unmeasured covariates, suggesting potential transitivity violations. Similar to NMA, prior to analysis, the exchangeability of study populations should be assessed in order to reduce potential transitivity violations. Visual predictive checks, a common method used to diagnose and evaluate pharmacometric models,Reference Bergstrand, Hooker, Wallin and Karlsson ⁹ provides a complementary assessment of whether the model captures observed data patterns across all study arms.

2.5 Response function specification and model selection

Parametric response functions are used to link the relationship between clinical outcome to time and dose. The response functions used in the paper are listed in Table 1. Throughout the article ‘T’ is used to denote a time–course parameter and ‘D’ denotes a dose–response parameter. The listed functions (Table 1) include the E _max function, which is widely used in pharmacology and drug development to understand and quantify the relationship between drug dose/concentration and effect, aiding in the determination of optimal dosages and the prediction of drug responses under different doses/concentrations. It has two parameters E _maxD and ED50, which are the maximum effect (as dose →∞) and the median effective dose. The log-linear and linear functions are obvious choices if efficacy endpoints are far from plateauing for incremental doses. The framework is not limited to the response functions considered here.

Time–course relationships are usually characterized by a nadir in treatment response and the time delay between treatment initiation and nadir. This is conveniently expressed using functional forms that are parametrized with these quantities. The E _max function has been applied,Reference Pedder, Dias, Bennetts, Boucher and Welton ³ but also the exponential function captures the relationship, where the time–course parameters are denoted E _maxT (as time →∞) and ET50 (time to 50% maximum response). The time-constant in the exponential function (i.e., kT) can be expressed in terms of a pharmacological half-life (i.e., log(2)/kT) and the effect, when time →∞ as E _maxT, representing the nadir in treatment response. Finally, we introduce the Bateman equation, which is a function of two exponentials with a biphasic pattern and is well known from pharmacokineticsReference Garrett ¹⁰ (i.e., one compartment model with linear absorption and elimination). The parameters can be associated with observable clinical events such as improvement in response and return to baseline phenomena (e.g., due to disease progression or trial fatigue). The time-constants can be expressed in terms of half-lives (i.e., log(2)/kT₁ and log(2)/kT₂) like with the exponential function. For model comparison, the expected log posterior density (ELPD) criterion is used as defined in the loo package in R,Reference Vehtari, Gelman and Gabry ¹¹ where larger criterion values reflect a better-performing model. It approximates the leave-one-out cross-validation and has been shown to have various advantages over simpler estimates of predictive error such as the deviance information criterion (DIC).Reference Vehtari, Gelman and Gabry ¹¹ Moreover, an approximate SE for estimated predictive errors can also be obtained, which is useful for comparison of predictive errors between two or more competing models.

2.6 Simulation study

The following section describes a simulation study undertaken to understand the robustness of MBNMA under different scenarios and to compare it to standard NMA and T-MBNMA in terms of the gain in precision of relative treatment effects. We emulate the situation where phase 2 data are included in the evidence synthesis.

A network of three agents (q = [A, B, C]) is simulated. We assume a typical development program, where one phase 2 dose-finding trial is conducted, having a placebo/standard-of-care arm and four arms with increasing dose levels. This is followed by two phase 3 confirmatory trials with a placebo/standard-of-care arm and two treatment arms. Finally, a phase 3 trial is conducted that provides direct evidence of A. against B. A table with design features is given in Table S1 in the Supplementary Material, and the corresponding network is drawn in Figure S1 in the Supplementary Material. The clinical outcome is modelled from an exponential function (i.e., ${E}_{max}T\left(1-\exp \left(-\exp (kT)\times t\right)\right)$ , Table 1), having two parameters E _maxT and kT. Further, it is assumed that these parameters change in dose-dependent manner, where E _maxT will change according to an E _max function and kT will change according to a log-linear function. Hence, the simulation model for the study-by-arm specific parameters can be written as

(4) $$\begin{align}{log}\left({kT}_{ij}\right)=\left\{\begin{array}{l@{\qquad}l}{\mu}_i, & j=1\;(\textit{reference}),\\ {\mu}_i+{b}^{(q)}\mathit{\log}\left({d}_{ij}+1\right)+{\gamma}_{ij}^1, & j>1,\end{array}\right.\end{align}$$

(5) $$\begin{align}\quad{E}_{max}{T}_{ij}=\left\{\begin{array}{l@{\qquad}l}{\mu}_i, & j=1\;(\textit{reference}),\\ {\mu}_i+\frac{E_{max}{D}^{(q)}\times {d}_{ij}}{\exp \left( ED{50}^{(q)}\right)+{d}_{ij}}+{\gamma}_{ij}^2, & j>1.\end{array}\right.\end{align}$$

The study-specific treatment effects are assumed random and are generated from equation (3). The numerical values used in the simulations are given in Table S2 in the Supplementary Material. In the simulations, the long-term treatment effect is given as E _maxT in the exponential function (i.e., effect as time →∞). Therefore, indirect treatment comparisons (ITCs) are inferred from E _maxT in the T-MBNMA and DT-MBNMA, whereas NMA uses the week 52 study results, as NMA cannot extrapolate. The true mean time–course and dose–response profiles are presented in Figure S2A and S2B in the Supplementary Material. In Figure S2A in the Supplementary Material, the true mean time–course profile for higher doses of agent C, the time–course has not reached its nadir (i.e., stabilized) and hence by design the full effect size cannot be quantified unless the time–course treatment profile is used and modelled. This design feature is introduced to showcase the limitations of standard NMA (i.e., the week 52 timepoint) for handling a clinical endpoint that is evolving over time.

Three scenarios are studied, where treatment duration and study selection are varied. These are:

1. 1. A case, where all studies and treatment arms are observed for 52 weeks. This is considered the base case.

2. 2. A case, where phase 2 for agent C is truncated to 28 weeks. This is considered the realistic case, where the phase 2 study for agent C is conducted with a shorter treatment duration.

3. 3. A case, where phase 2 is truncated to 28 weeks and phase 3 for agent C is excluded. This is a scenario that mimics an ITC feasibility assessment, where the proposed statistical framework is used to predict the long-term efficacy of agent C against agents A and B. Hence, the statistical framework is used as part of dose setting for phase 3 development for agent C, i.e., model-informed drug development.

Two ITCs are assumed of interest, comparing ‘agent A – dose 1 mg’ vs. ‘agent C – dose 2 mg’ and ‘agent B – dose 3 mg’ vs. ‘agent C – dose 2 mg’. For simplicity and ease of notation, ITCs are denoted ‘A vs. C’ and ‘B vs. C’, going forward.

The IPD are simulated at timepoints 2, 4, 8, 12, 16, 20, 28, 36, 44 and 52 weeks after randomization from a multivariate normal distribution with the mean given by the exponential response function and arm-level parameters specified by equations (4) and (5). The individual patient covariance matrix is presented in Table S3 in the Supplementary Material. Similar to published results, each treatment arm IPD within a clinical study is aggregated into a mean time–course profile and corresponding SE. In addition, the within-arm correlation matrix is also derived and used in the parameter estimation. Three analysis models are considered:

1. 1. Ignore the within-arm correlation structure by specifying a univariate normal likelihood. Here, the within-arm correlation structure is set to zero.

2. 2. Account for within-arm correlation structure by specifying a multivariate normal likelihood. Here, the within-arm correlation structure is approximated with an AR-1 process.

3. 3. Account for within-arm correlation structure by specifying a multivariate normal likelihood. Here, the observed within-arm correlation matrix is used in the estimation.

We estimate parameters using the random effects DT-MBNMA with an exponential time–course function and the dose–response specified in equations (4) and (5). As a comparator estimation method, we also apply standard NMA, where absolute aggregated response at week 52 is used as the dependent variable. A T-MBNMA used as an additional comparator model, where absolute aggregated response at all timepoints is used as the dependent variable. For both analysis models, phase 2 is excluded. Simulation results are summarized in standard measures: Median absolute deviation (MAD), root mean square error (RMSE), and bias.

2.7 Illustrative example: A network of studies in obesity

Glucagon-like peptide-1 receptor agonists (GLP-1 RAs), used to treat type 2 diabetes, have been shown to be effective in generating weight loss in preclinical and clinical studies. Thus, the methodology is illustrated using a dataset of RCTs investigating obesity medication or obesity therapy with GLP-1 RAs. Studies used in the current analysis are published in peer-reviewed journals, comparing the effect of a GLP-1 RA to another GLP-1 RA or placebo (both as adjunct to lifestyle interventions) for weight reduction in patients with obesity. Although individual patient data were available, we analysed only aggregate summary statistics (means and SEs at protocol-defined timepoints) to illustrate the methodology as it would be applied to published trial reports where IPD is unavailable. GLP-1 RAs approved for weight management include Saxenda^{® 12} (active agent: liraglutide) and Wegovy^{® 13} (active agent: semaglutide). For Saxenda, the approved dosage is liraglutide 3.0 mg once-daily^{® 12} and for Wegovy it is semaglutide 2.4 mg and 1.7 mg once-weekly^®. ¹³ The primary endpoint in weight management studies is the percent change in body weight (BW) from baseline to end of treatment, which is collected at selected timepoints during treatment. Neither the end-of-treatment timepoint nor the selected timepoints during treatment are aligned across studies (Table S4 in the Supplementary Material). A phase 2 dose-finding studyReference Astrup, Rossner and Van Gaal ⁶ for liraglutide was conducted, investigating the dose–response of incremental liraglutide doses on weight and a large phase 3 trialReference Pi-Sunyer, Astrup and Fujioka ¹⁴ establishing safety and efficacy of liraglutide 3.0 mg once-daily. For semaglutide, a phase 2 dose-finding studyReference O'Neil, Birkenfeld and McGowan ¹⁵ investigating the dose–response of increasing once-daily doses of semaglutide on weight, is available. Semaglutide is a GLP-1 analogue that has an elimination half-life of approximately 1 week,Reference Kapitza, Nosek, Jensen, Hartvig, Jensen and Flint ¹⁶ supporting once-weekly dosing, and therefore a once-weekly dose regimen was selected for the phase 3 development program. Hence, the semaglutide once-weekly dose equivalents are computed by multiplying the once-daily doses in phase 2 with a factor of 7, which is supported by pharmacokinetic data collected in the trials.Reference Strathe, Horn and Larsen ¹⁷ Four phase 3 studiesReference Wilding, Batterham and Calanna ^{18 –} Reference Rubino, Greenway and Khalid ²¹ are included, investigating safety and efficacy of 2.4 mg semaglutide and additionally one of the studies investigated safety and efficacy of 1.7 mg once-weekly.Reference Kadowaki, Isendahl and Khalid ²⁰

The treatment network is shown in Figure 1, where Figure 1a is the full treatment network and Figure 1a is reduced to phase 3 studies.

Figure 1

Network of treatment comparisons from seven RCTs, where each treatment is represented by a node. Where direct RCT evidence exists for a particular comparison, the nodes are connected by a line, the thickness of which is proportional to the number of comparisons. Network ‘a’ is the full network, including phase 2 and ‘b’ includes only phase 3 data.

Figure 1. Long description

Panel A at the top displays a dense network diagram with eleven nodes labeled placebo, sema F 0.40 mg, sema F 0.30 mg, sema 2.4 mg, sema 1.7 mg, sema 0.40 mg, sema 0.30 mg, sema 0.20 mg, sema 0.10 mg, sema 0.05 mg, lira 1.2 mg, lira 1.8 mg, lira 2.4 mg, and lira 3.0 mg. The largest node is placebo, followed by sema 2.4 mg and lira 3.0 mg, indicating higher total sample sizes. Nodes are connected by lines of varying thickness, with the thickest lines between placebo and sema 2.4 mg, and between placebo and lira 3.0 mg, representing the highest number of direct studies. Thinner lines connect other treatment pairs, forming a highly interconnected web. A legend at the right shows line thickness for 1 to 4 studies and node size for total sample sizes of 1000, 2000, and 3000. Panel B at the bottom shows a simplified network with only four nodes: sema 2.4 mg, placebo, sema 1.7 mg, and lira 3.0 mg. The thickest line connects sema 2.4 mg and placebo. Other lines connect each node to the others, forming a square with diagonals. The same legends for number of studies and sample size appear at the right.

The modelling workflow is defined as follows:

1. 1. Time–course models (Table 1) not including a dose–response component are fitted and ranked according to the goodness of fit using the ELPD criterion. First, common-effects models are fitted, followed by the introduction of exchangeable treatment effects, and finally an AR-1 process is introduced for the arm-correlation structure. The dose–response is specified as a set of discrete treatment effects (i.e., a treatment-level MBNMA), making no assumption about the shape of the dose–response, but does assume consistency.

2. 2. Based on the best performing time–course model, parametric dose–response models are introduced into the structural parameter of the time–course functions. These models are compared to the treatment-level MBNMA and the UME, using graphical methodsReference Pedder, Dias, Boucher, Bennetts, Mawdsley and Welton ⁸ and information criteria. For the selected dose–response model in this step, the observed within-arm correlations are introduced and used in the estimation of dose–response parameters.

Two standard NMAs are conducted for comparison purposes. The first uses primary statistical analysis (including imputation of missing data) that was defined in the individual study protocols and evaluated at the per protocol specified end-of-treatment timepoint (spanning from week 56 to week 104). The second uses the aggregated data (as observed) at the per protocol specified end-of-treatment timepoint, which is introduced to align with aggregated time–course data used in the DT-MBNMAs.

2.8 Implementation

Estimation and inference is carried out using Stan (Stan Development Team, version 2.26.1), where the R package brms provides an interface to R (version 2.21.1).Reference Brkner ²² All computations are conducted with R version 4.3.1 (a GitHub repo is available at https://github.com/agorstras/dt-mbnma). Models that included the observed within-arm covariance matrix are implemented directly in Stan because these cannot be fitted with the current version of brms. Choice of prior distributions and convergence diagnostic are presented below and in-detail in the Supplementary Material.

2.9 Convergence monitoring and prior specification

In stochastic simulation and estimation part, two independent chains were run for 1,000 burn-in iterations followed by 1,000 iterations. Gelman’s R statistic was monitored for the model parameters during the simulations, and it appeared that models converged as all parameters had an R < 1.05. In the analysis of the illustrative dataset, two independent chains were run for 2,500 burn-in iterations followed by 2,500 iterations. Gelman’s R statistic was monitored for the model parameters during the simulations, and it appeared that models converged as all parameters had an R < 1.05. If the models did not converge, the number of iterations were tripled.

Vague normal prior distributions with a SD of 100 were given to the nuisance and drug effect parameters (e.g., b or x50, Table 1). For parameters that only take positive values (e.g., time-constants, median effective dose), priors were expressed on a log-scale. The between-study SDs on exchangeable parameters were given half-normal prior distributions with an SD of 5. The prior for correlation matrix for the exchangeable parameters uses the Lewandowski-Kurowicka-Joe prior (LKJ) is used with a degree freedom of 1. In models with a multivariate likelihood, the correlation parameter was given a uniform prior distribution (uniform(–1,1)) to reflect that outcomes at different timepoints within the same study arm are correlated.