# Point and interval estimation of the target dose using weighted regression
modelling

Saswati Saha
TAGC, Theories and Approaches of Genomic Complexity,
Aix-Marseille Université, Inserm,
Marseille, France.
<EMAIL_ADDRESS>
&Werner Brannath
FB03 Mathematics / Computer Science - Institute for Statistics, Universität,
Bremen, Bremen 28359, Germany.
<EMAIL_ADDRESS>

###### Abstract

In a Phase II dose-finding study with a placebo control, a new drug with
several dose levels is compared with a placebo to test for the effectiveness
of the new drug. The main focus of such studies often lies in the
characterization of the dose-response relationship followed by the estimation
of a target dose that leads to a clinically relevant effect over the placebo.
This target dose is known as the minimum effective dose (MED) in a drug
development study. Several approaches exist that combine multiple comparison
procedures with modeling techniques to efficiently estimate the dose-response
model and thereafter select the target dose. Despite the flexibility of the
existing approaches, they cannot completely address the model uncertainty in
the model-selection step and may lead to target dose estimates that are
biased. In this article, we propose two new MED estimation approaches based on
weighted regression modeling that are robust against deviations from the dose-
response model assumptions. These approaches are compared with existing
approaches with regard to their accuracy in point and interval estimation of
the MED. We illustrate by a simulation study that by integrating one of the
new dose estimation approaches with the existing dose-response profile
estimation approaches one can take into account the uncertainty of the model
selection step.

##  1 Introduction

Proper understanding and characterization of the dose-response relationship is
a fundamental step in the clinical drug development process. A clear
understanding of the dose-response relationship is crucial for addressing two
primary questions in a drug development study: 1) Is there an overall effect
in the drug, which is typically measured by a clinical endpoint of interest?
This is often termed as the Proof of Concept (PoC) step. (ii) If yes, which
dose should be selected for further development? There exist an extensive
literature [Tamhane et al., 1996, Ruberg, 1989, Saha and Brannath, 2018]
addressing the PoC step in dose-response studies. Our focus is more on the
target dose estimation that can be reliably used in the later stages of drug
development. In Phase II clinical trials, often the aim is to obtain the
target dose that produces a clinically relevant effect over the placebo. This
is known as the minimum effective dose (MED). We focus mainly on MED
estimation in this article.

Two common statistical approaches that address the above questions are 1)
Multiple comparison procedures (MCP), 2) Modeling techniques. Under a MCP
approach, the new drug with several active levels is compared with placebo,
and PoC is usually established if at least one of the dose is more effective
than the placebo control [Hochberg and Tamhane, 1987]. Once the PoC is
established, the minimum active dose which gives a statistically significant
and clinically relevant effect compared to the placebo is selected as the
minimum effective dose (MED) [Ruberg, 1995, Tamhane et al., 1996, Tamhane and
Logan, 2002]. Such procedures treat the dose as a qualitative parameter and
make very few, if any, assumptions about the underlying dose-response model. A
drawback of these procedures is that they cannot extrapolate information
beyond the observed dose levels and thereby target dose estimation and
inference are restricted to the observed dose-levels under investigations.
Modeling techniques, on the other hand, establishes the PoC by first fitting a
parametric dose-response model and then comparing it with the ANOVA model with
no dose-response effect. If the PoC is established then the underlying
regression is utilized for MED estimation [Pinheiro et al., 2006, Saha and
Brannath, 2018]. However, the true underlying dose-response shape is not known
before the start of the trial, and specifying one parametric model for dose-
response characterization using modeling techniques often leads to inaccurate
and biased inferences.

Bretz et al., [Bretz et al., 2005] proposed an integrated approach namely MCP-
Mod (Multiple Comparison Procedures-Modeling) that combines multiple
comparisons and modeling techniques to address the above goals in clinical
drug development. The MCP-Mod approach begins with a set of potential dose-
response models, tests for a significant dose-response effect (PoC) using
multiple linear contrasts tests and selects the “best” model using some model
selection criteria among those with a significant contrast test. The “best”
model is then utilized for target dose estimation using regression techniques.
A major disadvantage of the MCP-Mod approach is that it makes prior
assumptions about the values of the model parameters used in the candidate set
which may further lead to a loss in power and unreliable model selection.
Several variations of MCP-Mod approaches [Dette et al., 2015, Gutjahr and
Bornkamp, 2017, Baayen et al., 2015] are also explored over the past few years
which address the apriori parameter value selection issue in the MCP-Mod
approach by replacing the multiple contrast test with the log-likelihood ratio
tests. A brief comparison study on the above approaches [Saha and Brannath,
2018] shows that all the aforementioned approaches can establish the PoC with
strong control of FWER. Note that the above approaches only differ in the PoC
testing step. The model selection and dose selection steps are similar. Even
in the MCP-Mod method, the parameter values are re-estimated in the model
selection step even though they are pre-specified for the contrast test in the
PoC testing step. The above approaches are quite flexible and can address the
model uncertainty to some extent. However, it often happens that they lead to
the selection of an incorrect model. This model may be adequate for an
approximate estimation of the dose-response shape but can be more misleading
when estimating the target dose. Moreover, the true dose-response model may
not at all belong to the candidate set from which the model is selected.

In this article, two new approaches for estimating MED are proposed, which are
robust against model misspecifications. The intention is to fit a weighted
regression model where the weights are so designed such that residuals around
the target dose are given more importance than the residuals away from the
target dose. By doing so, we primarily focus on reducing the bias in
estimating the target dose with a misspecified model at the cost of slightly
inflating the variance. However, it is important to note that a weighted
regression approach is a model-based approach and its’ performance will still
depend on the choice of the dose-response model. Hence, it will often be
advantageous to apply a dose-response analysis strategy with several
parametric candidate models, similar to the ones described in the previous
paragraph. By applying one of the following dose-response analysis approaches;
the MCP-Mod approach [Bretz et al., 2005] or the approach proposed by Dette et
al.,[Dette et al., 2015] or the approach proposed by Baayen et al., [Baayen et
al., 2015] to the dataset, one can simultaneously test the PoC (no dose-
related effect) and detect the dose-response shape that best describes the
data. When the PoC is rejected, the ‘best’ model can be utilized for target
dose estimation using the estimation methods discussed in this article. The
intuition lies in the fact that when the model is misclassified, one can
leverage the idea of using weights to obtain less biased and thereby better
estimates of target doses.

We also propose new confidence intervals for the MED in this article.
Confidence intervals contain a wealth of clinically relevant information that
is generally not available from P values and usual significance testing
[Visintainer and Tejani, 1998]. Confidence intervals estimation for dose-
finding studies have not yet been thoroughly investigated.

An obvious option for interval estimation around the target dose based on
modeling techniques is to use the delta method for obtaining asymptotic
confidence bounds (Wald confidence limits) around the dose-response curve
[Bretz et al., 2005, Dette et al., 2008, Seber and Wild, 1989]. This is
further discussed later in this article, calling this the classical approach
in the rest of the article. The computation of the $100(1-\alpha)\%$ pointwise
Wald confidence interval is simple for linear models. But, dose-response
relationships are usually best described by models that are non-linear in the
parameter as well as the doses. Traditionally, approximate Wald confidence
intervals are used for such non-linear models. But often they yield coverage
rates that deviate from the nominal level [Baayen and Hougaard, 2015]. Often
bootstrap approaches are suggested to improve the coverage probability
[Bornkamp et al., 2009, Pinheiro et al., 2014]. However, fitting non-linear
models are computationally intensive and may result in convergence issues, so
that bootstrap intervals based on a large number of model fits can be complex
and time-consuming. Baayen and Hougaard [Baayen and Hougaard, 2015] proposed a
nice alternative for computing pointwise CI using a profile likelihood
approach. In this article, we propose confidence interval estimates for the
MED with a more robust weighted estimation approach. We compare the resulting
confidence intervals with the classical approach and the two approaches
proposed by Baayen and Hougaard [Baayen and Hougaard, 2015], the percentile
bootstrap and the profile likelihood approaches.

The rest of the article is organized as follows: In Section 2, we will
introduce the parametric set-up, in Section 3, we will give a brief overview
of the classical approach and give details of the new approaches for MED
estimation and inference. Section 3.1 discusses in detail the weight functions
used by the simulations presented later in this article and Section 3.2 and
Section 3.3 describes the two new MED estimation approaches. Section 4
provides a detailed summary of the findings of our extensive simulation study.
The article concludes with a discussion in Section 6.

##  2 Parametric set-up

A parallel dose-group design is considered here with increasing dose levels
$d_{0},\,\ldots,\,d_{k}$ where $d_{0}$ is the placebo dose. The response $Y$
(which can be an efficacy or a safety variable) is assumed to follow a certain
dose-response shape with normally distributed errors:

$\displaystyle Y_{ij}=\mu(d_{i},\bm{\theta})+\epsilon_{ij}=\alpha+\beta
x_{\bm{\gamma}}(d_{i})+\epsilon_{ij},$ (1)
$\displaystyle\quad\epsilon_{ij}\sim N(0,\sigma^{2}),\quad
i=0,\ldots,k,\,j=0,\ldots,n_{i},$

where $Y_{ij}$ refers to the response of the $j\textsuperscript{th}$ patient
in dose group $i$, $n_{i}$ is the number of patients in dose group $i$ and
$n=\sum_{i=0}^{k}n_{i}$ denotes the total sample size.
$\mu(d_{i},\bm{\theta})$ denotes the non-linear regression function and
$\bm{\theta}=\\{\alpha,\beta,\bm{\gamma};\alpha,\beta\in\mathbbm{R},\bm{\gamma}\in\Gamma\subseteq\mathbbm{R}^{p-2}\\}$
refers to the vector of corresponding model parameters, where $\alpha$ and
$\beta$ are the intercept and slope parameters respectively, and $\bm{\gamma}$
is the non-linear parameter capturing the non-linear shape of the model,
$x_{\bm{\gamma}}(d)$. $\Gamma$ denote the parameter space for $\bm{\gamma}$.
Table 1 shows examples of linear or non-linear transformations of the dose
variable $d$ which can be considered in $x_{\bm{\gamma}}(d)$.

Table 1: Linear or non-linear regression functions considered for the dose-
response modelling

Model |  $\mu(\alpha,\beta,{\gamma},d)$ |  $x_{\gamma}(d)$
---|---|---
Linear |  $\alpha+\beta d$ |  $d$
Linear in log-dose |  $\alpha+\beta log(d+c)$ |  $log(d+c)$
Emax |  $\alpha+\beta\frac{d}{ED_{50}+d}$ |  $\frac{d}{ED_{50}+d}$
Exponential |  $\alpha+\beta[exp(\frac{d}{\delta}-1)]$ |  $[exp(\frac{d}{\delta}-1)]$
Quadratic |  $\alpha+\beta(d+\frac{\beta_{2}}{\mid\beta\mid}d^{2})$ |  $d+\frac{\beta_{2}}{\mid\beta\mid}d^{2}$ for $\beta_{2}<0$
Sigmoid Emax |  $\alpha+\beta\frac{d^{h}}{ED_{50}+d^{h}}$ |  $\frac{d^{h}}{ED_{50}+d^{h}}$

##  3 Estimation methods

In this section we introduce two new model-based estimation approaches for the
minimum effective dose (MED).

The MED is defined as,

$MED=\inf\\{d\in(d_{0},d_{k}]\,|\,\mu(d_{i},\bm{\theta})>\mu(d_{0},\bm{\theta})+\Delta\\},$

where $\Delta>0$ denote the clinically relevant effect over the placebo. The
classical approach of estimating the MED by modeling techniques is shown in
Bretz et al., and Dette et al., [Bretz et al., 2005, Dette et al., 2008]. An
estimator of the MED [Bretz et al., 2005] is given by :

$\widehat{MED}=\inf\\{d\in(d_{0},d_{k}]\,|\,\mu(d,\widehat{\bm{\theta}})>\mu(d_{0},\widehat{\bm{\theta}})+\Delta,L_{d}>\mu(d_{0},\widehat{\bm{\theta}})\\},$
(2)

where $\widehat{\bm{\theta}}$ is the least squares estimate of $\bm{\theta}$
and $L_{d}$ is the lower $1-2\alpha$ confidence limit of the estimated mean
value $\mu(d,\widehat{\bm{\theta}})$ at dose $d$ based on the model in
equation (1). Following Dette et al., [Dette et al., 2008], for large sample
sizes $\widehat{MED}$ can be approximated as:

$\displaystyle\widehat{MED}\approx MED(\widehat{\bm{\theta}}),$ (3)
$\displaystyle\textit{ where
}MED(\bm{\theta})=x_{\bm{\gamma}}^{-1}\bigg{(}x_{\bm{\gamma}}(d_{0})+\frac{\Delta}{\beta}\bigg{)}$

