## Accounting for parameter uncertainty in two-stage designs for Phase II dose-response studies

## Emma McCallum 1 and Bj¨rn Bornkamp o 2

1 MRC Biostatistics Unit, Cambridge, United Kingdom

2 Novartis Pharma AG, Basel, Switzerland

## Abstract

In this paper we consider two-stage adaptive dose-response study designs, where the study design is changed at an interim analysis based on the information collected so far. In a simulation study, two approaches will be compared for these type of designs; (i) updating the study design by calculating the maximum likelihood estimate for the dose-response model parameters and then calculating the design for the second stage that is locally optimal for this estimate, and (ii) using the complete posterior distribution of the model parameter at interim to calculate a Bayesian optimal design (i.e. taking into account parameter uncertainty). In particular, for an early interim analysis respecting parameter uncertainty seems more adequate, on the other hand for a Bayesian approach dependency on the prior is expected and an adequately thought-through prior is required. A computationally efficient method is proposed for calculating the Bayesian design at interim based on approximating the full posterior sample using k-means clustering. The sigmoid Emax dose-response model and the D-optimality criterion will be used in this paper.

Keywords: adaptive dose-finding, Bayesian optimal design, D-optimality, k-means clustering, Hill equation, sigmoid Emax

## 1 Introduction

In the past few years there has been substantial interest in adaptive dose-finding studies, see for example, the white papers of the PhRMA Working Group for Adaptive Dose-Ranging Designs (Bornkamp et al., 2007; Dragalin et al., 2010). Other references include Grieve and Krams (2005); Dragalin et al. (2007); Miller et al. (2007); Berry et al. (2010); Bornkamp et al. (2011); Jones et al. (2011) among many others. As fixed designs are more often used in clinical practice, one obvious question is to characterize situations when an adaptive design will outperform a non-adaptive design from a statistical efficiency perspective. In practice, of course also non-statistical (or alternative statistical) considerations will play a role at the trial design stage, when one decides for or against using an adaptive design, but we will focus on statistical estimation efficiency in this paper.

To that end, Dette et al. (2013) investigated two-stage designs, where an interim analysis is performed and the maximum likelihood (ML) estimate of the model parameters is used to calculate a locally optimal design for the second stage of the trial. Based on analytical approximations they identified key factors, when adaptive designs are beneficial, including signal to noise ratio (the larger the better for adaptive designs), timing of the interim analysis (if the interim estimate is too early the interim estimate might be too noisy to produce a good design, whereas if the interim analysis is too late there might be not enough patients left to allocate into the trial for the calculated optimal design), and adequacy of the starting design (if an adequate starting design is used, there might not be much to improve).

A problem of the ML approach discussed in Dette et al. (2013) is that the study design is optimized based on the ML point estimate at the interim. By using only a point estimate, the uncertainty in the parameter estimate is ignored and

potentially plausible parameter values discarded, ultimately leading to designs optimized for the 'wrong' parameter values. Uncertainty can be accounted for at the interim analysis by adopting a Bayesian approach with prior distributions for the model parameters specified at the beginning of the trial and the posterior distribution based on the first stage data being utilized at the interim analysis to update the study design.

The primary aim of this work is to investigate whether there are any gains to be made by accounting for parameter uncertainty at the interim by using Bayesian methods or whether ML updating of the optimal design is preferable. As analytical considerations such as those discussed in Dette et al. (2013) are infeasible for the Bayesian designs considered here, we conducted a large-scale simulation study in the realistically relevant settings of the PhRMA working group. As we concentrate on parameter uncertainty in this paper and not model uncertainty, we only use the sigmoid Emax model (MacDougall, 2006): It is commonly used for modeling the dose-response relationship and its main assumption is monotonicity which is a desirable feature for modeling many clinical trial dose-response curves (Kirby et al., 2011). In addition the sigmoid Emax model has been found to give a good approximation to a number of other possible dose-response models as well (Thomas, 2006). The design criterion used will be D-optimality. The type of adaptation under consideration involves changing the randomisation ratios to existing dose-levels based on efficacy.

## 2 Methodology

## 2.1 Statistical model

Let Y be a clinical outcome observed at a dose level x ∈ [0 , D ] where 0 and D are the placebo and maximum dose, respectively. A parallel group design is

used where several active dose levels plus a placebo are considered, with patients being randomized to one of the treatments. There are q discrete doses available, with q being a fixed number and the grid of available doses being denoted x = ( x , . . . , x 1 q ).

Consider the following regression model

$$Y _ { i j } = \eta ( x _ { i }, \boldsymbol \theta ) + \epsilon _ { i j }, \quad \epsilon _ { i j } \sim N ( 0, \sigma ^ { 2 } )$$

where Y ij denotes the response of patient j at dose x i ( i = 1 , . . . , q , j = 1 , . . . , n i ) and glyph[epsilon1] ij is the residual error, which is assumed to be independent and normally distributed with common variance σ 2 . The true, but unknown, r -parameter doseresponse model is denoted by η x , ( i θ ) with parameter vector θ = ( θ , . . . , θ 1 r ). The total sample size is ∑ q i =1 n i = N where n i denotes the number of patients assigned to dose x i .

## 2.1.1 Sigmoid Emax model

