Emax Model

Consider a clinical dose finding study for an anti-anxiety drug (Dette et al (2010) for details) with the endpoint being the change from baseline score at study end. Assume that the placebo effect is at \(\theta_0=0\) and the maximum treatment effect within the dose range [0, 150]mg is 0.4. The goal is to estimate the smallest dose achieving 100p% of the maximum effect in the observed dose range.

A dose-response relationship that is increasing and has a maximum effect which is achieved asymptotically at large dose levels can be modelled by the Emax model. This is reflected by modelling the response \(y(x)\) at a dose \(x\) by the EMAX model

\[ y(x)=\eta(x,\theta)=\theta_0+\frac{\theta_1x}{x+\theta_2}+\varepsilon(x), \]

where \(\theta_0\) denotes the placebo effect at dose 0, \(\theta_1\) is the asymptotic maximum treatment benefit over placebo and \(\theta_2\) is the dose that gives half of the asymptotic maximum effect.

From previous studies, the initial parameter values that will be used to compute the locally optimal designs, are \(\theta=(0, 0.467, 25)'\).

Using this onformation, we can create the matrix of candidate regressors using the function Fx_dose

F.emax <- Fx_dose(0:150, c(0,0.467,25))
## [1] Call of the function:
## Fx_dose(dose.levels = 0:150, theta0 = c(0, 0.467, 25))

\(ED_p\)-optimal Design

Now, our focus will be on computing approximate \(ED_p\)-optimal design for in the Emax model.

The \(ED_p\), \(0 < p < 1\), is the smallest dose achieving 100p% of the maximum effect in the dose range \([a, b]\) and formally, it can be defined as \[ ED_p=\mathrm{arg}\min_{x\in[a,b]}\left\{ \frac{\eta(x,\theta)-\eta(a,\theta)}{\eta(x_{\max}, \theta)-\eta(a,\theta)}\geq p\right\}, \] where \(x_{\max}\) denotes the dose at which the maximum expected response is observed. Using the reasoning in Dette et al (2010), we have \(x_{\max}=b=150\) and \(a=0\), and we can express the \(ED_p\) as a function of \(\theta\): \[ ED_p=\frac{ab+\theta_2((1-p)a+pb)}{\theta_2+pa+(1-p)b}=\frac{150p\theta_2}{\theta_2+150(1-p)}=:\beta(\theta). \]

For a maximum likelihood estimate \(\hat{\theta}\), \(\beta(\hat{\theta})\) is an estimator of \(ED_p\) and its asymptotic variance is \[ Var(\beta(\hat{\theta}))\approx c'(\theta)M^{-}(\xi,\theta)c(\theta), \] where \(c'(\theta)=\frac{\partial}{\partial\theta}\beta(\theta)(0,0,k)\) for some \(k\in\mathbb{R}\).

Using the expression for \(\beta(\theta)\) in our setting, we can see that \[ c(\theta)=\frac{150^2p(1-p)}{(\theta_2-150p+150)^2}(0,0,1)'. \]

Thus, we have expressed the problem of finfind \(ED_p\)-optimal design as a problem of finding \(c\)-optimal design with the vector \(c\) given by \(c(\theta)\) above, which can be computed using the function od_REX:

p <- 0.5
k <- (1-p) * p * 150^2 / (25 + p*150 + 150)^2
app.emax <- od_REX(F.emax, crit="c", h=c(0,0,k))
## [1] Call of the function:
## od_REX(Fx = F.emax, crit = "c", h = c(0, 0, k))
## Linear programming reformulation will be used for c-optimality.
## [1] Method od_c_LP finished after 0.05 seconds at 2020-03-15 10:58:50 .
## [1] c-criterion value: 4.70821994067889e-09

In practise, the so-called standard design is often used \[ \xi_s=\left\{ \begin{array}{cccccc} 0&10&25&50&100&150\\ \frac{1}{6}&\frac{1}{6}&\frac{1}{6}&\frac{1}{6}&\frac{1}{6}&\frac{1}{6} \end{array} \right\} \]

xi.s <- rep(0,151)
xi.s[1] <- 1/6; xi.s[11] <- 1/6; xi.s[26] <- 1/6; xi.s[51] <- 1/6; xi.s[101] <- 1/6; xi.s[151] <- 1/6

We can easily compute the efficiency of this design with respect to the \(ED_p\)-optimal design:

optcrit(F.emax, xi.s, crit="c", h=c(0,0,k), echo=FALSE) / app.emax$Phi.best
## [1] 0.6416512

Now, consider a practical situation where 30 patients are available for our trial. In this case, we need to compute exact \(ED_p\)-optimal design of size 30.

ex.emax <- od_MISOCP(F.emax, 30, crit="c", h=c(0,0,k))

Now, the efficiency of the standard design with respect to the exact 30-point \(ED_p\)-optimal design is

optcrit(F.emax, 30*xi.s, crit="c", h=c(0,0,k), echo=FALSE) / ex.emax$Phi.best
## [1] 0.6427764

Homework

Repeat the previous analysis for a log-linear model \(y(x)=\theta_0+\theta_1\log(x+\theta_2)+\varepsilon(x)\) and the initial parameter values \(\theta=(0, 0.0797, 1)\).

References

Dette H, Kiss C, Bevanda M, Bretz F (2010). Optimal designs for the EMAX, log-linear and exponential models. Biometrika, 97(2), 513-518.