with $x_{\bm{\gamma}}^{-1}$ denoting the inverse of the function
$x_{\bm{\gamma}}$ with respect to the dose variable $d$. Asymptotic variance
and asymptotic confidence intervals of $\widehat{MED}$ (equation 2) using the
classical $\delta$-method [Van der Vaart, 1988] are derived and detailed in
Dette et al., [Dette et al., 2008], where the problem of deriving efficient
designs for the estimation of MED are investigated. In this article, we use
the approximate confidence bounds described in Dette et al., [Dette et al.,
2008] as a benchmark for the confidence intervals developed here. The
approximate confidence bounds described by them [Dette et al., 2008] are
referred to as the classical approach in the rest of the article. We want to
obtain robust estimates for the MED using modeling techniques that are less
biased under a misspecified model. With this intention, we make the following
propositions.

In the regression model stated earlier in equation (1), the following sum of
squares:

$\displaystyle SSE=\sum_{i,j}(Y_{ij}-\alpha-\beta x_{{\gamma}}(d_{i}))^{2},$
(4)

are minimized to obtain the least squares estimates. Instead of minimizing the
classical $SSE$ we propose to minimize the weighted SSE ($wSSE$),

$wSSE=\sum_{i,j}w_{ij}(d_{i},MED(\alpha,\beta,\gamma))(Y_{ij}-\alpha-\beta
x_{{\gamma}}(d_{i}))^{2},$ (5)

where the weights depend on the regression estimates, and the MED function
given in equation (LABEL:MED). By doing this we are in a way increasing the
sample size around the target dose with the hope of getting better estimates
of the dose-response curve around the target dose. This can be helpful under
model misspecification when the true model is approximated by a similar model
from the class of dose-response model families. Hence, if the estimation is
more focused on a smaller part of the dose-response curve instead of the
entire curve, one can more accurately get estimates around the point of
interest.

Weighted least squares regression are generally implemented for locally
weighted scatterplot smoothing [Cleveland and Devlin, 1988] or robust
estimation methods in heteroscadestic regression models [Carroll and Ruppert,
1982]. However, the weights considered in the above approaches are either
fixed numbers known apriori or some parametric functions of the predictors
that are not related to the regression parameters. The weights in our approach
are function of the regression parameters and hence the above approaches
cannot be applied to our weighted least squares problem. We consider two
estimation approaches in this article. The first approach performs an iterated
re-weighted least squares regression (IRWLS), where the solution depends on
the optimization algorithm (Newton approach or alternative optimization
algorithm). In the second approach, we consider a robust regression (RR)
approach [Fraiman, 1983], where the weighted least squares estimates are
obtained by solving for $\bm{\theta}$ in the following normal equations:

$\sum_{i,j}(Y_{ij}-\mu(d_{i},\bm{\theta}))\cdot
w_{ij}(d_{i},\bm{\theta})\frac{\partial\mu(d_{i},\bm{\theta})}{\partial\bm{\theta}}=0.$
(6)

This is analogous to the weighted score regression [Antonijevic et al., 2010]
where generalized estimating equations (GEE’s) [Liang and Zeger, 1986] are
used to solve for the above normal equations. We are instead using
’M’-estimation approach proposed by [Huber et al., 1967]. This is further
elaborated in Section 3.3.

Note that our theory is developed based on the assumption that the data
follows a model from the assumed model family. More precisely, with the second
approach we show consistency and asymptotic normality of the weighted least
squares estimates under the assumption that the model class is correctly
specified. We show in Section 4 via simulations that our estimates perform
better than the classical least squares estimates if the model family is
misspecified. We further give a heuristic justification in Appendix 7.1 on how
the estimates obtained by solving the above normal equations (6) can lead to
the estimates which minimize the weighted SSE in equation (5).

###  3.1 Weights

The choice of the weight function is crucial in the above estimation
approaches. The aforementioned RR approach makes some asymptotic inference
that relies on smooth and continuously differentiable weight functions. This
is elaborated in the next section. Even within the class of smooth functions,
a wide range of options are available and the choice of the weights is at the
discretion of the user. We have considered 6 possible weight functions, some
of which are taken from Cleveland and Devlin, [Cleveland, 1979, Cleveland and
Devlin, 1988] and others are experimental. The six possible weight functions
considered in this article are:

$\displaystyle w_{1}(d,d_{MED})$ $\displaystyle=1-|z|^{2}$ $\displaystyle w_{2}(d,d_{MED})$ $\displaystyle=(1-|z|^{2})^{2}$ |  $\displaystyle\text{where }z=\begin{cases}\frac{d_{MED}-d}{d_{MED}}&\quad\text{if }|\frac{d_{MED}-d}{d_{MED}}|<0.9999\\\ 0.9999&\quad\text{o.w}\\\ \end{cases},$
---|---
$\displaystyle w_{3}(d,d_{MED})$ $\displaystyle=1-|z|^{2}$ $\displaystyle w_{4}(d,d_{MED})$ $\displaystyle=(1-|z|^{2})^{2}$ |  $\displaystyle\text{where }z=\begin{cases}\frac{d_{MED}-d}{\underset{d\in\\{d_{1},\ldots,d_{k}\\}}{\min}|d-d_{MED}|}&\,\text{if }|\frac{d_{MED}-d}{\underset{d\in\\{d_{1},\ldots,d_{k}\\}}{\min}|d-d_{MED}|}|<0.9999\\\ 0.9999&\quad\text{o.w}\\\ \end{cases}$
$\displaystyle w_{5}(d,d_{MED})$ $\displaystyle=1-|z|^{2}$ $\displaystyle w_{6}(d,d_{MED})$ $\displaystyle=(1-|z|^{2})^{2}$ |  $\displaystyle\text{where }z=\begin{cases}\frac{d_{MED}-d}{\underset{d\in\\{d_{(2)},\ldots,d_{(k)}\\}}{\min}|d-d_{MED}|}&\,\text{if }|\frac{d_{MED}-d}{\underset{d\in\\{d_{(2)},\ldots,d_{(k)}\\}}{\min}|d-d_{MED}|}|<0.9999\\\ 0.9999&\quad\text{o.w}.\\\ \end{cases}$

Here $d_{MED}$ is the MED estimate ($MED(\widehat{\bm{\theta}})$) obtained
using equation (LABEL:MED) and $d_{(1)},\ldots,d_{(k)}$ are arranged such that
$d_{(1)}$ is closest dose in $\\{d_{1},\ldots,d_{k}\\}$ to $d_{MED}$ and
$d_{(k)}$ is furthest away from $d_{MED}$. The above weight functions are
illustrated in Figure 1. The rationale behind considering the weight functions
$w_{3}$ to $w_{6}$ is to provide a narrow window around the MED where the
residuals in $wSSE$ (equation (5)) are given maximum weights.

Figure 1: Different weight functions plotted when the underlying true dose-
response model is the following emax model: $0.2+0.7\frac{d}{d+0.2}$.The
weight functions peak at the MED (0.26 here) and fades away, as one moves away
from the MED. This is the MED required to achieve the clinical threshold of
$\Delta=0.4$ in the aforementioned emax model.

We have also considered a discrete weight function. The discrete weight
function is used only in the IRWLS method because the RR approach is only
applicable for continuous weight functions. The discrete weight function
considered in this article is:

$w_{7}(d,d_{MED})=\begin{cases}k_{1}&\quad\text{if
}d=\underset{d\in\\{d_{1},\ldots,d_{k}\\}}{\operatorname*{arg\,min}}|d-d_{MED}|\\\
k_{2}&\quad\text{o.w}\\\ \end{cases}$

where $k_{1}$ and $k_{2}$ are integers with $k_{1}>k_{2}$ and $d_{MED}$ is
same as the estimated MED ($MED(\widehat{\bm{\theta}})$) given by equation
(LABEL:MED).

###  3.2 Iterated re-weighted least squares regression (IRWLS)

Weighted non-linear regression using unconstrained optimization is a common
problem and a broad variety of standard tools are available. Numerical
optimization based on simplex methods like the Nelder-Mead algorithm [Nelder
and Mead, 1965], pattern search algorithms [Torczon, 1997] or gradient descent
algorithm can be used to solve the above. However, the weights involved in a
non-linear regression problem are typically some linear or non-linear
functions that are not dependent on the regression parameters[Lim et al.,
2012]. So far, we have not encountered a situation in literature where the
weights in the regression model are treated as functions of the regression
parameters. In this article, we propose an iterated re-weighted least squares
regression (IRWLS) method by which one can naively obtain estimates for the
weighted least squares estimates, where the weights are expressed as a
function of the regression parameters. The IRWLS method can be outlined as
follows:

  * Step 1.

Start with the un-weighted non-linear least squares estimates of model
equation (1) obtained by minimizing $SSE$ in equation (4). Denote them by
$\widehat{\bm{\theta}}_{0}$ and $\hat{\sigma}$.

  * Step 2.

The estimates for the subsequent iterations ($i=2,\ldots,100$) can be computed
using the follow steps:

    * a.

The parametric form of the weight functions are known apriori. Examples of
such parametric forms are illustrated in Section 3.1. For the $i^{th}$
iteration, the weights can be obtained based on the estimates of
$\widehat{\bm{\theta}}$ obtained in the $(i-1)^{th}$ iteration. We organize
the weights for all the doses in the form of a diagonal matrix
$\bm{W_{\nu\times\nu}}$.

$\bm{W_{\nu\times\nu}}=\begin{pmatrix}w(d_{0},d_{MED})&0&\dots&0\\\
\vdots&\ddots&\vdots&\vdots\\\ 0&0&w(d_{,}d_{MED})&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\dots&w(d_{k},d_{MED})\end{pmatrix},$

where $\nu=\sum_{i}n_{i}=n$ and $w$ are weight functions, examples of which
are provided in Section 3.1 and

$d_{MED}=\inf\\{d\in(d_{0},d_{k}]\,|\,\mu(d_{i},\bm{\theta})>\mu(d_{0},\bm{\theta})+\Delta\\},$

as defined in the previous section. In the subsequent sections we are dropping
the ranks of $\bm{W_{\nu\times\nu}}$ for the sake of convenience.

    * b.

Solve $\tau$ for

$\bm{WY}=\bm{WX}\bm{\tau},$

where $\bm{\tau}=(\alpha,\beta)$ are the slope and intercept parameters
introduced with equation (1) which are updated to the estimate
$\bm{\widehat{\tau}_{i}}$ by the weighted least squares regression estimates
obtained by solving the above. Note that, $X=[1;dose]$ for a linear model and
$X=[1;log(dose)]$ for a log linear model. But in a non-linear regression model
the design matrix $X$ is a function of the non-linear parameters
$\bm{\gamma}$. An algorithm implemented in the R package DoseFinding [Bornkamp
et al., 2014] is used here to solve for the non-linear weighted least squares
regression. A grid of non-linear parameters $\bm{\gamma}$ is first selected
and for each possible choice of $\bm{\gamma}$ in the selected grid, $\alpha$
and $\beta$ estimates are obtained as a solution of the following equation:

$\bm{WY}=\bm{WX}_{\bm{\gamma}}\bm{\tau}\text{ where
}X_{\bm{\gamma}}=[1;\mu(dose,\bm{\gamma})]$

Then the final set of parameters which provide the smallest $SSE$ are selected
for the final estimates:
$\widehat{\bm{\theta}_{i}}=(\widehat{\bm{\tau_{i}}},\widehat{\bm{\gamma_{i}}})$.
$\widehat{\sigma_{i}}$ is calculated with the new estimate
$\widehat{\bm{\theta}_{i}}$.

    * c.

Obtain the estimate of $\widehat{MED_{i}}$ using equation (LABEL:MED) and the
expected mean response at the estimated $\widehat{MED_{i}}$.

  * Step 3.

Check for convergence of the regression parameters $\widehat{\bm{\theta}_{i}}$
based on one of the following convergence criteria ($tol$ denotes the
tolerance limit here):

    * a.

$\mathrel{\bigg{(}}\frac{(\widehat{MED}_{new}-\widehat{MED}_{old})}{\widehat{MED}_{old}}\mathrel{\bigg{)}}^{2}\leq
tol$ or if $\widehat{MED}<=0$

    * b.

$\mathrel{\bigg{(}}\frac{(Resp(\widehat{MED}_{new})-Resp(\widehat{MED}_{old}))}{Resp(\widehat{MED}_{old})}\mathrel{\bigg{)}}^{2}\leq
tol$ or if $\widehat{MED}<=0$

where $tol$ (tolerance) is set to $0.001$ (similar to the tolerance limit used
in obtaining convergence in the least squares estimates using R package
nlstools [Baty et al., 2015]). Note that a and b are two alternative options
for checking convergence. We conduct simulations firstly with criteria a and
then with criteria b. Simulation results in Section 4 are shown with criteria
a only, because it performed better than criteria b in our simulations. If the
method converges we stop, else we go to Step 4.

  * Step 4.

Repeat Step 2-3 till convergence or the maximum iteration is reached. If the
algorithm fails to converge or ends in a non-positive MED value, the initial
un-weighted non-linear least squares estimates in Step 1 are assigned as the
final estimates.