The four parameter sigmoid Emax model is defined as

$$\eta ( x _ { i }, \boldsymbol \theta ) = \theta _ { 1 } + \theta _ { 2 } \frac { x _ { i } ^ { \theta _ { 4 } } } { \theta _ { 3 } ^ { \theta _ { 4 } } + x _ { i } ^ { \theta _ { 4 } } }$$

where θ 1 denotes the placebo effect, θ 2 is the asymptotic maximum treatment effect (often called Emax), θ 3 is the dose that gives half of the maximum treatment effect (often called ED 50 ), and θ 4 is the so-called Hill parameter that determines the steepness of the curve.

The gradient of the sigmoid Emax model, g ( x, θ ), is given as

$$g ( x, \theta ) = \left ( 1, \frac { 1 } { 1 + ( \theta _ { 3 } / x ) ^ { \theta _ { 4 } } }, \frac { - \theta _ { 4 } \theta _ { 2 } \theta _ { 3 } ^ { ( \theta _ { 4 } - 1 ) } } { x ^ { \theta _ { 4 } } ( 1 + ( \theta _ { 3 } / x ) ^ { \theta _ { 4 } } ) ^ { 2 } }, \frac { - \theta _ { 2 } \log ( \theta _ { 3 } / x ) ( \theta _ { 3 } / x ) ^ { \theta _ { 4 } } } { ( 1 + ( \theta _ { 3 } / x ) ^ { \theta _ { 4 } } ) ^ { 2 } } \right ) \quad ( 3 )$$

Here there are discrete dose levels and so the Fisher information matrix for the statistical model given in (1) for a sigmoid Emax model is defined as

$$M ( x, \mathbf w, \theta ) = \sum _ { i = 1 } ^ { q } w _ { i } g ( x _ { i }, \theta ) ^ { T } \mathbf g ( x _ { i }, \theta )$$

where w i = n /N i is the proportion of patients allocated to dose x i and w = ( w , . . . , w 1 q ). If a dose from the grid of available doses x , . . . , x 1 q is not allocated any patients then w i = 0. Note that the Fisher information matrix for the nonlinear sigmoid Emax model (Equation 4) is dependent on the unknown true parameter values θ and the experimental design ξ = ( x w , ).

## 2.2 Two-stage D-optimal design

A fixed (or non-adaptive) design is a design in which observations are taken at pre-specified doses. In a two-stage adaptive design, the data accrued in the first stage is used to determine the design for the second stage, as follows;

- 1. Stage 1: N (1) observations at starting design ξ (1) .
- 2. Interim updating: fit the sigmoid Emax model to the first stage data and use the obtained information about the model parameters θ to calculate the optimal design for the second stage, ξ (2) .
- 3. Stage 2: N (2) = N -N (1) observations based on updated design ξ (2) .

Step 2 in the procedure above will be different for the ML-based updating and the Bayesian updating. In this paper, the entire dose-response curve is of interest

and so the D -optimal design is calculated at the interim analysis. A design is D -optimal if it maximizes the determinant of the Fisher information matrix or, equivalently, minimizes the generalized variance of the parameter estimates. Such an optimal design minimizes the volume of the confidence ellipsoid of the parameters, and minimizes the maximum predicted variance around the doseresponse curve (Atkinson et al., 2007, Ch.10,11).

## 2.3 Interim updating

For both types of updating the same grid of pre-specified discrete doses are available at the interim and so the weights are optimised whilst the doses remain fixed. A continuous optimal design is calculated and then rounded to a fixed number of patients so that the total number of patients N is achieved (Pukelsheim, 2006, Chapter 12).

## 2.3.1 Maximum Likelihood updating

At the interim, the sigmoid Emax model is fitted to the data from the first stage using the "fitMod" function from the DoseFinding R package (Bornkamp et al., 2013) which implements nonlinear least squares regression. Bounds will be used for the nonlinear parameters θ 3 and θ 4 ; this is a sufficient condition for existence of the ML estimate, see Jennrich (1969). The bounds [0 001 1 5 . , . D ] for the θ 3 parameter and [0 5 10] for the . , θ 4 will be used (here D is the maximum used dose). Both boundaries are chosen relatively wide to allow essentially all shapes of the underlying sigmoid Emax shape in the considered dose-range [0 , D ].

After the ML estimate has been calculated the estimates are then input to the Fisher information matrix (Equation 4). The optimal design for the second stage, ξ (2) , is found by maximising the vector of weights, w , in the following expression:

$$\xi ^ { ( 2 ) } = \arg \max _ { w } | M ( x, w, \hat { \theta } ^ { ( 1 ) } ) | = \arg \max _ { w } \det \left [ \sum _ { i = 1 } ^ { q } w _ { i } g ( x _ { i }, \hat { \theta } ^ { ( 1 ) } ) ^ { T } g ( x _ { i }, \hat { \theta } ^ { ( 1 ) } ) \right ] \ ( 5 )$$

This optimisation is carried out using the "optDesign" function in the DoseFinding package in R , which performs a nonlinear optimization using an augmented Lagrange method. The remaining N (2) patients are then assigned to the optimal design, ξ (2) , for the second stage.

## 2.3.2 Bayesian updating

