Characterising the dose-response relationship and finding the right dose are important but challenging in the pharmaceutical drug development process. Nearly half of failures in Phase III trials result in part from a lack of understanding of the dose-response relationship in Phase II trials (Sacks et al., 2014). Over the last decade, the Multiple Comparison Procedure-Modelling (MCP-Mod) method, developed by Bretz et al. (2005), has been increasingly popular for Phase II trials as it can provide superior statistical evidence for dose selection.

MCP-Mod is a two-step approach that combines MCP principles and modelling techniques. In the first MCP step, it establishes a dose-response signal (Proof of Concept, PoC), and in the second Mod step, it estimates the dose-response curve and target doses of interest. MCP-Mod overcomes shortcomings of traditional approaches for dose-finding studies (see e.g. Ting, 2006, for an excellent introduction). MCP-Mod requires the pre-specification of plausible candidate models and model parameters to capture model uncertainty, which however is primarily based on the limited knowledge of the agent/compound being studied, if it is available at all (Chen and Liu, 2020). Misspecification of candidate models and model parameters in MCP-Mod may cause a loss in power and unreliable model selection (see Saha and Brannath, 2019, and references therein).

Motivated by the limitations of MCP-Mod, non-parametric methods have gained popularity for detecting dose-response trends and estimating dose-response relationships in Phase II trials. These methods offer flexibility and adaptability to capture complex patterns in the data without imposing strong assumptions on the underlying functional form of the dose-response curve. West and Harrison (2006) introduced the normal dynamic linear model (NDLM), which leverages the flexibility of dynamic linear models to estimate the dose-response curve. By accommodating time-varying coefficients, NDLM can capture the dynamic nature of the dose-response relationship, allowing for more accurate and interpretable estimates. Kirby et al. (2009) applied cubic smoothing splines with generalised cross-validation for the smoothing parameter to estimate the dose-response curve.

However, these approaches still face challenges in incorporating prior knowledge about the curvature or shape of the dose-response curve. Determining the optimal level of smoothness or choosing the appropriate model can be subjective and dependent on the specific data at hand. To address these limitations, in this work, we develop a model-free Bayesian approach, which is a novel Bayesian hierarchical framework incorporating the total (in the L 2 sense) curvature of the dose-response curve as a prior parameter. Our approach avoids the requirement of a set of pre-specified candidate models. The responses at the given set of doses are estimated through maximum a posteriori (MAP), with which we construct a test statistic to establish PoC through simulations. We can then estimate the dose-response relationship with simple interpolation.

The remainder of this work is organised as follows. We describe in detail our MAP approach with a curvature prior, abbreviated LiMAP-curvature, in Section 2. In Section 3, we evaluate the operating characteristics of LiMAP-curvature through simulations and compare its performance to that of MCP-Mod and smoothing spline. We present concluding remarks and future directions in Section 4.

We consider a trial with a total of M + 1 distinct doses x 0 , x 1 , … , x M , where x 0 represents placebo. We let N i be the number of patients in dose group i and f ( x ) be the true dose-response function at dose x. We assume that f ( x ) is defined on the interval [ 0 , 1 ] to align with the typical range of doses in pharmaceutical drug development where doses are often scaled or normalised within this range for convenience and comparability. For i = 0 , 1 , … , M and j = 1 , 2 , … , N i , we let Y i j be the response observed for patient j allocated to dose x i . We assume that (2.1) Y i j = μ i + ϵ i j , where μ i = f ( x i ) denotes the mean response at dose x i , and ϵ i j ∼ i i d N ( 0 , σ 2 ) denotes the error term for patient j in dose group i. For the sake of simplicity, we assume that the variance σ 2 is known and constant across all dose groups. This assumption is commonly used in different areas, including MCP-Mod (Bretz et al., 2005), BMCPMod (Fleischer et al., 2022) and Bayesian meta-analysis (Burke et al., 2018). See Fleischer et al. (2022) for further discussion on the assumption on the standard deviation σ.