A more formal description of the algorithm can be found in Appendix 7.3. In
the IRWLS approach, elaborated above, the weights are treated as if they are
fixed and known before each iteration. Several authors [Carroll and Ruppert,
1982, Box and Hill, 1974] have proposed robust approaches using iterated least
squares algorithm where the first three steps are similar to the algorithm
proposed by us. But their objectives were to fit a heteroscedastic model to a
regression set-up and propose robust estimates for the parameters under such a
set-up. Our objective is different in the sense that we assume a homoscedastic
model and we intend to estimate the MED that is robust against model
misspecification.

Confidence interval estimates can also be obtained with the above approach.

$\widehat{MED}_{W}\pm
u_{1-a}\frac{\hat{\sigma}_{W}}{\sqrt{n}}b^{T}(\hat{\alpha}_{W},\hat{\beta}_{W},\bm{\hat{\gamma}}_{W}){M_{W}}^{-}b(\hat{\alpha}_{W},\hat{\beta}_{W},\bm{\hat{\gamma}}_{W}).$
(7)

where $\hat{\alpha}_{W},\,\hat{\beta}_{W},\,\bm{\hat{\gamma}}_{W}$ are the
weighted estimates from the IRWLS approach, $\widehat{MED}_{W}$ is the MED
estimated by the IRWLS approach, $\hat{\sigma}_{W}$ is the corresponding
standard deviation obtained with the weighted estimates
($\hat{\alpha}_{W},\,\hat{\beta}_{W},\text{ and }\bm{\hat{\gamma}}_{W}$) and
$b(\alpha,\beta,\bm{\gamma})=\frac{\partial{}}{\partial\bm{\theta}}MED(\bm{\theta})$
denotes the gradient of the function $MED$ with respect to $\bm{\theta}$.
${M_{W}}$ is the information matrix which can be defined as follows:

$M_{W}=\sum_{j=0}^{j=k}w_{j}(d_{j},\theta)g(d_{j},\bm{\theta})g^{T}(d_{j},\bm{\theta})\in\mathbbm{R}^{{p}\times{p}}$

where the weights $w_{j}$ are from the weight functions (illustrated in
Section 3.1) used and
$g^{T}(d,\widehat{\bm{\theta}})=\frac{\partial{\mu(d,\widehat{\bm{\theta}})}}{\partial\bm{\theta}}=\bigg{(}1,x_{\hat{\gamma}}(d),\hat{\beta}\frac{\partial{x_{\hat{\gamma}}(d)}}{\partial\gamma_{1}},\ldots\bigg{)}$
is the gradient of the response function $\mu$ with respect to parameter
$\bm{\theta}$. It should be noted that the above confidence interval can only
be obtained when we are considering continuous weight functions and the
estimation method is converging. Moreover, note that with step 4 in the above
algorithm we are safeguarding the method from non-convergence; we propose to
use the least squares estimates when we fail to observe any gain in the
estimation using the weighted least squares algorithm.

###  3.3 Robust regression (RR)

In this section we will describe a MED estimation approach using robust
regression following the works of Ricardo Fraiman [Fraiman, 1983].

Huber et al.,[Huber et al., 1967] considered the class of M-estimates
($\widehat{\bm{\theta}}_{n}$) of a parameter $\bm{\theta}\in\mathbb{R}^{k}$
given by the solution of the system of equations:

$\sum_{i=1}^{n}\phi(Z_{i},\widehat{\bm{\theta}}_{n})=0$ (8)

where $Z_{1},\ldots,Z_{n}\in\mathbb{R}^{q+1}$ are i.i.d variables with common
distribution $F$,
$\phi:\mathbb{R}^{q+1}\times\mathbb{R}^{k}\rightarrow\mathbb{R}^{k}$ is a
function such that the expected value
$\lambda_{F}(\bm{t})=E_{F}(\phi(Z,\bm{t}))$ exists for all $\bm{t}$, and has a
unique root at $\bm{t}=\bm{\theta_{0}}$. Huber et al., [Huber et al., 1967]
showed that $\widehat{\bm{\theta}}_{n}$ converges to the solution of the
equation:

$\lambda_{F}(\bm{t})=E_{F}(\phi(Z,\bm{t}))=0$ (9)

Fraiman [Fraiman, 1983] extends the results shown by Huber et al., [Huber et
al., 1967] for a general class of estimators in a non-linear regression set
up:

$Y_{i}=g(X_{i},\bm{\theta_{0}})+\epsilon_{i}$ (10)

where $Z_{i}=(Y_{i},X_{i})\in\mathbb{R}^{q+1}$ is a sequence of $i.i.d.$
random vectors with common distribution $F$, $X_{i}\in\mathbb{R}^{q}$,
$Y_{i}\in\mathbb{R}$,
$g:\mathbb{R}^{q}\times\mathbb{R}^{k}\rightarrow\mathbb{R}$ is the regression
function, $\bm{\theta_{0}}\in\mathbb{R}^{k}$ is the true unknown parameter,
$\epsilon_{i}\,$ is independent of $X_{i}$ and has distribution
$\mathbbm{N}(0,\sigma^{2})$. They consider the same class of estimators given
by equation (8) but where equation (9) can have more than one solution. They
derive robust estimators in a non-linear regression model and their objective
was to estimate $\bm{\theta_{0}}$ efficiently subject to a bound on the
estimator’s sensitivity to aberrant data.

We want to apply their method on the dose-response models generalized by
equation (1) but our objective is to derive estimators of $\theta_{0}$ that
will ultimately lead to efficient estimator of MED. In order to describe the
M-estimate of $\bm{\theta_{0}}$ and their asymptotic distribution, Fraiman
[Fraiman, 1983] made the following propositions and assumptions:
Define:

$P(\bm{\theta})=\sum_{i=1}^{n}\phi(Y_{i},X_{i},\bm{\theta})$

Then, the Newton-Raphson (NR) algorithm for (8) is given by:

$\widehat{\bm{\theta}}_{n,k+1}=\bm{\widehat{\theta}}_{n,k}-{\bigg{(}\sum_{i=1}^{n}\frac{\partial\phi(Y_{i},X_{i},\bm{\widehat{\theta}}_{n,k})}{\partial\bm{\theta}}\bigg{)}}^{-1}\sum_{i=1}^{n}\phi(Y_{i},X_{i},\bm{\widehat{\theta}}_{n,k})$
(11)

Further define the estimates

$\widehat{\bm{\theta}}_{n}=\begin{cases}\lim_{k\to\infty}\widehat{\bm{\theta}}_{n,k}&\quad\text{
if limit exists }\\\ \widehat{\bm{\theta}}_{n,0}&\quad\text{o.w}\\\
\end{cases}$

where the iteration is started in some consistent initial estimate
$\widehat{\bm{\theta}}_{n,0}$ where
$\widehat{\bm{\theta}}_{n,0}\to\bm{\theta_{0}}$.

#### Assumptions:

  * H1:

For each value of $Z=(Y,X),t\mapsto\phi(Y,X,\bm{t})$ is a continuous function
and the set $V_{Z}=\\{\bm{t}\in\mathbb{R}^{k}|\phi(Y,X,\bm{t})\,\textnormal{is
not differentiable}\\}$ is a nowhere dense subset of the real line, and if $L$
is a line, the intersection $V\cap L$ is a union of intervals or isolated
points.

  * H2:

$F$ is continuous, $E_{F}(\phi(Y,X,\bm{\theta_{0}}))=0$, and
$E_{F}(\frac{\partial\phi(Y,X,\bm{\theta_{0}})}{\partial\theta})$ is non
singular.

  * H3:

There exist $d_{0}>0$ such that