The conceptual framework behind the Bayesian optimal design is the assumption that the information of the parameters of interest and their uncertainty can be adequately captured in the probability distribution. This probability density averages out the parameter values and the criterion is no longer dependent on the parameters. The Bayesian D-optimality criterion is defined in Atkinson et al. (2007) as

$$\Phi ^ { \dagger } _ { B } ( \xi ) = E _ { \theta } \left [ \log | \mathbf M ( \xi, \theta ) | \right ] = \int _ { \theta } \log | \mathbf M ( \xi, \theta ) | p ( \theta ) d \theta$$

where M ( ξ, θ ) is the Fisher information matrix and p ( θ ) is the density capturing the information on θ . This type of design is more robust than a locally optimal design in the sense that the performance will be adequate for parameter values that have relevant probability mass, whereas for a locally optimal design one only considers a single parameter value. In many situations p ( θ ) might be a prior distribution; in this paper the distribution used in the design criterion will be the posterior distribution p ( θ y | (1) ) based on the data observed in the first stage of the trial. Note that this distribution will contain less information ( e.g. it will have a larger posterior variance) if an early interim analysis is performed. When the posterior is calculated on a late interim the density will be more peaked. The

'robustness' of the design criterion hence automatically adapts to the amount of information available. The ML based approach in contrast always uses the point estimate without acknowledging the uncertainty associated with it.

For the sample size and signal-to-noise ratio of interest in dose-finding studies there will be a dependence on the prior distribution specified at the beginning of the trial (Temple (2012); Bornkamp (2012)), in particular if the interim analysis is performed early in the trial. For the purpose of generic simulations, outside a concrete real example, it is difficult to come up with examples of prior information one might have, so the focus in this paper will be on the situation of little prior information. It is clear, that a fully Bayesian approach will be better than the maximum-likelihood approach if relevant historical prior information is available, so the situation considered in the simulations might be considered conservative for the performance of the Bayesian approach.

## The functional uniform prior distribution

One potential weakly informative prior distribution is the Jeffreys prior distribution, as for example discussed in (Gelman et al., 2003, p.63). Here one uses

$$p ( \theta ) \in | \mathsf M ( \mathsf x, \mathsf w, \theta ) | ^ { 1 / 2 },$$

where M ( x w , , θ ) is the Fisher information matrix defined in Equation (4). So the Jeffreys prior depends on the Fisher information matrix, which in turn depends on the observed design ( x w , ). However, in a two-stage trial the design at the interim will be different to the design at study end, so that the Jeffreys prior will be different at interim and study end, which is why the prior violates the likelihood principle (Ghosh et al., 2006).

Functional uniform priors (Bornkamp, 2012, 2014) strive to derive a distribution that is uniformly distributed on the potential different shapes of the under-

lying nonlinear model function. These priors are invariant with respect to parameterization of the model function and typically result in rather non-uniform prior distributions on the parameter scale. When using the L 2 distance at the grid of available doses x , . . . , x 1 q with equal weight to define functional uniformity, one obtains a functional uniform distribution that is essentially a modified version of Jeffreys prior (Bornkamp, 2012). Instead of the actual observed design one uses the x values and weights for w used in the L 2 distance. Here we will use the grid x of all doses available and equal weights w . This prior will then no longer violate the likelihood principle (it is completely specified before trial start) and has the clear interpretation of being uniformly distributed on the different shapes available in the chosen distance metric. The same parameter bounds as for ML estimation will be used for the Bayesian approach. This is the prior used in this paper.

## Numerical calculations to update Bayesian designs

Two numerically challenging tasks are involved in calculation of the updated design at the interim analysis based on the Bayesian approach; (i) one needs to perform the integration in Equation (6) to evaluate the design criterion once, and (ii) one needs to optimize the design criterion with respect to the design, requiring multiple evaluations of Equation (6) for different candidate designs.

One approach to perform the integration is to use Monte Carlo. That means at the interim analysis JAGS will be used to create a Markov chain Monte Carlo sample of size 10000 from the posterior distribution. The integral can then be estimated as the corresponding arithmetic mean. So for 10000 samples from the posterior distribution this requires calculation of 10000 evaluations of the determinant of the Fisher information to evaluate the criterion for one design.

This can become time consuming, in particular because evaluation of the design criterion is embedded in the numerical optimization. In context of a simulation study this approach would be computationally prohibitive.

A more efficient alternative to pure Monte Carlo is to use k-means clustering to approximate the full posterior distribution for the parameters by a weighted discrete distribution. The k-means cluster algorithm partitions n observations into k clusters so as to minimise the within-cluster sum of squares (Wu, 2012). So instead of using all posterior samples, each weighted equally, a weighted mean is used, based on the identified cluster centers (and corresponding cluster weights), so that the determinant of the Fisher information needs to be evaluated less often for one evaluation of the design criterion. The Hartigan and Wong (1979) k-means clustering algorithm implemented in R in the function kmeans will be used.

Theoretically for a large number of cluster centers (and a large posterior sample), the k-means approximation of the integral should be close to the correct integral, in Appendix A we outline some heuristic mathematical arguments that support this. In practice it is unclear how many cluster centers will provide an adequate approximation of the integral and so a small simulation study was conducted, where k = 10 cluster centers were used. Actual calculation of the optimal design was performed with the "optDesign" function, where the cluster centers were used for the model parameters and the corresponding cluster weights specified. The cluster weights correspond to the proportion of posterior samples associated with each of the cluster centres.