In MCP-Mod, a set of plausible candidate models is required. This constrains the possible set of dose-response curves. However, with limited knowledge of the agent/compounds in Phase II trials, it is more likely to mis-specify the set of candidate models. The procedure for specifying candidate models is also somewhat cumbersome. So we would like to avoid the pre-specification of possible models beforehand, but still want to impose a certain degree of smoothness on the dose-response curve. To this end, we introduce the L 2 -total curvature S f = ( ∫ x 0 x M f ″ ( x ) 2 d x ) 1 / 2 to measure how far the dose-response curve f ( x ) is from being a straight line. The L 2 -total curvature is chosen as the prior parameter because it provides a mathematically tractable and interpretable measure of the smoothness of the dose-response curve. Specifically, it quantifies the global curvature by integrating the squared second derivative of the function, ensuring that the method captures overall trends in the data rather than being overly influenced by local deviations. We will impose a half-normal H N ( γ 2 ) prior on S f to give low prior probabilities to the dose-response relationship that is very curved, where the standard deviation γ controls the trade-off between the L 2 -total curvature of f ( x ) and fidelity to data Y. To capture an appropriate level of curvature and allow for optimal model performance, we assign a hyperprior for γ with a half-normal distribution H N ( τ 2 ), which reflects our prior beliefs about the likely range of values for γ before observing any data. Smaller values of τ favour smoother curves, which may be appropriate when the dose-response relationship is expected to be close to linear. Larger values of τ allow for more curvature, which may be necessary when the dose-response relationship is highly non-linear. By incorporating a hyperprior, we introduce additional uncertainty and flexibility into the prior distribution of S f, influencing the estimation of the dose-response curve. The specification of the standard deviation τ will be discussed in Section 4.

Given the dose-response function f ( x ) being available at M + 1 distinct doses, for i = 1 , 2 , … , M − 1, we have f ″ ( x i ) ≈ 2 ( μ i + 1 − μ i ( x i + 1 − x i ) ( x i + 1 − x i − 1 ) − μ i − μ i − 1 ( x i − x i − 1 ) ( x i + 1 − x i − 1 ) ) through the second-order central difference scheme (Burden et al., 2015), therefore the L 2 -total curvature S f being approximated through numerical integration with S μ = 2 ( ∑ i = 1 M − 1 ( μ i + 1 − μ i ( x i + 1 − x i ) ( x i + 1 − x i − 1 ) − μ i − μ i − 1 ( x i − x i − 1 ) ( x i + 1 − x i − 1 ) ) 2 Δ x i ) 1 / 2 , where Δ x i = ( x i + 1 − x i − 1 ) / 2 for i = 2 , 3 , … , M − 2 with Δ x 1 = ( x 2 + x 1 ) / 2 − x 0 and Δ x M − 1 = x M − ( x M − 1 + x M − 2 ) / 2.

We let Y = { Y i j , i = 0 , 1 , … , M , j = 1 , 2 , … , N i } and define our Bayesian hierarchical model to be (2.2) p ( μ , γ ∣ Y ) ∝ p ( μ , γ ) p ( Y ∣ μ ) , where the prior p ( μ , γ ) = p ( γ ) ∏ i = 0 M p ( μ i ) p ( S μ ∣ γ ) with S μ ∣ γ ∼ H N ( γ 2 ) , γ ∼ H N ( τ 2 ) , μ i ∼ i i d U ( 0 , 1 ) for i = 0 , 1 , … , M , and the likelihood p ( Y ∣ μ ) = ∏ i = 0 M ∏ j = 1 N i p ( Y i j ∣ μ i ) .

Suppose the standard deviation τ in the hyperprior for γ is pre-specified. The MAP estimates of all parameters in the model defined in Eq. (2.2) can be obtained by maximising the corresponding log-likelihood μ ˆ , γ ˆ = arg max μ , γ log p ( μ , γ ∣ Y ) ∝ arg max μ , γ { − ∑ i = 0 M ∑ j = 1 N i ( Y i j − μ i σ ) 2 − 2 log γ − ( S μ γ ) 2 − ( γ τ ) 2 } through a numerical optimisation algorithm like the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method and its variants (see, e.g. Nocedal and Wright, 1999, for more details).