$\displaystyle\sup_{\|\bm{a}-\bm{\theta_{0}}\|\leq
d_{0}}\|(\phi(Y,X,\bm{a})\|\,,F\textnormal{ is integrable}.$
$\displaystyle\sup_{\|\bm{a}-\bm{\theta_{0}}\|\leq
d_{0}}\|\frac{\partial(\phi(Y,X,\bm{a})}{\partial\bm{\theta}}\|\,,F\textnormal{
is integrable}.$

  * H4:

$\widehat{\bm{\theta}}_{n,0}\to\bm{\theta_{0}}$ a.s.

  * N1:

There exist $d_{0}>0$ such that

$\displaystyle\sup_{\|\bm{a}-\bm{\theta_{0}}\|\leq
d_{0}}\|(\phi(Y,X,\bm{a})\|^{2}\,,F\textnormal{ is integrable}.$
$\displaystyle\sup_{\|\bm{a}-\bm{\theta_{0}}\|\leq
d_{0}}\|\frac{\partial(\phi(Y,X,\bm{a})}{\partial\bm{\theta}}\|^{2}\,,F\textnormal{
is integrable}.$

######  Theorem 1.

Assume H1 to H4. Then we can conclude that

  * (i)

$\widehat{\bm{\theta}}_{n}\to\bm{\theta_{0}}$ a.s.

  * (ii)

$\exists\,n_{0}>0\text{ such that
}\widehat{\bm{\theta}}_{n}=\lim_{k\to\infty}\widehat{\bm{\theta}}_{n,k}$ if
$n\geq n_{o}$ a.s.

######  Theorem 2.

Assume H1 to H4 and N1.
Let $V=E_{F}(\phi(Y,X,\bm{\theta_{0}})\phi(Y,X,\bm{\theta_{0}})^{t})$ and
$A=E_{F}(\frac{\partial\phi(Y,X,\bm{\theta_{0}})}{\partial\bm{\theta}})$.
Then, we have

$\sqrt{n}(\widehat{\bm{\theta}}_{n}-\bm{\theta_{0}})\xrightarrow{D}\mathrel{N}(0,A^{-1}V{A^{-1}}^{t})$

In order to apply the above theorem to our situation we substitute $Y_{i}$ by
$Y_{ij}$ and $X_{i}$ by dose $d_{ij}$ and define:

$\phi(Y_{ij},d_{ij},\bm{\theta})=(Y_{ij}-\mu(d_{ij},\bm{\theta}))\cdot
w_{d_{ij}}(d_{ij},\bm{\theta})\frac{\partial\mu(d_{ij},\bm{\theta})}{\partial\bm{\theta}}$
(12)

where $Y_{ij}$ are the responses and $d_{ij}(=d_{i})$ are the dose levels from
the set-up (1). Furthermore, to align with the model framework of
Fraiman[Fraiman, 1983], it is assumed that the dose is a random variable $D$
and doses are so allocated such that $\frac{n_{i}}{n}\to p_{i},$ where
$p_{i}=P(D=d_{i})$ when $n\to\infty\,\forall\,i\in\\{0,\ldots,k\\}$. The
distribution of $\epsilon_{ij}$ is assumed to be independent from the dose.
For the loss function in equation (12) and smooth weight functions in Section
3.1, one can easily verify that the assumptions $H1$ to $H4$ and $N1$ are
satisfied for the dose-response model in equation (1). Appendix 7.2 gives
further elaborations of how the above assumptions are verified with the above
choice of the loss function in equation (12) and with smooth weight functions
shown in Section 3.1.

The ordinary least squares estimates are treated as our initial estimates and
then a sequence $\widehat{\bm{\theta}}_{n}$ are obtained satisfying equation
(11). Further justifications of why the least squares estimates can be used as
initial estimates in our method are given in the Appendix 7.1. Following
Theorem 1 and Theorem 2, the sequence of $\widehat{\bm{\theta}}_{n}$ is
distributed as follows:

$\sqrt{n}(\widehat{\bm{\theta}}_{n}-\bm{\theta_{0}})\xrightarrow{D}\mathrel{N}(0,A^{-1}V{A^{-1}}^{t})$
(13)

where

$\begin{split}V&=E_{F}[(Y-\mu(d,\bm{\theta_{0}}))^{2}\cdot
w^{2}(d,\bm{\theta_{0}})\cdot h(d,\bm{\theta_{0}})h^{t}(d,\bm{\theta_{0}})]\\\
\end{split}$ (14)

with
$h(d_{i},\bm{\theta})=\frac{\partial\mu(d_{i},\bm{\theta})}{\partial\bm{\theta}}$
and

$\begin{split}A&=E_{F}\bigg{(}\frac{\partial\phi(Y,d,\bm{\theta_{0}})}{\partial\bm{\theta}}\bigg{)}\\\
&=E_{F}\bigg{(}-w(d,\bm{\theta_{0}})\frac{\partial\mu(d,\bm{\theta_{0}})}{\partial\bm{\theta}}{\bigg{(}\frac{\partial\mu(d,\bm{\theta_{0}})}{\partial\bm{\theta}}\bigg{)}}^{t}+(Y-\mu(d,\bm{\theta_{0}}))\cdot\frac{\partial[{w(d,\bm{\theta_{0}})\frac{\partial\mu(d,\bm{\theta_{0}})}{\partial\bm{\theta}}}]}{\partial\bm{\theta}}\bigg{)}\end{split}$
(15)

Finally, the MED estimate with the RR approach can be obtained using the
equation (LABEL:MED) where estimate of $\gamma$ and $\beta$ are components of
the estimates $\widehat{\bm{\theta}}$ obtained by solving for equation (9)
using NR algorithm (equation 11). The asymptotic distribution of
$\widehat{\bm{\theta_{n}}}$ shown in equation (13) can be approximated by
plugging-in the empirical estimates of $V$ and $A$.

Furthermore, the asymptotic distribution of $\widehat{MED}$ can also be
derived by applying the delta method to equation (13) similar to Dette et al.,
[Dette et al., 2008]:

$\sqrt{n}(MED(\widehat{\bm{\theta}}_{n})-MED(\bm{\theta_{0}}))\xrightarrow{\mathbb{D}}\mathrel{N}(0,\nabla
MED(\bm{\theta_{0}})^{{}^{\prime}}\cdot(A^{-1}V{A^{-1}}^{t})\cdot\nabla
MED(\bm{\theta_{0}}))$ (16)

where $\nabla MED(\bm{\theta})$ is the gradient of the function
$MED(\bm{\theta})$ (equation LABEL:MED). We can use the empirical estimate of
the above variance to obtain the confidence intervals for MED with the RR
approach.

To better understand the advantages of using the proposed estimation
approaches in this article, we present here a graphical illustration in Figure
2. We have simulated a dataset from an emax model with parameters shown in
Table 5 and then plotted the effect curve, the curve that models the response
curve with the placebo response subtracted from it. The effect-curve is
estimated with the classical approach and the RR approach (with weight
$w_{5}$) in Figure 2. Since we are interested in the minimum effective dose,
the effect curve makes more sense than the full model curve, which models the
observed response. Confidence intervals for the effect curve are also plotted
in Figure 2. For the desired effect of $\Delta=0.2$, the minimum effective
dose is $0.052$. We observe that the effect curve from the classical,
unweighted least squares regression shows a uniform deviation from the true
curve across the full dose range, however, is less precise than the RR method
in the estimation for the effect close to the true MED. Furthermore, the
confidence interval of the effect curve shrinks for the RR method
(Effect_Curve_CI (Weighted Rgn)) as compared to the classical approach effect
curve (Effect_Curve_CI (Wald)) around the actual MED and widens as the dose
increases away from the MED. Thus our objective of getting more accurate dose-
response estimates around the target dose is fulfilled in the illustrated data
with the weighted regression RR approach.

Figure 2: Fitted Emax model on a data simulated from the Emax model in Table 5
using both the classical approach and RR approach (with weight $w_{5}$). The
solid red line labeled as Effect_Curve (Wald) denote the effect curve and the
dotted black curve labeled as Effect_Curve_CI (Wald) denote the confidence
bounds for the effect curve using the classical approach and the solid blue
line labeled as Effect_Curve (Weighted Rgn) and the dotted green curve labeled
as Effect_Curve_CI (Weighted Rgn) denote corresponding effect curve and
confidence interval for the effect curve using the RR approach. The black
solid line labeled as Effect Curve (Actual) is the effect curve of the emax
model from which the data is simulated. The desired threshold $\Delta=0.2$
needed to obtain the MED is shown by a dotted vertical black line and the
actual $MED=0.052$ is indicated by a dotted horizontal black line.

In the following section, we illustrate the implementation of the above
methods and compare their performance with the help of simulations.

##  4 Simulation Study

An extensive simulation study is conducted to evaluate the performance of the
new methods proposed in this article and benchmark them against existing
approaches. To evaluate the estimation accuracy of the MED, we have compared
the MED estimates from the RR and IRWLS approach with the classical MED
estimates proposed by Bretz et al., [Bretz et al., 2005] (in equation (2)). As
already stated in Section 1, the new estimation approach is proposed to get
MED estimates that are robust not only under the ‘true’ model assumptions but
also under model misspecification. We have conducted two sets of simulations
to evaluate the MED estimation accuracy. In the first set, we have simulated
data from two popular data-generating shapes in clinical trials, the Emax and
sigmoidal Emax models, and estimated the MED both under the true model and a
misspecified model. The choice of the misspecified model is elaborated later
in this section.

In a second set of simulations, we generate data from several data-generating
shapes, apply the MCPMod method to simultaneously test the PoC (no dose-
related effect) and detect the dose-response shape that best describes the
data. If the PoC is rejected, the ’best’ model is utilized for MED estimation
using both the classical MED estimates proposed by Bretz et al., [Bretz et
al., 2005] and the RR approach. MCPMod is considered here instead of other
dose-response strategies [Gutjahr and Bornkamp, 2017, Baayen et al., 2015,
Dette et al., 2015] because it was observed in Saha and Brannath, [Saha and
Brannath, 2018] that the model selection performance of all the dose-response
strategies were similar and the implementation of MCPMod is faster than the
other strategies. Since our primary motivation is to enhance the dose-finding
step of the existing dose-response analysis strategies [Saha and Brannath,
2018], we want to evaluate the performance of our approach when integrated
with the MCPMod approach. RR approach integrated with MCPMod approach is
referred to as the MCPMod-RR approach in the rest of the article. The original
MED estimates used in Bretz et al., [Bretz et al., 2005] is referred to as the
classical MCPMod approach in the rest of the article. In the previous
sections, we have also proposed confidence interval estimates for the MED with
the RR approach and the IRWLS approach. In this section, after evaluating the
MED estimation accuracy, we compare our confidence intervals with the
classical approach [Dette et al., 2008] and the two approaches proposed by
Baayen and Hougaard, [Baayen and Hougaard, 2015], the percentile bootstrap and
profile likelihood approaches. The coverage of the $95\%$ confidence interval
for the MED estimates is compared both under the true model and the
misspecified model for the aforementioned methods.

The study design used for the simulations is similar to the case study
mentioned in Bretz et al., [Bretz et al., 2005]. We investigate five dose
levels $(d=0,\,0.05,\,0.2,\,0.6,\,1)$, with a single endpoint measured per
patient, $Y_{i}\sim\mathcal{N}(\mu(d),\sigma^{2})$. A standard deviation of
$\sigma=0.65$ is used in our simulations. This is consistent with the
estimated residual SD observed for the case study [Bretz et al., 2005]. MED is
estimated corresponding to $\Delta=0.4$ in both the scenarios. The following
section describes the two simulation scenarios for evaluating the accuracy of
the different MED estimation approaches:

#### Simulation Scenario 1:

Two different data generating dose-response shapes are investigated for the
mean responses $\mu$. They are defined in Table 2. They have the property that
at $d=0$ the response value is about $0.2$ and the maximum response is about
$0.8$ i.e., maximum dose effect of about $0.6$. For each dose-response shape
in Table 2, $5000$ simulated datasets were generated. For each dose-response
shape in Table 2, we first fit the true model using the classical and new
approaches. We apply our methods with all weights mentioned in Section 3.1 to
check which weights are performing well. After analyzing the true model fit we
analyze the methods under model misspecifications. For data generated from the
sigmoidal Emax model in Table 2, we fit the data also with a linear and an
emax model using both, the classical and RR approach. The IRWLS approach is
not considered in these simulations because it is not giving good results if
the model is not correctly specified. The linear model is chosen because it
was observed in Saha and Brannath, [Saha and Brannath, 2018] that the
sigmoidal Emax model was often misclassified as linear model in the simulation
results, especially for small sample sizes. We have considered the emax model
as it is a commonly used dose-response model and generally considered in the
candidate set of the dose-response estimation methods[Bretz et al., 2005,
Dette et al., 2015, Baayen et al., 2015, Gutjahr and Bornkamp, 2017]. We
analyze the performance of the RR approach under the true model for sample
size per dose group: $n=25$ and under model mis-specifications for two sample
sizes per dose group: $n=25$ and $n=50$.

Table 2:  Data generating dose-response shapes Model |  Simulated from |  Fitted on
---|---|---
Emax |  $0.2+0.7\frac{d}{0.2+d}$ |  emax
Sigmoidal |  $0.2+0.615\frac{d^{4}}{0.4^{4}+d^{4}}$ |  sigEmax, linear, emax

#### Simulation Scenario 2:

In the second simulation scenario, we have used a similar design as in the
first but the data are generated from seven non-linear models as tabulated in
Table 4. Since the IRWLS approach did not perform as good as the RR approach
in the earlier simulations, it is not considered in the second simulation
scenario. We apply our methods with selected weights, $w_{4}$, $w_{5}$ and
$w_{6}$ from Section 3.1 because the above weights show consistent results in
Simulation Scenario 1. Note that results are shown only for weight $w_{6}$ in
Section 5 as the other weights are giving comparable results. We have also
shown the impact on the MED estimation approaches when the data-generating
shapes included in the candidate set of MCPMod deviate somewhat from the
models in Table 4. To achieve this, we have added some noise to the parameters
of the linear, emax and sigmoidal Emax model. The new set of data generating
shapes with added noise are illustrated in Table 4. We have introduced two
noise variables, $\epsilon_{1}$ and $\epsilon_{2}$ in Table 4 that follows
uniform distribution with ranges $(0,0.01)$ and $(0,0.06)$, respectively. To
summarize, we have applied the classical MCPMod approach and the MCPMod-RR
approach to the data generated from both Table 4 and Table 4. The results are
shown in Section 5.

Table 3:  Data generating dose-response shapes Model |  $\mu_{l}(\alpha,\beta,{\gamma}_{l},d)$
---|---
Linear |  $0.2+0.6d$
Linlog |  $0.74+0.33{log(d+0.2)}$
Emax1 |  $0.2+0.7\frac{d}{0.2+d}$
Emax2 |  $0.1+0.3\frac{d}{0.01+d}$
Sigmoidal Emax |  $0.2+0.6\frac{d^{4}}{0.4^{4}+d^{4}}$
Power |  $0.2+0.6d^{0.5}$
Truncated Logistic |  $0.2+\frac{0.682}{(1+exp(10(0.8-d)))}$

Table 4:  Data generating dose-response shapes Model |  $\mu_{l}(\alpha,\beta,{\gamma}_{l},d)$
---|---
Linear |  $0.2+(0.6+\epsilon_{2})d$
Linlog |  $0.74+0.33{log(d+0.2)}$
Emax1 |  $0.2+(0.7+\epsilon_{1})\frac{d}{0.2+\epsilon_{1}+d}$
Emax2 |  $0.1+0.3\frac{d}{0.01+d}$
Sigmoidal Emax |  $0.2+(0.6+\epsilon_{1})\frac{d^{4}}{(0.4+\epsilon_{1})^{4}+d^{4}}$
Power |  $0.2+(0.6+\epsilon_{2})d^{0.5+\epsilon_{1}}$
Truncated Logistic |  $0.2+\frac{0.682}{(1+exp(10(0.8-d)))}$

In order to evaluate the coverage of the new confidence bounds for the MED, we
conduct a simulation study similar to the ones mentioned in Baayen and
Hougaard, [Baayen and Hougaard, 2015]. Following them, we simulate data for
the same five dose levels ($0,\,0.05,\,0.2,\,0.6,\,1$) from the Emax and
sigEmax models, shown in Table 5.

Table 5:  Data generating dose-response shapes for coverage probability evaluation Model |  Simulated from
---|---
Emax |  $0.32+0.74\frac{d}{0.14+d}$
Sigmoidal |  $0.32+0.66\frac{d^{4}}{0.3^{4}+d^{4}}$

The standard deviation of $\sigma=0.65$ is used like in the earlier scenarios.
For each dose-response shape in Table 5, $2000$ simulations are considered.

The coverage of the $95\%$ confidence interval for the MED (corresponding to
$\Delta=0.3$) are evaluated, both under the true model and misspecified model,
for the following methods: confidence interval using a) the classical approach
[Dette et al., 2008], b) the bootstrap approach, c) the profile likelihood
approach introduced by Baayen and Hougaard, [Baayen and Hougaard, 2015] and d)
the weighted regression approach with the RR and IRWLS method introduced in
this article. A different scenario from earlier is considered in Table 5 for
the evaluation of the coverage probability because we wanted to keep the
simulation set-up the same as the simulation set-up used in Baayen and
Hougaard, [Baayen and Hougaard, 2015]. The bootstrap approach for confidence
intervals can be parametrical or non-parametrical. Under the parametric
bootstrap approach, maximum likelihood estimates from the dose-response model
are first obtained and then samples are generated under the model assumptions
in (equation 1) using these maximum likelihood estimates. The MED is then
obtained from these bootstrap samples using the estimates proposed by Bretz et
al., [Bretz et al., 2005] in equation (2). The lower and upper
$100(1-\alpha)\%$ quantiles of the MED obtained from the bootstrapped samples
are used for the confidence intervals of the MED. However, the above bootstrap
approach did not perform well in our simulation scenarios. So we used the
percentile bootstrap approach suggested in Baayen and Hougaard, [Baayen and
Hougaard, 2015]. It can be outlined as follows: Select a grid of doses
$d_{1},\,\ldots,\,d_{G}$ and for the selected grid, the aim is to obtain
pointwise confidence intervals for the expected outcomes
$\hat{y}=\mu(d_{g},\widehat{\bm{\theta}})$ based on a certain dose-response
model $\mu$ (equation 2). For each dose group $i$, $n_{i}$ samples are drawn
with replacement from the observed data. This is repeated $B$ times to
generate $B$ bootstrap samples ($y_{b}^{*},\,b=\\{1,\ldots,B\\}$). For each
bootstrap sample $b$: 1) fit the dose-response model parameters with the
bootstrap data and the given dose-response model and 2) compute the expected
outcomes at pre-selected grid of doses
($\widehat{y}_{b}^{*}(d_{g}),g=1,\ldots,G$). For each bootstrap sample,
$100(1-\alpha)\%$ pointwise confidence intervals for the expected outcomes
$(\widehat{y}_{(\alpha)}^{*}(d_{g}),\widehat{y}_{(1-\alpha)}^{*}(d_{g}))$ are
computed for all doses. Finally, inverse regression techniques are applied to
the above confidence bounds to obtain the confidence intervals for the MED. We
use the above percentile bootstrap approach to benchmark the performance of
our approach with regard to the coverage probability.