We used one of the simulation scenarios studied in the simulation study described later to evaluate the loss of efficiency when using the k-means approximation instead of the full posterior, see Appendix B for the exact simulation scenario. The relative performance is defined as

$$\text{RelativePerformance} = \frac { \Psi _ { B } ( \xi _ { \text{km} } ^ { * } ) } { \Psi _ { B } ( \xi _ { \text{fp} } ^ { * } ) }$$

where ξ ∗ km is the optimal design calculated using the k-means approximation, ξ ∗ fp is the optimal design calculated using the full posterior sample and Ψ B the Bayes design criterion from Equation (6).

Figure 1 shows one example of the distribution of efficiencies for 1000 datasets simulated using a sigmoid Emax model with 15 observations. The sample size was chosen to be small, as in these situations it is assumed that the posterior has the most complex shape, and is most difficult to approximate.

Figure 1: Relative performance of k-means approximation to full posterior

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-666119f5-9d45-4027-a916-496eb5c1b555.jpg)

The k-means approximation for the full posterior is fairly good with an average efficiency of 0 915 for this scenario. . The approximation improves as the sample size in the first stage of the trial increases with an average efficiency of 0 961 and . 0 976, for sample sizes 60 and 150, respectively. .

In a scenario with a total sample size of 250 the k-means approximation required on average 7.7 seconds to calculate the optimal design, compared to 124.9

seconds for the calculation of optimal design using the full posterior distribution (refer to Appendix B for the simulation scenario). The distributions of the running time for 100 simulations are shown in Figure 2, when using (a) the k-means approximation and (b) the full posterior distribution. So a more than 10-fold reduction in computing time can be achieved by using this approximation.

Figure 2: Distribution of running times for calculation of the second stage design

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-5e31841e-a34c-4735-91d9-7f20b7175bbc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-d85d0d60-f612-4851-b362-3c42cddcbf95.jpg)

## 3 Simulations

In total, 448 different scenarios were considered in the simulations. Each scenario correspond to a different combination of starting design, dose-response profile, total sample size, timing of the interim analysis and type of interim updating. The scenarios selected are similar to those conducted by the PhRMA Adaptive Dose-Ranging Designs Working Group (Bornkamp et al., 2007). In the PhRMA Working Group paper, a comprehensive simulation study based on a hypothetical neuropathic pain study was conducted to evaluate existing adaptive dose-ranging methods. The primary endpoint considered was the change from baseline at 6 weeks in a visual analog scale (VAS) of pain and this is taken to be the primary endpoint in the following simulation study. The VAS measurements are on a

continuous scale from 0 (no pain) to 10 (highest pain). If the VAS measurement at the k th week is denoted by V AS k with k = 0 1 , , . . . , 6 and k = 0 representing the baseline, then the primary endpoint is defined as

$$y = V A S _ { 6 } - V A S _ { 0 }$$

Hence, negative values of y indicate efficacy as there is a reduction in neuropathic pain.

## 3.1 Starting designs

In order to investigate the impact of number and spacing of doses, four different starting designs are considered in the simulations (Table 1). For the first stage, patients are allocated equally between the doses considered.

Table 1: Starting designs

| Design A   | Five equally spaced doses                       | ξ A = (0 , 2 , 4 , 6 , 8)                 |
|------------|-------------------------------------------------|-------------------------------------------|
| Design B   | Five unequally spaced doses (mostly low doses)  | ξ B = (0 , 1 , 2 , 4 , 8)                 |
| Design C   | Five unequally spaced doses (mostly high doses) | ξ C = (0 , 6 , 7 , 7 . 5 , 8)             |
| Design D   | Nine equally spaced doses                       | ξ D = (0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8) |

Designs A and B reflect what is most commonly seen in practice, in particular, design B is an equal distribution on the log 2 scale. Design C is more unusual but can occur if the dose range considered is too narrow and the lowest dose is not low enough. Design D is less likely to be seen in practice as it is unusual to see so many dose levels utilized.

At the interim analysis it is assumed that 17 equally spaced doses can be utilized ranging from 0 to 8 in 0 5 increments: . x = (0 0 5 1 1 5 2 , . , , . , , . . . , 7 5 . , 8).

## 3.2 Dose-response profiles

A total of four different dose-response profiles were used to simulate the primary endpoint, which are considered to span a range of dose-response profiles often observed in practice (Figure 3). As in Bornkamp et al. (2007), for all models the placebo effect was set to 0 and the maximum effect within the observed dose range was set to -1 65 units. .

$$\text{Linear} \colon \ \eta ( x, \boldsymbol \theta ) & = - ( 1. 6 5 / 8 ) x \\ \text{Quadratic} \colon \ \eta ( x, \boldsymbol \theta ) & = - ( 1. 6 5 / 3 ) x + ( 1. 6 5 / 3 6 ) x ^ { 2 } \\ \text{Emax} \colon \ \eta ( x, \boldsymbol \theta ) & = - 1. 8 1 x / ( 0. 7 9 + x ) \\ \text{Sigmoid Emax} \colon \ \eta ( x, \boldsymbol \theta ) & = - 1. 7 0 x ^ { 5 } / ( 4 ^ { 5 } + x ^ { 5 } )$$