To establish PoC, we propose a test statistic T = max { μ 1 − μ 0 , μ 2 − μ 0 , … , μ M − μ 0 } with hypotheses H 0 : μ 1 = ⋯ = μ M = μ 0 H 1 : max { μ 1 , μ 2 , … , μ M } > μ 0 . We define a significance level α for a dose-response signal, such that the corresponding critical value c satisfies P ( T > c ∣ H 0 ) < α . For a given α, we compute the critical value c via simulation P ( T > c ∣ H 0 ) ≈ 1 R ∑ r = 1 R 1 { T ( r ) > c } , where T ( r ) is the test statistic computed from the data Y ( r ) simulated under the null hypothesis H 0 , R is the total number of replicates, and 1 A is the indicator function equal to 1 if condition A holds and 0 otherwise.

Finally, the mean response estimates μ ˆ can be linearly interpolated (or using a more sophisticated interpolation scheme) to obtain an estimate of the dose-response curve f ( x ), which will also yield an estimate of the target dose of interest. In this work, we only perform simple linear interpolation as we shall focus on the PoC stage of dose-finding trials.

In this section, we assess the performance of LiMAP-curvature, in terms of power to detect dose-response signals as well as estimates of the dose-response curve and target doses of interest. We compare LiMAP-curvature to smoothing spline and MCP-Mod.

In this work, we have introduced LiMAP-curvature, a novel Bayesian approach for establishing PoC and estimating the dose-response curve alongside target doses of interest in Phase II trials. LiMAP-curvature is “model-free”, in the sense that it does not require pre-specification of candidate dose-response models, which can influence the performance of MCP-Mod. It is built on a Bayesian hierarchical framework incorporating prior information on the L 2 -total curvature of the dose-response curve.

We have shown through simulations that LiMAP-curvature has performance comparable to that of MCP-Mod in establishing PoC and estimating MED when the true dose-response model is included in the candidate model set of MCP-Mod, which is fairly rare in practice. When the true dose-response model deviates from the candidate model set of MCP-Mod, LiMAP-curvature has been demonstrated to outperform MCP-Mod. Furthermore, we note that LiMAP-curvature also outperforms smoothing spline in terms of power to detect dose-response signals in most cases. This indicates the advantages of LiMAP-curvature over the widely used non-parametric method of smoothing spline. The ability of LiMAP-curvature to consistently outperform both MCP-Mod and smoothing spline across various scenarios highlights its effectiveness and potential as a valuable tool in Phase II trial analysis.

To obtain optimal performance, MCP-Mod requires specialised expertise to pre-specify plausible candidate models and model parameters, but the knowledge of the agent/compounds being studied is usually limited. Compared to MCP-Mod, the only requirement for pre-specification in LiMAP-curvature is the standard deviation τ for the hyperprior γ ∼ H N ( τ 2 ), which encodes prior knowledge of how far the dose-response curve is from a straight line. With additional simulations with varying values of the standard deviation τ (see Supplemental Material, Figure S12 and Table S5), we recommend choosing

A number of relevant issues for LiMAP-curvature deserve further research. One such extension involves integrating the sigmoid Emax model as the default response function within the LiMAP-curvature framework. The sigmoid Emax model is noted for its ability to approximate most common monotonic dose-response relationships, providing a robust framework for modelling complex dose-response relationships. However, the challenge lies in selecting appropriate prior distributions for the additional parameters introduced by the sigmoid Emax model, which requires careful consideration to ensure accurate and reliable dose-response estimation. Another important extension involves adapting LiMAP-curvature to handle a wider variety of endpoints and trial designs. The current version of LiMAP-curvature is limited to the analysis Phase II dose-finding trials with continuous endpoints, i.e. a single normally distributed homoscedastic response measured at the end of the trial. To expand its applicability, extensions to other common types of endpoints (e.g. binary, counts and survival endpoints) and trials (e.g. longitudinal dose-finding trials) require further investigation. Other directions of future research include the investigation of trial designs tailored to LiMAP-curvature and the development of statistical software implementing LiMAP-curvature.

The authors report there are no competing interests to declare.

The opinions expressed in this paper are solely those of the authors and not those of their affiliations. The authors’ affiliations do not guarantee the accuracy or reliability of the information provided herein.