Baayen and Hougaard, [Baayen and Hougaard, 2015] also suggested another
approach where one first reformulate the regression model in equation (1) in
terms of the difference of the expected response to the placebo, i.e.
$\beta^{*}=y^{*}-\theta_{0}$ at a specific known dose $d^{*}$ $>0$, where
$y^{*}$ is the expected response at $d^{*}$. Then they showed how to obtain
profile likelihood confidence intervals for the parameter $\beta^{*}$.
Followed by this, they select a grid of increasing doses
$d_{1},\,\ldots,\,d_{G}$, similar to the earlier approach, and for each dose
$d_{g}\text{ where }g\in\\{1,\ldots,G\\}$ obtain the profile likelihood
confidence interval of $\beta_{g}$, i.e, the corresponding value of
$\beta^{*}$ at dose $d^{*}=d_{g}$. This leads to the confidence interval for
the effect curve using a profile likelihood approach for the selected grid of
doses. To obtain an estimate of the MED, they proposed to take the smallest
dose amongst the pre-selected grid, which achieves an effect of $\Delta$ over
the placebo. To obtain a confidence interval for the MED using the profile
likelihood approach they proposed to take the inverse of the confidence bounds
of the effect curve at response value $\Delta$. More details can be found in
their article [Baayen and Hougaard, 2015].

##  5 Simulation Results

In this section, we present the results of the simulation studies described
above and provide a detailed discussion on the pros and cons of the different
MED estimation methods mentioned in this article.

The boxplot distribution of the estimated MED under the different methods and
different scenarios gives an idea of the MED estimation accuracy.
Additionally, a table is added to each boxplot with measures for the distance
of the estimated MED to the true MED (MED: the minimum dose producing an
improvement of $\Delta=0.4$ over the placebo). Since the true MED changes with
the model, the performance of the estimator $\widehat{MED}_{i}$ is measured in
terms of its relative deviation $R_{i}$ from the true MED, where
$R_{i}=\frac{100(\widehat{MED}_{i}-MED)}{MED}$. The mean, median and
interquartile range (IQR) of $R_{i}$ from $5000$ simulation runs characterize
the relative bias and variability of $\widehat{MED}_{i}$.

For the first set of simulations (data generated from Table 2), the results
are shown for $n=25$. Figure 3 and Figure 4 give the boxplots of estimated MED
and the summary statistics of the relative deviation, where the RR and IRWLS
approach are respectively applied with the true dose-response shape. Figure
3(a) and Figure 4(a) show the corresponding distribution of the MED estimates
for the emax model in Table 2 and Figure 3(b) and Figure 4(b) for the
sigmoidal Emax model. From Figure 3, it is evident that the MED estimates,
under the true model fit with the new approaches, do not deviate much from the
classical approach. However, some of the weights performed better compared to
the others in obtaining the MED estimates. From Figure 3(a) it is observed
that the weights $w_{1}$, $w_{4}$, $w_{5}$ and $w_{6}$ are giving similar or
better estimates than the classical approach when the emax model is fitted
with the RR approach. From Figure 3(b) one can conclude that the weights
$w_{3}$, $w_{4}$, $w_{5}$ and $w_{6}$ are giving similar estimates as the
classical approach when the sigEmax model is fitted with the RR approach. For
the IRWLS approach, the MED estimates have more or less identical distribution
as the classical approach for the sigmoidal Emax model (see Figure 4(b)).
However, in Figure 4(a), we observe that the IRWLS approach leads to biased
estimates of MED for the weights: $w_{1}$, $w_{2}$, $w_{5}$ and $w_{6}$.
Hence, in the evaluation of the performance of the different approaches under
model misspecification we only applied the RR approach. The results are shown
in Figure 5 and 6. Additionally, only the weights $w_{4}$, $w_{5}$ and $w_{6}$
are considered under model mis-specification because these weights show
consistent behavior under the true model fit as can be seen from Figure 3(a)
and 3(b).

Figure 5 shows the distribution of the MED estimates when the RR approach with
weight $w_{6}$ is applied with a linear model to the data simulated from the
sigmoidal Emax model in Table 2. Similarly, Figure 6 shows the distribution of
the MED estimates when the RR approach with weight $w_{5}$ is applied with an
emax model to the data simulated from the sigmoidal Emax model in Table 2. It
is interesting to note from Figure 5 that the weighted linear fit estimates
the MED as accurate as the MED estimated under the true model (sigmoidal
Emax). Furthermore, as the sample size increases from $25$ patients to $50$
patients per dose group, the precision in estimating the correct MED with the
weighted linear fit increases. This is evident from Figure 5(a) and 5(b).
Moreover, with the true model fit (sigmoidal Emax), the performance of the RR
approach is similar to the classical approach (see Figure 5). The weighted
linear model was also fitted with $w_{4}$ and $w_{5}$. But with $w_{4}$, the
RR approach is performing identical to the classical approach and with
$w_{5}$, it is giving similar performance as $w_{6}$. Hence, the respective
plots corresponding to weights $w_{4}$ and $w_{5}$ are added to Appendix 7.3.
From Figure 6 we observe that the RR approach can improve the bias in
estimating the true MED over the classical approach under both the
misspecified (emax) and true (sigmoidal Emax) models. For the latter case,
this is contradictory to the presumed theoretical assumptions that the
classical approach, which is based on the maximum likelihood-based estimates,
already provides the UMVUE estimates for the model parameters. Hence, it is
non-intuitive that the RR approach achieves any improvement over the classical
approach under the true model assumptions. However, our observations can be
explained by the fact that for non-linear models the predictions can be mean-
and median-biased even when the model family has been correctly specified and
that a weighted least squares fit has still the potential to reduce this bias
for the doses of interest (i.e. those close to the target dose). The weighted
emax model is also fitted with weights $w_{4}$ and $w_{6}$ and provided
similar results as $w_{5}$. The corresponding plots are added to the Appendix
7.4. We can infer from above that for mis-specified models, one can reduce the
bias considerably and may inflate the variance negligibly. For the true model
fit, the RR approach either performs similar to the classical approach or
achieves a small reduction in bias at the cost of marginal inflation of the
variance. In summary, there is practically no harm in applying the new method
for MED estimation with an appropriate choice of the weight function.

So far, we have observed that the RR approach succeeds in estimating the MED
more accurately than the existing approaches for data generated from two
popular data generating shapes. With the next set of simulations, we intend to
test if we can enhance the dose-finding step of the existing dose-response
analysis strategies by integrating them with the RR approach.

As stated in Section 4, for simulation scenario 2, data are simulated from
seven non-linear models shown in Table 4 and Table 4 to evaluate the
effectiveness of the RR approach when combined with the MCPMod approach. In
this simulation scenario, the linear, emax and the sigEmax model highlighted
in blue in Table 4 is considered in the candidate set of MCPMod. Since the
simulation scenarios in the two tables are leading to similar results, we are
showing the results only for the scenarios from Table 4 in this section. Note
that the second emax model, emax2, is not desirable as it does not reach the
anticipated efficacy. It is introduced in Saha and Brannath, [Saha and
Brannath, 2018] and here, to evaluate the performance of MCPMod under
unfavourable situations. Data are also simulated from the power and truncated
logistic (tlog) models apart from the models in the candidate set because we
want to check how accurately the MCPMod-RR approach estimates the true MED
under model misspecification. Figure 7 to Figure 10 shows the performance of
the aforementioned approaches for data simulated from the models in Table 4.
As evident from Figure 7, we do not gain much with the MCPMod-RR approach for
data simulated from the linear model. We get the same accuracy as the
classical MCPMod approach but the variance is slightly inflated under the
MCPMod-RR approach. However, amongst models included in the candidate set, we
see a considerable improvement in the MED estimation for the sigmoidal Emax
model with the MCPMod-RR approach, as evident from Figure 9(a). The
improvement is significant both in terms of accuracy and precision. For data
generated from emax model, the MCPMod-RR approach is showing comparable
performance as the original MCP-Mod approach (see Figure 8(a)). Particularly
interesting is the performance for the second emax model (emax2) in Figure
8(b). It is important to mention here that in Saha and Brannath, [Saha and
Brannath, 2018] MCPMod fails to detect the emax2 dose-response shape well and
so the MED estimation with MCPMod was also very poor for the data generated
from the emax2 model. Note that the second emax model (emax2) is quite close
to a constant model with no dose-response effect and the MED occurs at the
placebo dose, $d=0$. We observe in Figure 8(b) that though the MCPMod-RR
approach is not able to accurately estimate the true MED for small sample
sizes, but as the sample size increases the accuracy and precision in
estimating the true MED drastically improves compared to the classical MCPMod
approach. For the models not included in the candidate set; while the accuracy
in estimating the true MED improves for the linear log model (see Figure 7(b))
and the power model (see Figure 9(b)) with the MCPMod-RR approach, the
performance of the classical MCP-Mod approach and the MCPMod-RR approach
remains comparable for the truncated logistic model (see Figure 10). It is
important to point out here that for large sample sizes (sample size $=50$ and
$75$ per dose group), the MCPMod-RR approach is showing improvement or similar
performance as the classical MCP-Mod approach both in terms of accuracy and
precision across all the shapes. For small sample sizes, this is not always
true. Often the MCPMod-RR method estimates the true MED at the cost of
inflating variance for small sample sizes.

After analyzing the estimation performance of the new approaches, the next
goal is to analyze the coverage probability of the confidence intervals for
the MED. In order to evaluate the coverage probability of the different
methods, we conduct simulation studies whose set-up was already discussed
earlier in Section 4. Table 6 gives the coverage probabilities of the
different methods for the Emax model, and Table 7 and 8 give the corresponding
results for the sigmoidal Emax model. Note that the RR and IRWLS methods are
applied with weights $w_{5}$ and $w_{6}$ for evaluating the coverage
probability. However, with the RR approach, we were unable to derive the
confidence bound for the MED with weight $w_{6}$ for the sigEmax data
generating shape in $10\%$ simulations. This is mainly because in many
situations the covariance matrix is not estimable due to singularity issues.
Hence, we show results for the RR implementation only with $w_{5}$ in Table 7
and 8. For IRWLS, results are only shown for weight $w_{6}$ for both the emax
and sigmoidal Emax model because $w_{5}$ is giving similar performance as
$w_{6}$. From the results in Table 6, it is clear that the profile likelihood
and percentile bootstrap approach perform superior to all the other methods
for the Emax model. For sample sizes larger than $25$, the percentile
bootstrap performs even better than the profile likelihood approach. The
Weighted RR approach like the classical approach fails to attain the nominal
value. The IRWLS method does not perform well here. They lead to undercoverage
in most of the situations under the emax model. The poor performance of the
IRWLS approach may be due to the fact that theoretically the asymptotic
covariance of the weighted least squares estimates is given by the covariance
matrix in equation (16) and the approximate asymptotic covariance given by the
IRWLS approach (equation (7)) in Section 3.2 deviates from the covariance in
equation (16). The results in Table 7 show that the percentile bootstrap
approach is attaining the nominal level but the classical approach and RR
approaches show overcoverage for the sigmoidal Emax model for large sample
sizes. This may be due to the flexibility of the sigmoid Emax curve when a
limited number of doses are included in the study design. For a limited number
of dose groups, the study data cannot discriminate very well between different
possible model fits and thereafter end up in very wide confidence intervals in
most of the cases. The IRWLS and profile likelihood approaches are prone to
undercoverage for all sample sizes. However, the coverage of the confidence
interval also depends on the dose allocation design. As discussed in Baayen
and Hougaard, [Baayen and Hougaard, 2015] and Dette et al., [Dette et al.,
2008] the coverage of the confidence bounds is expected to improve if a more
optimal dose range is investigated (optimal in terms of dose allocation). To
investigate this, we ran few more simulations for the sigmoidal Emax model in
Table 5 with dose levels $0,\,0.25,\,0.5,\,0.75,\,\text{and }1$. The coverage
probabilities based on $2000$ simulated datasets are shown in Table 8. The
results clearly improve. The classical approach and the weighted RR approach
both attain the nominal level with the new dose allocation design. The
percentile bootstrap approach is attaining the nominal level similar to the
earlier scenario. But the profile likelihood method is prone to undercoverage
in this simulation scenario as well. The undercoverage of the profile
likelihood confidence interval might be because of the convergence issues
arising in the implementation of the method. We need to investigate more to
understand this properly. The performance of IRWLS method improves but it
fails to attain the nominal level even in this simulation scenario. The
simulations in Table 8 consider an ideal case where the dose-response
relationship is quite clear, and an informative dose range is included in the
study.