The residual error is assumed to be independently normally distributed with mean 0 and variance 4 5. .

Note that despite the fact that the linear and quadratic model are used here as true simulation scenarios, we will always only use the sigmoid Emax model to update the design and analyse the simulated studies. So for the linear and quadratic scenarios we investigate the behaviour of the procedure under misspecification of the dose-response model.

## 3.3 Sample size and timing of interim analysis

Two total sample sizes are used in the simulations: 150 and 250 patients. These sample sizes are consistent with those commonly used in neuropathic pain dosefinding studies (Bornkamp et al., 2007). The timing of the interim analysis was varied in order to consider the impact of sample size in the first stage. Seven different timings for the interim analysis were considered for each total sample size (Table 2).

Figure 3: Dose-response profiles

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-32709a02-d2cb-4a83-a6c0-9ffd2b0a3814.jpg)

Table 2: First stage sample size

|          |   Sample size in first stage N (1) |   Sample size in first stage N (1) |   Sample size in first stage N (1) |   Sample size in first stage N (1) |   Sample size in first stage N (1) |   Sample size in first stage N (1) |   Sample size in first stage N (1) |
|----------|------------------------------------|------------------------------------|------------------------------------|------------------------------------|------------------------------------|------------------------------------|------------------------------------|
| N = 150: |                                 15 |                                 38 |                                 60 |                                 83 |                                105 |                                128 |                                150 |
| N = 250: |                                 15 |                                 54 |                                 93 |                                133 |                                172 |                                211 |                                250 |

A fixed design is also included for comparison to the adaptive designs. This corresponds to the scenario when the sample size for the first stage is 150 or 250 for N = 150 and N = 250 respectively. For the fixed designs, all study participants are assigned to pre-specified doses according to the starting design and there is no updating of the design. Two types of updating were considered

at the interim analysis; ML and Bayesian updating. For each of the 448 scenarios outlined above, 5 , 000 trials were simulated.

## 3.4 Simulation metrics

To ensure comparability of the two approaches for design updating we used ML estimation at study end to estimate the dose-response curve in both situations. That means also in the situation where Bayesian updating was used at the interim analysis, at study end ML was used to estimate the dose-response model parameters. The parameter bounds for θ , θ 3 4 as discussed in Section 2.3.1 will be used for fitting the dose-response model.

## 3.4.1 Efficiency of the second stage design

The relative efficiency of the second stage design is defined as

$$\text{RelEff} & = \frac { | \mathbf M ( \xi ^ { ( 2 ) }, \theta ) | } { | \mathbf M ( \xi ^ { * }, \theta ) | }$$

This is the ratio of the determinant of the information matrix evaluated at the second stage design ξ (2) to the determinant of the information matrix evaluated at the locally optimal design ξ ∗ , assuming that the true parameter values θ are known. The relative efficiency is bound between 0 and 1. The closer the relative efficiency is to 1, the closer the second stage design is to truly optimal.

Whilst the true sigmoid Emax model parameter values are known for the Emax and the sigmoid Emax model scenario, best fitting 'pseudo-true' parameter values must be calculated for the linear and quadratic model scenarios. This is achieved by using least-squares regression to fit the sigmoid Emax model to the true mean function from the linear and quadratic dose-response models at the doses x . The best fitting sigmoid Emax model parameters for the

linear model were θ = ( -0 0396 . , -4 305 12 1 349) and for the quadratic model . , , . θ = ( -0 06617 . , -1 661 1 823 1 948). . , . , .

## 3.4.2 Mean absolute error

As a second metric we will investigate the mean absolute error (MAE) for estimating the dose-response curve (as in Bornkamp et al. (2007)). This is defined as

$$\text{MAE} & = \frac { 1 } { q } \sum _ { i = 1 } ^ { q } | \eta ( x _ { i }, \hat { \boldsymbol \theta } ) - \eta _ { t r u e } ( x _ { i } ) |,$$

where ˆ θ is the ML estimate at study end, and η true ( x ) the true dose-response model function in this scenario. The mean absolute error is calculated for both the adaptive (MAE adapt ) and non-adaptive designs (MAE fixed ). For the nonadaptive design the timing of the 'interim' occurs after all N patients have been allocated, i.e. there is only one stage.

The MAE ratio is defined as

$$R \text{at} _ { \text{MAE} } = \frac { \text{MAE} _ { \text{fixed} } } { \text{MAE} _ { \text{adapt} } }$$

and can be used to compare the adaptive design to the non-adaptive, with a ratio greater than 1 suggesting that the adaptive design performs better.

## 4 Results

## 4.1 Efficiency of the second stage design

The distribution of efficiencies of the calculated second stage design are illustrated in Figure 4. One can observe that in most of the situations the mean efficiency of Bayesian updating outperforms the ML updating, in particular for early interim

analyses, for the sigmoid Emax model scenario and the 'bad' starting design C.