We also investigate the coverage probabilities of the different approaches
under model mis-specification. The data are simulated from the sigmoidal Emax
model in Table 5 and fitted using an unweighted and weighted linear model
(with weight ($w_{6}$) as in Figure 5) model and an unweighted and weighted
emax model (with weight $w_{5}$ as in Figure 6). For the unweighted linear and
emax model, confidence intervals are estimated using the classical approach,
bootstrap approach, and profile likelihood approach. For the weighted linear
and emax model, confidence intervals are estimated using the RR approach.
Coverage probabilities of the $95\%$ confidence interval from the linear and
emax fit are shown in Table 9 and 10, respectively. We observe in Table 9 that
the classical approach tend to give overcoverage, whereas the percentile
bootstrap and the profile likelihood approach tends to give undercoverage. The
RR approach attains the nominal level of $95\%$ with the weighted linear fit.
Similarly, Table 10 shows that the RR approach with the weighted emax fit is
performing better than the classical approach, percentile bootstrap approach
and profile likelihood approach with the unweighted emax fit. However, for
large sample sizes even the RR approach is giving slight overcoverage as can
be seen from Table 9 and 10.

(a) MED estimation for data simulated from the Emax Model and fitted with the
Emax model

(b) MED estimation for data simulated from the sigEmax Model and fitted with
the sigEmax model

Figure 3: Boxplot distribution of MED estimated under the classical and RR
approach for the different weights given in Figure 1. This is followed by a
table which summarizes the mean, median and IQR of $R_{i}$ for the 5000
simulation trials. Figure (a) and Figure (b) show the MED estimated with the
true model for the data simulated from the emax model and sigEmax model in
Table 2, respectively. MED_Fitted denotes the MED estimated using the
classical approach and MED_Weights1,$\cdots$,MED_Weights6 denote the MED
estimated using RR approach under then different choice of weights
$w_{1},\cdots,w_{6}$ given in Figure 1.

(a) MED estimation for data simulated from the Emax Model and fitted with the
Emax model

(b) MED estimation for data simulated from the sigEmax Model and fitted with
the sigEmax model

Figure 4: Boxplot distribution of MED estimated under the classical and IRWLS
approach for the different weights given in Figure 1. This is followed by a
table which summarizes the mean, median and IQR of $R_{i}$ for the 5000
simulation trials. Figure (a) and Figure (b) show the MED estimated with the
true model for the data simulated from the emax model and sigEmax model in
Table 2, respectively. MED_Fitted denotes the MED estimated using the
classical approach and MED_Weights1$,\cdots,$MED_Weights6, MED_Weights7 denote
the MED estimated using IRWLS approach under then different choice of weights
$w_{1},\cdots,w_{7}$ given in Figure 1.

(a) MED estimation for data simulated from the sigEmax model and fitted with
the linear and sigEmax model for sample size 25

(b) MED estimation for data simulated from the sigEmax model and fitted with
the linear and sigEmax model for sample size 50

Figure 5: Boxplot distribution of MED estimated under the classical and RR
approach for weight $w_{6}$. This is followed by a table which summarizes the
mean, median and IQR of $R_{i}$ for the 5000 simulation trials. Figure (a) and
Figure (b) show the MED estimated with the true model (sigEmax) and the mis-
specified model (linear) for the data simulated from the sigEmax model in
Table 2 with sample size $25$ and $50$, respectively. MED_Ft_Orig and
MED_Ft_Mis denote the MED estimated using the classical approach with the true
and mis-specified model, respectively. MED_Wt_Orig and MED_Wt_Mis denote the
MED estimated using the RR approach with the true and mis-specified model,
respectively.

(a) MED estimation for data simulated from the sigEmax model and fitted with
the emax and sigEmax model for sample size 25

(b) MED estimation for data simulated from the sigEmax model and fitted with
the emax and sigEmax model for sample size 50

Figure 6: Boxplot distribution of MED estimated under the classical and RR
approach for weight $w_{5}$. This is followed by a table which summarizes the
mean, median and IQR of $R_{i}$ for the 5000 simulation trials. Figure (a) and
Figure (b) show the MED estimated with the true model (sigEmax) and the mis-
specified model (emax) for the data simulated from the sigEmax model in Table
2 with sample size $25$ and $50$, respectively. MED_Ft_Orig and MED_Ft_Mis
denote the MED estimated using the classical approach with the true and mis-
specified model, respectively. MED_Wt_Orig and MED_Wt_Mis denote the MED
estimated using the RR approach with the true and mis-specified model,
respectively.

(a) Evaluating MED estimation accuracy between the classical MCP-Mod approach
and the RR approach combined with MCP-Mod approach for data simulated from the
linear model

(b) Evaluating MED estimation accuracy between the classical MCP-Mod approach
and the RR approach combined with MCP-Mod approach for data simulated from the
linear log model

Figure 7: Boxplot distribution of MED estimated under the classical MCP-Mod
and RR combined with the MCP-Mod approach with weight $w_{6}$. This is
followed by a table which summarizes the mean, median and IQR of $R_{i}$ for
the $2000$ simulation trials. Figure (a) and Figure (b) show the comparison of
MED estimation accuracy between the two approaches for data simulated from the
linear and linear log model in Table 4, respectively, with sample size $10$,
$25$, $75$ and $50$. RR(w6) and MM denote the MED estimated using the RR
integrated with MCP-Mod approach and the classical MCP-Mod approach
respectively.

(a) Evaluating MED estimation accuracy between the classical MCP-Mod approach
and the RR approach combined with MCP-Mod approach for data simulated from the
emax model

(b) Evaluating MED estimation accuracy between the classical MCP-Mod approach
and the RR approach combined with MCP-Mod approach for data simulated from the
emax2 (second emax) model

Figure 8: Boxplot distribution of MED estimated under the classical MCP-Mod
and RR combined with the MCP-Mod approach with weight $w_{6}$. This is
followed by a table which summarizes the mean, median and IQR of $R_{i}$ for
the 2000 simulation trials. Figure (a) and Figure (b) show the comparison of
MED estimation accuracy between the two approaches for data simulated from the
emax and emax2 in Table 4, respectively, with sample size $10$, $25$, $75$ and
$50$. RR(w6) and MM denote the MED estimated using the RR integrated with MCP-
Mod approach and the classical MCP-Mod approach respectively.

(a) Evaluating MED estimation accuracy between the classical MCP-Mod approach
and the RR approach combined with MCP-Mod approach for data simulated from the
sigmoidal Emax model

(b) Evaluating MED estimation accuracy between the classical MCP-Mod approach
and the RR approach combined with MCP-Mod approach for data simulated from the
power model

Figure 9: Boxplot distribution of MED estimated under the classical MCP-Mod and RR combined with the MCP-Mod approach with weight $w_{6}$. This is followed by a table which summarizes the mean, median and IQR of $R_{i}$ for the $2000$ simulation trials. Figure (a) and Figure (b) show the comparison of MED estimation accuracy between the two approaches for data simulated from the sigmoidal Emax and power model in Table 4, respectively, with sample size $10$, $25$, $75$ and $50$. RR(w6) and MM denote the MED estimated using the RR integrated with MCP-Mod approach and the classical MCP-Mod approach, respectively.  Figure 10: Boxplot distribution of MED estimated under the classical MCP-Mod and RR combined with the MCP-Mod approach with weight $w_{6}$ for data simulated from the truncated logistic model. This is followed by a table which summarizes the mean, median and IQR of $R_{i}$ for the $2000$ simulation trials. The figure compares the accuracy in estimating the true MED between the two approaches data simulated from the truncated logistic model in Table 4 with sample size $10$, $25$, $75$ and $50$. RR(w6) and MM denote the MED estimated using the RR integrated with MCP-Mod approach and the classical MCP-Mod approach respectively.  Table 6: Table showing the coverage of $95\%$ confidence interval under the different methods and across different sample sizes per dose group for data simulated from the emax model in Table 5 Sample Size | Methods
---|---
Classical | RR($w_{6}$) | RR($w_{5}$) | IRWLS ($w_{6}$) | PBootstrap | Prof Lik
25 | 0.8292 | 0.8466 | 0.8390 | 0.6720 | 0.9355 | 0.9375
50 | 0.8784 | 0.9207 | 0.9035 | 0.7420 | 0.9500 | 0.9415
75 | 0.8988 | 0.9265 | 0.9235 | 0.7510 | 0.9580 | 0.9525
100 | 0.9095 | 0.9378 | 0.9255 | 0.7675 | 0.9550 | 0.9545
Table 7: Table showing the coverage of $95\%$ confidence interval under the different methods and across different sample sizes per dose group for data simulated from the sigEmax model in Table 5 Sample Size | Methods
---|---
Classical | RR($w_{5}$) | IRWLS($w_{5}$) | PBootstrap | Prof Lik
25 | 0.9545 | 0.9070 | 0.6845 | 0.9550 | 0.9132
50 | 0.9735 | 0.9480 | 0.6375 | 0.9525 | 0.8993
75 | 0.9835 | 0.9715 | 0.6890 | 0.9520 | 0.9049
100 | 0.9865 | 0.9785 | 0.6570 | 0.9510 | 0.9118
Table 8: Table showing the coverage of $95\%$ confidence interval under the different methods and across different sample sizes per dose group for data simulated from the sigEmax model with dose groups $0,0.25,0.5,0.75$ and $1$ in Table 5 Sample Size | Methods
---|---
Classical | RR($w_{5}$) | IRWLS($w_{5}$) | PBootstrap | Prof Lik
25 | 0.9425 | 0.9290 | 0.8015 | 0.9370 | 0.9097
50 | 0.9595 | 0.9485 | 0.8315 | 0.9500 | 0.9410
75 | 0.9550 | 0.9500 | 0.8435 | 0.9500 | 0.9317
100 | 0.9580 | 0.9500 | 0.8565 | 0.9475 | 0.9276
Table 9: Table showing the coverage of $95\%$ confidence interval under the different methods and across two sample sizes per dose group for data simulated from the sigEmax model in Table 5 and fitted on linear model Sample Size | Methods
---|---
Classical | RR($w_{6}$) | PBootstrap | Prof Lik
25 | 0.9800 | 0.9530 | 0.903 | 0.9265
50 | 0.9865 | 0.9665 | 0.891 | 0.9180
Table 10: Table showing the coverage of $95\%$ confidence interval under the different methods and across two sample sizes per dose group for data simulated from the sigEmax model in Table 5 and fitted on emax model Sample Size | Methods
---|---
Classical | RR($w_{5}$) | PBootstrap | Prof Lik
25 | 0.9235 | 0.9470 | 0.8920 | 0.9705
50 | 0.9470 | 0.9690 | 0.8915 | 0.9720

The percentile bootstrap approach performs superior to the other approaches
for all the scenarios where the true underlying model is fitted. Unlike other
approaches, it is not so sensitive to the dose-allocation and attains the
nominal level in all the scenarios. Amongst the new approaches suggested in
this article, the RR approach performs better than the IRWLS approach. It
performs similar to the classical approach when the true model is fitted but
it performs better than the other approaches under model mis-specification.
One disadvantage of the percentile bootstrap approach and profile likelihood
approach is that they are grid-based algorithms and the inference from these
approaches are heavily dependent on the choice of grids. Furthermore, due to
the repeated optimization involved in these grid-based methods, they are quite
time-consuming as compared to the other methods. The advantage of the RR
approach is that it is comparatively much less time-consuming.

##  6 Discussion

While developing a drug, typically the choice of dose is based on a Phase II
dose‐finding trial, where selected doses are included along with the placebo.
The analysis of such dose-finding trial/studies can be classified into two
approaches: separate comparisons of each dose to placebo using the multiple
comparison procedure (MCP) or a model‐based strategy, where a dose-response
model is fitted to the data [Pinheiro et al., 2006, Bates and Watts, 1988].
The first approach is preferred when there are only a few dose levels in a
dose-finding study and typically limited knowledge is available about the
dose-response relationship. But if there exist several dose groups with
multiple observations and we want to infer beyond the observed doses, modeling
techniques should be adopted. In this article we have focused on target dose
estimation, particularly, MED, using model-based strategies. Our inspiration
behind the method proposed in this article lies in our earlier article, [Saha
and Brannath, 2018]. While comparing the dose-response profile estimation
strategies in the above article, we inferred that while the existing
approaches can very well establish the PoC of no dose-related effect and
approximately estimate the true dose-response shape by a sufficiently well-
approximated model, there is still some scope of improvement in the target
dose estimation approach.