For all dose-response scenarios it can be seen that the efficiency curves for the Bayesian designs only slightly increasing sometimes almost flat, while for ML-based updating the dependency on the size of the first stage is much larger, which might be due to the acknowledgment of parameter uncertainty.

Another interesting observation is that the variability of the efficiencies for the calculated design are considerably smaller for Bayesian designs compared to the designs calculated based on ML. So using a design criterion that acknowledges uncertainty, results in less variable second stage designs.

## Updating ML Bayes

Figure 4: Mean (solid line) with 10% and 90% quantiles (dashed lines) of the observed relative efficiencies of ML and Bayesian updating. N = 250

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-ec8f9808-0a82-4cb1-a5e5-d8979d7b1eae.jpg)

The results for N = 150 are similar and can be found in Appendix C. In Table 3 one can observe the efficiency of the used starting designs, compared to the locally optimal ones. One can see from the Table that, apart from design C, a few of the designs have an efficiency of &gt; 0 8, so that not too much is gained from .

the adaptive designs, when comparing the efficiencies obtained from the adaptive designs in Figure 4. In other words, if the starting design already performs well, then there is less efficiency to be gained by updating the design at an interim analysis.

Table 3: Efficiencies of the different starting designs compared to the corresponding local D optimal designs for the sigmoid Emax model

|          |   Linear |   Quadratic |   Emax |   Sig Emax |
|----------|----------|-------------|--------|------------|
| Design A |     0.91 |        0.61 |   0.62 |       0.73 |
| Design B |     0.89 |        0.92 |   0.79 |       0.58 |
| Design C |     0.22 |        0.03 |   0.19 |       0.12 |
| Design D |     0.81 |        0.76 |   0.63 |       0.86 |

## 4.2 Mean Absolute Error

The ratio of the mean absolute error of the fixed versus the adaptive design is illustrated for each starting design when the sample size N = 250, in Figure 5.

## Updating ML Bayes

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-594a3e06-f7a7-436a-ad52-ea95d6105b6d.jpg)

Number of Patients in First Stage

Figure 5: Ratio of the MAE: N = 250

As can be seen in Figure 5, Bayesian updating seems to outperform ML updating of the design, also in this metric. This happens in particular for early interim analyses. When the size of the first stage increases Bayesian and ML updating become more similar. This is a quite general conclusion across all scenarios.

It can be seen that the benefit of adaption depends strongly on the starting

design. For starting design C one gains in all scenarios, while in other situations the gain is less apparent. For example, for starting design A one gains something by adapting in the Emax scenario and not much for the sigmoid Emax scenario, while for starting design B it is the other way round. This is also consistent with the efficiencies of the starting design in Table 3, all scenarios where a benefit of adaption could be observed in Figure 5 correspond to scenarios, where the starting design has an efficiency of less than 0.7.

The MAE ratio plots corresponding to the smaller sample size of 150 can be found in Appendix C. The findings for N = 150 are similar to those presented above for N = 250.

## 5 Conclusions

From the simulation performed it seems that Bayesian updating of parameter estimates at interim analysis can be used to account for parameter uncertainty. Assigning functional uniform priors to the parameters has been shown to increase the efficiency of the second stage design when compared to ML updating methods, at least in the scenarios studied. The improvement was particularly pronounced for early interim analyses. It would be reasonable to assume that if accurate and valid informative priors could be derived for the parameters, then the improvement may increase further.

An interesting side-result of the simulations is that the benefit of including an interim analysis has not been shown to universally to improve the performance of a dose-finding study. This result is consistent with other simulation studies (Miller et al., 2007; Dragalin et al., 2010; Bornkamp et al., 2011; Jones et al., 2011). Factors such as the signal to noise ratio at interim, or the efficiency of the starting design (as discussed in Dette et al. (2013)) need to be taken into

account. For the rather realistic scenarios discussed in this paper, not much benefit could be expected, as the starting designs, with exception of design C, all already performed relatively well. So most benefit is expected for strongly misspecified starting designs, such as design C. Practical experience suggests that such a scenario, where the increasing part of the dose-response curve is missed in the starting design, is actually not so rare. So in cases of great uncertainty on the dose-response curve (at trial design stage) an adaptive design is a valuable option to consider.

A potentially promising modification of the interim updating strategy, could be 'update-on-demand' adaptive designs. Here one would evaluate the efficiency of the starting design at interim based on current information ( i.e. the current posterior distribution, or ML estimate). One would then only update the design, if the estimated efficiency of the starting design is smaller than some threshold. If the estimated efficiency of the starting design is larger than some threshold one would continue with the starting design until study end.

## A K-means approximation of the full posterior distribution

## A.1 Justification of k-means approximation of posterior distribution

Suppose one would like to approximate an integral

$$\int _ { \ } g ( \theta ) \pi ( \theta ) d \theta,$$

where π ( θ ) is the posterior density, with θ ∈ R d , and a posterior sample θ 1 , . . . , θ T is available.

Suppose also that the function g is expensive to evaluate so that the plain Monte Carlo estimate, 1 T ∑ T t =1 g ( θ t ), is not feasible to evaluate, as T is large.

One approach is to use the k-means clustering algorithm to obtain cluster centers a 1 , . . . , a k ∈ R d with k &lt;&lt; T . This approach works by minimizing

$$\sum _ { i = 1 } ^ { k } \sum _ { \theta _ { t } \in \Omega _ { i } } | | a _ { i } - \theta _ { t } | | ^ { 2 },$$

with respect to the cluster centers a 1 , . . . , a k ∈ R d , where Ω i is given by Ω i = { θ ∈ R d | || θ -a i || = min i ∈{ 1 ,...,k } || θ -a i ||} . That means that every point θ ∈ R d belongs to one set Ω . i The probability of the sets Ω i under π : α i = ∫ Ω i π ( θ ) d θ can be estimated by the proportion ̂ α i of samples within Ω . i

Using these considerations one can come up with a discrete approximation of the posterior distribution of form ˜( π θ ) = ∑ k i =1 ̂ α δ i a i ( θ ), where δ a ( θ ) denotes the function that is equal to 1 if θ = a and 0 otherwise.

Now the idea is to approximate the integral ∫ g ( θ ) π ( θ ) d θ by the corresponding integral with respect to ˜, being given by π ∑ k i =1 ̂ α g i ( a i ).

Now the integration error is given by

$$\text{Now the integration error is given by} \\ & \quad \left | \int g ( \theta ) \pi ( \theta ) d \theta - \sum _ { i = 1 } ^ { k } \widehat { \alpha } _ { i } g ( a _ { i } ) \right | \\ = \ & \left | \sum _ { i = 1 } ^ { k } \int _ { \Omega _ { i } } g ( \theta ) \pi ( \theta ) d \theta - \alpha _ { i } g ( a _ { i } ) + \alpha _ { i } g ( a _ { i } ) - \widehat { \alpha } _ { i } g ( a _ { i } ) \right | \\ = \ & \left | \sum _ { i = 1 } ^ { k } \int _ { \Omega _ { i } } ( g ( \theta ) - g ( a _ { i } ) ) \pi ( \theta ) d \theta + \alpha _ { i } g ( a _ { i } ) - \widehat { \alpha } _ { i } g ( a _ { i } ) \right | \\ \leq \ & \left | \sum _ { i = 1 } ^ { k } \int _ { \Omega _ { i } } ( g ( \theta ) - g ( a _ { i } ) ) \pi ( \theta ) d \theta \right | + | ( \alpha _ { i } - \widehat { \alpha } _ { i } ) g ( a _ { i } ) | \\ \text{The term on the right of Equation (15) will converge almost surely towards}$$

The term on the right of Equation (15) will converge almost surely towards zero when k-means is used to derive a 1 , . . . , a k , see Pollard (1981). As the number of samples T will be large in relation to k , this term will typically converge quite quick towards 0.