Note that in the process of obtaining the approximate dose-response shape in
the existing dose-response profile estimation methods, parameter estimates are
obtained that focus on accurately estimating the entire dose-response shape.
The classical approach uses these parameter estimates and it fails to provide
accurate target dose estimates when the true underlying model is misclassified
[Saha and Brannath, 2018]. This is mainly because, in the course of estimating
the entire dose-response curve accurately, it loses in precision and gains
bias in estimating the dose-response curve in a particular region around the
target dose. In this article, inspired by the above notion, two novel
approaches (IRWLS and RR) are introduced for estimating the MED and its
confidence interval. Both approaches aim to provide solutions to a weighted
least squares minimization, where the doses around the MED are given more
weights than the doses away from the MED. In the first approach, we apply an
iterated re-weighted least squares algorithm to obtain our estimates and in
the second approach, we propose to solve a weighted sum of scores function.
With the latter approach, we obtain consistent and asymptotic normal estimates
under specific regularity conditions, when the model family is correctly
specified. Since the RR approach performs much better in our simulations than
the IRWLS approach, we propose to use the second and discard the IRWLS
approach, which often fails to reach the target coverage probability. We also
observe from our simulations that the RR approach gives results similar to the
classical approach when the true model is fitted. Under model
misspecifications, the RR approach produces less biased estimates of the MED
compared to the classical approach. A limitation of the RR approach is that
the asymptotic theory relies on the assumption of a correctly specified model
and needs to be extended under model misspecifications which is a difficult
task, as, consistency cannot be expected in general. We have conducted
experiments where the RR approach is being applied for MED estimation after
the model selection by MCPMod approach. We have shown with our simulations
that the RR approach indeed enhances the dose-finding step of the MCPMod
approach and it might be advisable to use it for target dose estimation for
future stages of drug development.

Though the presentation of confidence intervals is quite mandatory in the area
of dose-finding, the methodology for this has not been well explored in the
literature. Helms et al., [Helms et al., 2015] studied design with active
control of a new drug with several dose levels and investigated the estimation
of the target dose that led to the same efficacy as the active control. The
authors suggested ways to construct confidence intervals under the assumption
of a linear dose-response curve and normal errors. Dilleen et al., [Dilleen et
al., 2003] investigated a similar design with a similar objective as Helms et
al., [Helms et al., 2015] but they also considered a placebo control group in
their study design. They assumed monotonically increasing dose-response
functions and presented different non-parametric methods to estimate the
monotonic dose-response curve. They also suggested ways to construct
confidence intervals for the dose-response curve using the isotonic approach
from Korn, [Korn, 1982] and bootstrap approaches. But both these approaches
are using a design different from ours. We do not have an active control in
our design and our objective is to estimate a dose that produces the desired
threshold (which is known apriori) over the placebo. Traditionally the
classical approach or MCP approaches were used to address this issue. The
approaches for confidence bounds by Baayen and Hougaard [Baayen and Hougaard,
2015] are something new in this area and they introduced a wide spectrum of
options by which one can estimate the confidence bounds for the target dose of
interest. We have only discussed the percentile bootstrap approach and profile
likelihood approach from their article but they have further implemented two
other bootstrap approaches. Unfortunately, they were not performing well, so
we have not included them in our simulations. To the best of our knowledge, we
did not find any article discussing and comparing the coverage of MED using
modeling techniques.

## Acknowledgements

This work was supported by funding from the European Union’s Horizon 2020
research and innovation program IDEAS (”Improving Design, Evaluation and
Analysis of early drug development Studies”; www.ideas-itn.eu) under the Marie
Sklodowska-Curie grant agreement No 633567.

##  7 Appendix

###  7.1 Justification behind the choice of loss function in the Robust
Regression Approach

In the weighted regression approach introduced in this article the primary
objective is to minimize the following weighted SSE:

$wSSE=\sum_{i,j}w_{ij}(d_{i},MED(\bm{\theta}))(Y_{ij}-\mu(d_{i},\bm{\theta}))^{2}$

Ideally such a least square problem is dealt using an iterated re-weighted
least squares algorithm (IRWLS). It starts with some good initial values
[Ruckstuhl, 2010, Cleveland, 1979]. Using the initial estimates, the weights
are obtained and then with these weights $wSSE$ is minimized. and thereafter
the subsequent estimates are obtained. Typically the IRWLS algorithm involves
an iterative procedure where at each iteration the weights are derived using
the estimates from the earlier iteration and then with these particular
weights the $wSSE$ is minimized to obtain the next(or final) regression
estimates. In detail, for the $n^{th}$ iteration step under the IRWLS
algorithm one need to minimize the following:

$\sum_{i,j}w_{ij}(d_{i},MED(\bm{\hat{\theta}_{n-1}}))(Y_{ij}-\mu(d_{i},\bm{\theta}_{n}))^{2}$
(17)

This is equivalent to solving the following normal equations:

$\sum_{i,j}w_{ij}(d_{i},MED(\bm{\hat{\theta}_{n-1}}))\frac{\partial(Y_{ij}-\mu(d_{i},\bm{\theta_{n}}))^{2}}{\partial\bm{\theta}}=0$
(18)

The IRWLS approach leads to a sequence of $\\{\bm{\hat{\theta}_{n}}\\}$ which
usually converges with respect to some criteria (already elaborated in section
3.2) by minimizing the above in the $n^{th}$ iteration. For the final
convergent step, $\bm{\hat{\theta}_{n}}$ is very close to
$\bm{\hat{\theta}_{n-1}}$ and obtaining solutions to (18) by IRWLS is fairly
similar to solving for the normal equations:

$\sum_{i,j}w_{ij}(d_{i},MED(\bm{\theta}))\frac{\partial(Y_{ij}-\mu(d_{i},\bm{\theta}))^{2}}{\partial\bm{\theta}}=0$
(19)

This justifies the choice of the loss function in (20) of Section 3.3 for the
robust regression approach.

###  7.2 Verifying the assumptions by [Fraiman, 1983] in the Robust Regression
Approach

In the robust regression (RR) approach we have considered the following
function:

$\phi(Y_{ij},d_{ij},\bm{\theta})=(Y_{ij}-\mu(d_{ij},\bm{\theta}))\cdot
w_{d_{ij}}(d_{ij},\bm{\theta})\frac{\partial\mu(d_{ij},\bm{\theta})}{\partial\bm{\theta}}$
(20)

We have substituted $d_{ij}$ by $d_{i}$ since $d_{ij}=d_{i}\forall\,j$. From
our model set-up in section 2,
$Y_{ij}\sim\mathcal{N}(\mu(d_{i},\bm{\theta_{0}}),\sigma^{2})$. Hence, we can
write

$\lambda_{F}(\bm{\theta_{0}})=E_{F}(\phi(Z,\bm{\theta_{0}}))=E_{F}\bigg{(}(Y_{ij}-\mu(d_{i},\bm{\theta_{0}}))\cdot
w_{ij}(d_{i},\bm{\theta_{0}})\frac{\partial\mu(d_{i},\bm{\theta_{0}})}{\partial\bm{\theta}}\bigg{)}=0$

We need to show that $\phi$ defined above satisfies the assumptions $H1-H4$
and $N1$ mentioned in Section 3.3. We verify the assumptions in the subsequent
steps.

  * H1:

The weights can be so chosen (for instance the weights $w_{1},\ldots,w_{6}$
shown in Section 3.1) such that $\phi$ is continuous. $\mu$ and $w_{ij}$ are
continuously differentiable over the range of parameters used in the model. So
assumption $H1$ is satisfied.

  * H2:

$F$ is continuous in the experimental set-up (1).
$E_{F}(\phi(Y,D,\theta_{0}))=0$ is evident, where $D$ is the random dose
variable, introduced in Section 2. Need to show;
$E_{F}(\frac{\partial\phi(Y,D,\theta_{0})}{\partial\theta})$ is non singular.
$E_{F}\bigg{(}\frac{\partial\phi}{\partial\theta}\bigg{)}=E_{F}\bigg{(}-\frac{\partial\mu(D,\theta_{0})}{\partial\theta}\psi_{2}(D,\theta_{0})\bigg{)}+E_{F}\bigg{(}\psi_{1}(Y-\mu(D,\theta_{0}))\frac{\partial\psi_{2}(D,\theta_{0})}{\partial\theta}\bigg{)}$

where $\psi_{1}(Y-\mu(D,\theta_{0}))=(Y-\mu(D,\theta_{0}))$ and
$\psi_{2}(D,\theta_{0})=w(D,\theta)\bigg{(}\frac{\partial\mu(D,\theta)}{\partial\theta}\bigg{)}^{\prime}$.
The second term is 0, when the model is correctly specified and close to 0 if
it is approximated by a fairly similar model. So,

$E_{F}\bigg{(}\frac{\partial\phi}{\partial\theta}\bigg{)}=-\sum_{i,j}\frac{\partial
g(d_{i},\theta_{0})}{\partial\theta}\cdot\psi_{2}(d_{i},\theta_{0})=-\sum_{i,j}w_{i,j}(d_{i},\theta_{0})\frac{\partial
g(d_{i},\theta_{0})}{\partial\theta}\cdot\bigg{(}{\frac{\partial
g(d_{i},\theta_{0})}{\partial\theta_{0}}}\bigg{)}^{\prime}$

Similar to non-linear least square method ([Seber and Wild,]) the product of
the gradient function and its transpose is assumed to be non-singular here.
The weights are also defined in Section 3.1 such that they are always non zero
for all dose groups. So from there it follows that
$E_{F}\bigg{(}\frac{\partial\phi}{\partial\theta}\bigg{)}$ is non-singular in
our set-up.

  * H3:

The third assumption claims that there exist a compact ball of radius
$r_{0}\geq 0$ around the true parameter $\bm{\theta_{0}}$ denoted by
$B(\bm{\theta_{0}},r_{0})$ s.t. $\underset{\bm{\theta}\in
B(\bm{\theta_{0}},r_{0})}{\sup}\|\phi(Y,d,\bm{\theta})\|$ and
$\underset{\bm{\theta}\in
B(\bm{\theta_{0}},r_{0})}{\sup}\|\frac{\partial(\phi(Y,d,\bm{\theta}))}{\partial\bm{\theta}}\|$
are both $F$\- integrable.

It is sufficient to show; $\phi(Y,D,\bm{\theta})$ and
$\frac{\partial(\phi(Y,D,\bm{\theta}))}{\partial\bm{\theta}}$ are bounded in
the compact neighbourhood of $\bm{\theta}$.
We have,

$\phi(Y_{i,j},d_{i},\bm{\theta})=\psi_{1}(Y_{i,j},d_{i},\bm{\theta})\psi_{2}(d_{i},\bm{\theta})$

Since $\mu$ is continuously differentiable with respect to $\bm{\theta}$, its
derivative is bounded in this compact neighbourhood. The weight functions are
so chosen that they are also bounded by 1. So $\psi_{2}$ is bounded (say by
$M$) and we have,

$\mid\phi(Y_{i,j},d_{i},\bm{\theta})\mid\leq
M.\mid\psi_{1}(Y_{i,j},d_{i},\bm{\theta})\mid$

If $\exists\,r_{0}$ s.t. $\|\bm{\theta}-\bm{\theta_{0}}\|\leq r_{0}$ and since
$\psi_{1}$ is a continuous in $\bm{\theta}$, therefore $\exists\,M_{r_{0}}$
s.t.
$\|\psi_{1}(Y_{i,j},d_{i},\bm{\theta})-\psi_{1}(Y_{i,j},d_{i},\bm{\theta_{0}})\|\leq
M_{r_{0}}$. Hence, each $\psi_{1}(Y_{i,j},d_{i},\bm{\theta})$ is also bounded
in that compact neighbourhood. So, $\phi$ is bounded. Now,

$\frac{\partial\phi}{\partial\bm{\theta}}=-\frac{\partial\mu(d_{i},\bm{\theta})}{\partial\bm{\theta}}\psi_{2}(d_{i},\bm{\theta})+\psi_{1}(Y_{i,j},d_{i},\theta)\frac{\partial\psi_{2}(d_{i},\bm{\theta})}{\partial\bm{\theta}}$

The weights functions are so chosen that $\psi_{2}$ is continuously
differentiable. Then $\psi_{2}$ and its’ derivative is bounded in the compact
interval. Using similar reasoning as above, one can also show the derivative
of $\phi$ is bounded. Thus $H3$ is satisfied for the above choice of $\phi$.

  * H4:

Under the true dose-response model, $\widehat{\bm{\theta}}_{n,0}$ are
estimates of the un-wegihted non-linear regression, so it follows,
$\widehat{\bm{\theta}}_{n,0}\to\bm{\theta_{0}}$ a.s.

  * N1:

Follows from the arguments for H3.

Thus Theorem 1 and Theorem 2 are satisfied.

###  7.3 Weighted Regression using IRWLS

The following algorithm elaborates how the IRWLS is executed in our approach
and how the optimal estimates are obtained.

Algorithm 1 Algorithm for estimating MED by IRWLS

1:$iteration\leftarrow 1$

2:$converged\leftarrow FALSE$

3:$Set\,tolerance\,for\,convergence\,in\,tol$

4:_loop_ :

5:while $iteration<maxit$ do

6: if $iteration=1$ then

7: $weights=1$ $\forall$ dose

8: else

9: $weights=w_{i},\,\text{for some weight function in Section \ref{Weights}}$

10: if model being fitted belongs to linear family then

11: if model=linear then X=[1;dose]

12: else if model=linlog then X=[1;log(dose)]

13: Solve for $WY=WX\beta$ using QR decomposition and obtain $\hat{\beta}$.

14: else

15: $bnds\leftarrow default\,band\,limits$ (Dose Finding Package)

16: Grid search (for $\gamma$ in equation (1) of Section 2) in the above

17: band limit to obtain the $\hat{\gamma}$ which minimizes the SSE. Refine
the

18: bounds of $\gamma$ around the estimate obtained by grid search

19: procedure How to refine bounds?

20: $strt\leftarrow estimate\,obtained\,from\,grid\,search$

21: $bnds\leftarrow default\,band\,limits$

22: $N\leftarrow 30$

23: $dif\leftarrow\frac{bnds[2]-bnds[1]}{N}$

24: $bnds[1]\leftarrow\max(c(strt-1.1*dif),bnds[1])$

25: $bnds[2]\leftarrow\min(c(strt+1.1*dif),bnds[2])$

26: With the refined bounds get the optimal estimates of the parameters

27: of the full model, renew weights and check for convergence

28: if
$\mathrel{\big{(}}\frac{(\widehat{MED}_{new}-\widehat{MED}_{old})}{\widehat{MED}_{old}}\mathrel{\big{)}}^{2}<=tol$
or if $\widehat{MED}<=0$ then

29: $converged\leftarrow TRUE$ $\triangleright$ Condition for convergence

30: return $\hat{\beta}$ and $\widehat{MED}$

31: else

32: $iteration\leftarrow iteration+1$, goto _loop_

33:if $iteration=maxit$ AND $converged=FALSE$ then

34: return $\hat{\beta}$= estimates obtained in the first iteration

###  7.4 Boxplot distribution of MED

(a) MED estimation for data simulated from the sigEmax model and fitted with
the linear and sigEmax model for sample size 25

(b) MED estimation for data simulated from the sigEmax model and fitted with
the linear and sigEmax model for sample size 50

Figure 11: Boxplot distribution of MED estimated under the classical and RR
approach for weight $w_{4}$. This is followed by a table which summarizes the
mean, median and IQR of $R_{i}$ for the 5000 simulation trials. Figure (a) and
Figure (b) show the MED estimated with the true model (sigEmax) and the
misspecified model (linear) for the data simulated from the sigEmax model in
Table 2 with sample size $25$ and $50$, respectively. MED_Ft_Orig and
MED_Ft_Mis denote the MED estimated using the classical approach with the true
and misspecified model, respectively. MED_Wt_Orig and MED_Wt_Mis denote the
MED estimated using the RR approach with the true and misspecified model,
respectively.

(a) MED estimation for data simulated from the sigEmax model and fitted with
the linear and sigEmax model for sample size 25

(b) MED estimation for data simulated from the sigEmax model and fitted with
the linear and sigEmax model for sample size 50

Figure 12: Boxplot distribution of MED estimated under the classical and RR
approach for weight $w_{5}$. This is followed by a table which summarizes the
mean, median and IQR of $R_{i}$ for the 5000 simulation trials. Figure (a) and
Figure (b) show the MED estimated with the true model (sigEmax) and the
misspecified model (linear) for the data simulated from the sigEmax model in
Table 2 with sample size $25$ and $50$, respectively. MED_Ft_Orig and
MED_Ft_Mis denote the MED estimated using the classical approach with the true
and misspecified model, respectively. MED_Wt_Orig and MED_Wt_Mis denote the
MED estimated using the RR approach with the true and misspecified model,
respectively.

(a) MED estimation for data simulated from the sigEmax model and fitted with
the emax and sigEmax model for sample size 25

(b) MED estimation for data simulated from the sigEmax model and fitted with
the emax and sigEmax model for sample size 50

Figure 13: Boxplot distribution of MED estimated under the classical and RR
approach for weight $w_{4}$. This is followed by a table which summarizes the
mean, median and IQR of $R_{i}$ for the 5000 simulation trials. Figure (a) and
Figure (b) show the MED estimated with the true model (sigEmax) and the
misspecified model (emax) for the data simulated from the sigEmax model in
Table 2 with sample size $25$ and $50$, respectively. MED_Ft_Orig and
MED_Ft_Mis denote the MED estimated using the classical approach with the true
and misspecified model, respectively. MED_Wt_Orig and MED_Wt_Mis denote the
MED estimated using the RR approach with the true and misspecified model,
respectively.

(a) MED estimation for data simulated from the sigEmax model and fitted with
the emax and sigEmax model for sample size 25

(b) MED estimation for data simulated from the sigEmax model and fitted with
the emax and sigEmax model for sample size 50

Figure 14: Boxplot distribution of MED estimated under the classical approach
and weighted regression method (RR) for weight $w_{6}$. This is followed by a
table which summarizes the mean, median and IQR of $R_{i}$ for the 5000
simulation trials. Figure (a) and Figure (b) show the MED estimated with the
true model (sigEmax) and the misspecified model (emax) for the data simulated
from the sigEmax model in Table 2 with sample size $25$ and $50$,
respectively. MED_Ft_Orig and MED_Ft_Mis denote the MED estimated using the
classical approach with the true and misspecified model, respectively.
MED_Wt_Orig and MED_Wt_Mis denote the MED estimated using the RR approach with
the true and misspecified model, respectively.

## References

  * [Antonijevic et al., 2010] Antonijevic, Z., Pinheiro, J., Fardipour, P., and Lewis, R. J. (2010).  Impact of dose selection strategies used in phase ii on the probability of success in phase iii.  Statistics in Biopharmaceutical Research, 2(4):469–486.
  * [Baayen and Hougaard, 2015] Baayen, C. and Hougaard, P. (2015).  Confidence bounds for nonlinear dose–response relationships.  Statistics in medicine, 34(27):3546–3562.
  * [Baayen et al., 2015] Baayen, C., Hougaard, P., and Pipper, C. (2015).  Testing effect of a drug using multiple nested models for the dose–response.  Biometrics, 71(2):417–427.
  * [Bates and Watts, 1988] Bates, D. M. and Watts, D. G. (1988).  Nonlinear regression: iterative estimation and linear approximations.  Nonlinear regression analysis and its applications, pages 32–66.
  * [Baty et al., 2015] Baty, F., Ritz, C., Charles, S., Brutsche, M., Flandrois, J.-P., Delignette-Muller, M.-L., et al. (2015).  A toolbox for nonlinear regression in r: the package nlstools.  Journal of Statistical Software, 66(5):1–21.
  * [Bornkamp et al., 2009] Bornkamp, B., Pinheiro, J., and Bretz, F. (2009).  Mcpmod: An r package for the design and analysis of dose-finding studies.  Journal of Statistical Software, 29(i07).
  * [Bornkamp et al., 2014] Bornkamp, B., Pinheiro, J., and Bretz, F. (2014).  Dosefinding: planning and analyzing dose finding experiments, 2014.  URL http://CRAN. R-project. org/package= DoseFinding [Accessed on 18 May 2016].
  * [Box and Hill, 1974] Box, G. E. and Hill, W. J. (1974).  Correcting inhomogeneity of variance with power transformation weighting.  Technometrics, 16(3):385–389.
  * [Bretz et al., 2005] Bretz, F., Pinheiro, J. C., and Branson, M. (2005).  Combining multiple comparisons and modeling techniques in dose-response studies.  Biometrics, 61(3):738–748.
  * [Carroll and Ruppert, 1982] Carroll, R. J. and Ruppert, D. (1982).  Robust estimation in heteroscedastic linear models.  The annals of statistics, pages 429–441.
  * [Cleveland, 1979] Cleveland, W. S. (1979).  Robust locally weighted regression and smoothing scatterplots.  Journal of the American statistical association, 74(368):829–836.
  * [Cleveland and Devlin, 1988] Cleveland, W. S. and Devlin, S. J. (1988).  Locally weighted regression: an approach to regression analysis by local fitting.  Journal of the American statistical association, 83(403):596–610.
  * [Dette et al., 2008] Dette, H., Bretz, F., Pepelyshev, A., and Pinheiro, J. (2008).  Optimal designs for dose-finding studies.  Journal of the American Statistical Association, 103(483):1225–1237.
  * [Dette et al., 2015] Dette, H., Titoff, S., Volgushev, S., and Bretz, F. (2015).  Dose response signal detection under model uncertainty.  Biometrics, 71(4):996–1008.
  * [Dilleen et al., 2003] Dilleen, M., Heimann, G., and Hirsch, I. (2003).  Non-parametric estimators of a monotonic dose–response curve and bootstrap confidence intervals.  Statistics in medicine, 22(6):869–882.
  * [Fraiman, 1983] Fraiman, R. (1983).  General m-esttmators and applications to bounded influence estimation for non-linear regression.  Communications in Statistics-Theory and Methods, 12(22):2617–2631.
  * [Gutjahr and Bornkamp, 2017] Gutjahr, G. and Bornkamp, B. (2017).  Likelihood ratio tests for a dose-response effect using multiple nonlinear regression models.  Biometrics, 73(1):197–205.
  * [Helms et al., 2015] Helms, H.-J., Benda, N., and Friede, T. (2015).  Point and interval estimators of the target dose in clinical dose-finding studies with active control.  Journal of biopharmaceutical statistics, 25(5):939–957.
  * [Hochberg and Tamhane, 1987] Hochberg, Y. and Tamhane, A. (1987).  Multiple comparison methods.
  * [Huber et al., 1967] Huber, P. J. et al. (1967).  The behavior of maximum likelihood estimates under nonstandard conditions.  In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1, pages 221–233. Berkeley, CA.
  * [Korn, 1982] Korn, E. (1982).  Confidence bands for isotonic dose-response curves.  Applied Statistics, pages 59–63.
  * [Liang and Zeger, 1986] Liang, K.-Y. and Zeger, S. L. (1986).  Longitudinal data analysis using generalized linear models.  Biometrika, 73(1):13–22.
  * [Lim et al., 2012] Lim, C., Sen, P. K., and Peddada, S. D. (2012).  Accounting for uncertainty in heteroscedasticity in nonlinear regression.  Journal of statistical planning and inference, 142(5):1047–1062.
  * [Nelder and Mead, 1965] Nelder, J. A. and Mead, R. (1965).  A simplex method for function minimization.  The computer journal, 7(4):308–313.
  * [Pinheiro et al., 2014] Pinheiro, J., Bornkamp, B., Glimm, E., and Bretz, F. (2014).  Model-based dose finding under model uncertainty using general parametric models.  Statistics in medicine, 33(10):1646–1661.
  * [Pinheiro et al., 2006] Pinheiro, J. C., Bretz, F., and Branson, M. (2006).  Analysis of dose–response studies—modeling approaches.  In Dose finding in drug development, pages 146–171. Springer.
  * [Ruberg, 1989] Ruberg, S. J. (1989).  Contrasts for identifying the minimum effective dose.  Journal of the American Statistical Association, 84(407):816–822.
  * [Ruberg, 1995] Ruberg, S. J. (1995).  Dose response studies i. some design considerations.  Journal of Biopharmaceutical Statistics, 5(1):1–14.
  * [Ruckstuhl, 2010] Ruckstuhl, A. (2010).  Introduction to nonlinear regression.  IDP Institut fur Datenanalyse und Prozessdesign, Zurcher Hochschule fur Angewandte Wissenschaften. See: http://www. idp. zhaw. ch.(Cited on p. 365.).
  * [Saha and Brannath, 2018] Saha, S. and Brannath, W. (2018).  Comparison of different approaches for dose response analysis.  Biometrical Journal.
  * [Seber and Wild, ] Seber, G. and Wild, C.  Nonlinear regression. 2003.  Hoboken: John Wiley & Sons Google Scholar.
  * [Seber and Wild, 1989] Seber, G. and Wild, C. (1989).  Nonlinear regression john wiley & sons.  New York.
  * [Tamhane et al., 1996] Tamhane, A. C., Hochberg, Y., and Dunnett, C. W. (1996).  Multiple test procedures for dose finding.  Biometrics, pages 21–37.
  * [Tamhane and Logan, 2002] Tamhane, A. C. and Logan, B. R. (2002).  Multiple test procedures for identifying the minimum effective and maximum safe doses of a drug.  Journal of the American Statistical Association, 97(457):293–301.
  * [Torczon, 1997] Torczon, V. (1997).  On the convergence of pattern search algorithms.  SIAM Journal on optimization, 7(1):1–25.
  * [Van der Vaart, 1988] Van der Vaart, A. (1988).  Estimating a real parameter in a class of semiparametric models.  The Annals of Statistics, pages 1450–1474.
  * [Visintainer and Tejani, 1998] Visintainer, P. F. and Tejani, N. (1998).  Understanding and using confidence intervals in clinical research.  The Journal of Maternal-fetal Medicine, 7(4):201–206.