The term on the left can be treated with arguments as in (Pag`s, 1997, Propoe sition 7). One can derive an upper bound as

$$\sum _ { i = 1 } ^ { k } \int _ { \Omega _ { i } } | g ( a _ { i } ) - g ( \theta ) | \pi ( \theta ) d \theta \leq L \sum _ { i = 1 } ^ { k } \int _ { \Omega _ { i } } | | a _ { i } - \theta | | \pi ( \theta ) d \theta$$

which holds under the assumption that g is Lipschitz continuous with bound L . Minimizing this upper bound with respect to a 1 , . . . , a k will hence minimize an upper bound of the integration error. In practice of course (16) cannot be optimized directly as π ( θ ) is not available in a closed form. Only a sample θ 1 , . . . , θ T is available. The empirical version of the upper bound in (16) is proportional to

$$\sum _ { i = 1 } ^ { k } \sum _ { \theta _ { t } \in \Omega _ { i } } | | a _ { i } - \theta _ { t } | |,$$

which is very similar to the k-means objective function, see Equation (14). Graf

and Luschgy (2002) consider the speed of convergence of the empirical version. But again as T will typically rather large this will be less of an issue.

## B Simulation scenario

The simulation scenario used to investigate the efficiency of the k-means approximation in terms of relative performance and running time, used as starting design five unequally spaced doses, ξ = (0 1 2 4 8) (Design B) with balanced alloca-, , , , tions. The sigmoid Emax dose-response model was used to simulate the response data with

$$\eta ( x, \boldsymbol \theta ) = \frac { - 1. 7 0 x ^ { 5 } } { ( 4 ^ { 5 } + x ^ { 5 } ) }$$

and an independent normally distributed residual error term with mean 0 and variance 4 5. .

## C Simulation results for N = 150

## Updating ML Bayes

Figure 6: Mean (solid line) with 10% and 90% quantiles (dashed lines) of the observed relative efficiencies of ML and Bayesian updating. N = 250

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-7414c496-d747-4940-9987-c442ffebd7aa.jpg)

## Updating ML Bayes

Figure 7: Ratio of the MAE: N = 150

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0534v2-figure-6769fc5a-2df0-493b-ae08-ca5ebca82a5a.jpg)

## References

Atkinson, A., Donev, A., and Tobias, R. (2007). Optimum Experimental Design, with SAS . Oxford Statistical Science Series. Oxford University Press.

Berry, S., Spinelli, W., Littman, G. S., Liang, J. Z., Fardipour, P., Berry, D. A., Lewis, R. J., and Krams, M. (2010). A Bayesian dose-finding trial with adaptive dose expansion to flexibly assess efficacy and safety of an investigational drug. Clinical Trials , 7:121-135.

Bornkamp, B. (2012). Functional uniform priors for nonlinear modeling. Biometrics , 68:893-901.

Bornkamp, B. (2014). Practical considerations for using functional uniform prior distributions for dose-response estimation in clinical trials. Biometrical Journal , 0:0-0.

Bornkamp, B., Bretz, F., Dette, H., and Pinheiro, J. (2011). Response-adaptive dose-finding under model uncertainty. The Annals of Applied Statistics , 5:16111631.

Bornkamp, B., Bretz, F., Dmitrienko, A., Enas, G., Gaydos, B., Hsu, C., Konig, F., Krams, M., Liu, Q., Neuenschwander, B., Parke, T., and Pinheiro, J. (2007). Innovative approaches for designing and analyzing adaptive dose-finding trials. Journal of Biopharmaceutical Statistics , 17:965-995.

Bornkamp, B., Pinheiro, J., and Bretz, F. (2013). DoseFinding: Planning and Analyzing Dose Finding Experiments , r package version 0.9-6 edition.

Dette, H., Bornkamp, B., and Bretz, F. (2013). On the efficiency of two-stage response-adaptive designs. Statistics in Medicine , 32:1646-1660.

Dragalin, V., Bornkamp, B., Bretz, F., Miller, F., Padmanabhan, S., Patel, N., Perevozskaya, I., Pinheiro, J., and Smith, J. (2010). A simulation study to compare new adaptive dose-ranging designs. Statistics in Biopharmaceutical Research , 2:487-512.

- Dragalin, V., Hsuan, F., and Padmanabhan, S. (2007). Adaptive designs for dosefinding studies based on sigmoid emax model. Journal of Biopharmaceutical Statistics , 17:1051-1070.
- Gelman, A., Carlin, J., Stern, H., and Rubin, D. (2003). Bayesian Data Analysis . Chapman &amp; Hall/CRC.
- Ghosh, J., Delampady, M., and Samanta, T. (2006). An Introduction to Bayesian Analysis: theory and methods . Springer.
- Graf, S. and Luschgy, H. (2002). Rates of convergence for the empirical quantization error. Annals of Probability , 30:874-897.
- Grieve, A. P. and Krams, M. (2005). Astin: a bayesian adaptive dose-response trial in acute stroke. Clinical Trials , 2:340-351.
- Hartigan, J. A. and Wong, M. A. (1979). A k-means clustering algorithm. Journal of the Royal Statistical Society: Series C (Applied Statistics) , 28:100-108.
- Jennrich, R. I. (1969). Asymptotic properties of nonlinear least squares estimation. Annals of Mathematical Statistics , 40:633-643.
- Jones, B., Layton, G., Richardson, H., and Thomas, N. (2011). Model-based bayesian adaptive dose-finding designs for a phase ii trial. Statistics in Biopharmaceutical Research , 3:276-287.

Kirby, S., Brain, P., and Jones, B. (2011). Fitting emax models to clinical trial dose-response data. Pharmaceutical Statistics , 10(2):143-149.

MacDougall, J. (2006). Dose Finding in Drug Development , chapter Analysis of Dose-Response Studies - Emax Model, pages 127-145. Statistics for Biology and Health. Springer.

Miller, F., Guilbaud, O., and Dette, H. (2007). Optimal designs for estimating the interesting part of the dose-effect curve. Journal of Biopharmaceutical Statistics , 17:1097-1115.

Pag` es, G. (1997). A space quantization method for numerical integration. Journal of Computational and Applied Mathematics , 89:1-38.

Pollard, D. (1981). Strong consistency of k -means clustering. Annals of Statistics , 9:135-140.

Pukelsheim, F. (2006). Optimal Design of Experiments . SIAM.

Temple, J. (2012). Adaptive Designs for Dose-Finding Trials . Thesis (doctor of philosophy (phd)), University of Bath.

Thomas, N. (2006). Hypothesis testing and bayesian estimation using a sigmoid emax model applied to sparse dose-response designs. Journal of Biopharmaceutical Statistics , 16:657-677.

Wu, J. (2012). Advances in k-means clustering: a data mining thinking . Springer.
Emma McCallum Björn Bornkamp
Accounting for parameter uncertainty in two-stage designs for Phase II dose-response studies
In this paper we consider two-stage adaptive dose-response study designs,
where the study design is changed at an interim analysis based on the
information collected so far. In a simulation study, two approaches will be
compared for these type of designs; (i) updating the study design by
calculating the maximum likelihood estimate for the dose-response model
parameters and then calculating the design for the second stage that is locally
optimal for this estimate, and (ii) using the complete posterior distribution
of the model parameter at interim to calculate a Bayesian optimal design (i.e.
taking into account parameter uncertainty). In particular, for an early interim
analysis respecting parameter uncertainty seems more adequate, on the other
hand for a Bayesian approach dependency on the prior is expected and an
adequately thought-through prior is required. A computationally efficient
method is proposed for calculating the Bayesian design at interim based on
approximating the full posterior sample using k-means clustering. The sigmoid
Emax dose-response model and the D-optimality criterion will be used in this
paper.