Sage Journals: Discover world-class research

Abstract

The Binary Emax model is widely employed in dose–response analysis during drug development, where missing data often pose significant challenges. Addressing nonignorable missing binary responses—where the likelihood of missing data is related to unobserved outcomes—is particularly important, yet existing methods often lead to biased estimates. This issue is compounded when using the regulatory-recommended ‘‘imputing as treatment failure’’ approach, known as non-responder imputation (NRI). Moreover, the problem of separation, where a predictor perfectly distinguishes between outcome classes, can further complicate likelihood maximization. In this paper, we introduce a penalized likelihood-based method that integrates a modified expectation-maximization (EM) algorithm in the spirit of Ibrahim and Lipsitz to effectively manage both nonignorable missing data and separation issues. Our approach applies a noninformative Jeffreys’ prior to the likelihood, reducing bias in parameter estimation. Simulation studies demonstrate that our method outperforms existing methods, such as NRI, and the superiority is further supported by its application to data from a Phase II clinical trial. Additionally, we have developed an R package, ememax (https://github.com/Celaeno1017/ememax), to facilitate the implementation of the proposed method.

Keywords

Dose–response Emax model EM algorithm noignorable missing Firth correction separation

1. Introduction

The dose–response relationship is a fundamental aspect of research in various applied statistics fields, particularly in clinical trials and bioinformatics. In real-world studies, incomplete data records are common due to reasons such as nonresponse to questionnaires, typographical errors, loss to follow-up, and data contamination. While missing data is often unavoidable, a common approach when the number of observations is large is to perform analysis using only complete records—a method known as complete case analysis (CC). However, in many scenarios, the application of CC yields biased estimates when analyzing dose–response relationships with binary outcomes.^1,2 Moreover, small samples—common in randomized controlled trials—can exacerbate its variance and instability.

The challenge of missing data has been a significant area of research for decades, with numerous statistical methods developed to address it. Various missing data mechanisms, such as missing completely at random (MCAR) and missing at random (MAR), have been extensively studied, and corresponding methodologies have been proposed. It is well known that if the missing data mechanism is not appropriately modeled, parameter estimates may be biased.³ While most related research has focused on scenarios where missing covariates and response values are assumed to be ignorable, limited attention has been given to scenarios with nonignorable missing response values, particularly within the context of modeling dose–response relationships. When the likelihood of missingness conditional on covariates depends on the unobserved values of the response, this is referred to as nonignorable missing data.⁴ In this article, we explore the implications of nonignorable missing data within the context of a dose–response model with binary outcomes.

The sigmoid Emax model is commonly used in clinical trials to explore the binary dose–response relationship. Let $n$ denote the total sample size, and suppose $y_{i}$ denotes the binary outcome for the $i$ -th patient randomized to a dose ${Dose}_{i}$ , where ${Dose}_{i}$ is one of the predefined dose levels from a set of $J$ levels ${D_{1}, \dots, D_{J}}$ , for $i = 1, \dots, n$ . Without loss of generality, let $P (y_{i} = 1 | {Dose}_{i}) = π_{i}$ be the probability of success for patient $i$ after dosage. Under this setup, the sigmoid four-parameter Emax model can be written as the following

l o g (\frac{π_{i}}{1 - π_{i}}) = E_{0} + \frac{E_{\max} \times {Dose}_{i}^{λ}}{{ED}_{50}^{λ} + {Dose}_{i}^{λ}}

(1)

where

E_{0}

is the expected logit of dose effect at

{Dose}_{i} = 0

, with

{Dose}_{i} = 0

often being considered as placebo;

E_{\max}

is the expected logit of the maximum achievable effect (at infinite dose);

{ED}_{50}

is the dose that produces half-maximal effect

E_{\max} / 2

; and

λ

is the slope factor or Hill parameter, determining the steepness of the dose–response curve. Even though the four-parameter Emax model is available in the literature, according to meta-analysis studies,⁵ oftentimes in practice the three-parameter Emax model fits well on most real-world data, where the Hill parameter is assumed to be

1

.^6,7 Hence by putting

λ = 1

the Emax model is reduced to

\log (\frac{π_{i}}{1 - π_{i}}) = E_{0} + \frac{E_{\max} \times {Dose}_{i}}{{ED}_{50} + {Dose}_{i}} .

(2)

Define

θ = (E_{0}, {ED}_{50}, E_{\max})^{⊤}

, then from (2), the success probability

π_{i}

for patient

i

can be written as

\begin{aligned} P (y_{i} = 1 | θ, {Dose}_{i}) & = \exp (E_{0} + \frac{E_{\max} \times {Dose}_{i}}{{ED}_{50} + {Dose}_{i}}) / (1 + \exp (E_{0} + \frac{E_{\max} \times {Dose}_{i}}{{ED}_{50} + {Dose}_{i}})) . \end{aligned}

(3)

Given a dataset, finding the estimate

\hat{θ}

of the parameter

θ

is the main point of interest in clinical trials for exploring dose–response characteristics. Once

\hat{θ}

is obtained, one can easily estimate the success rates of different doses and draw statistical inferences of the corresponding dose populations. The estimate

\hat{θ}

of the parameter

θ

can be obtained by maximizing the likelihood (actually log-likelihood) given in the following:

L (θ, Dose) = \prod_{i = 1}^{n} f (y_{i} | θ, {Dose}_{i}) = \prod_{i = 1}^{n} {p (y_{i} = 1 | θ, {Dose}_{i})}^{y_{i}} {1 - p (y_{i} = 1 | θ, {Dose}_{i})}^{1 - y_{i}} .

(4)

In clinical trials, it is common for binary outcomes to be missing for some patients receiving either an experimental or placebo dose. Given the typically small sample sizes in these studies, complete case analysis is often impractical. A standard practice, endorsed by the FDA and other regulatory agencies, is to impute all missing values as ”treatment failures.” This method, known as non-responder imputation (NRI),⁸ is widely used but has well-documented limitations. NRI can introduce substantial bias into estimates and inferences, particularly when the missing data mechanism is nonignorable.⁹

Another commonly used method for addressing missing data is multiple imputation (MI). This approach involves a two-stage process: first, generating several plausible imputed datasets based on assumed data distributions, and second, combining the results using a predefined pooling rule.¹⁰ However, MI typically assumes that the missingness mechanism does not depend on the unobserved data, an assumption that is violated under nonignorable missingness.¹¹ Consequently, MI can also result in biased estimates and inferences. Some modified MI methods have been proposed to address missing not at random (MNAR) scenarios by model selection approach^12,13 or pattern-mixture model approach,¹⁴ but these methods often rely on assumptions of imputation model, especially the difference between distributions of respondent and non-respondent. This limits their flexibility in practical applications, where variables associated with missingness can vary widely. Moreover, the performance of MI heavily depends on the chosen pooling methods after generating and analyzing imputed datasets, which can be challenging to determine in practice.¹⁵ Therefore, developing new approaches to fit the Emax model with incomplete data, particularly when the missing data mechanism is nonignorable, remains a critical issue in medical research.

In small or medium sample size settings for fitting binary outcome models, a common issue that can arise is the nonconvergence of estimates, a phenomenon known as “separation.”¹⁶ Complete separation occurs when a single covariate or a linear combination of covariates perfectly predicts the outcome, leading to divergence in the estimation process. Even if complete separation does not occur, the presence of quasi-complete separation—where a subset of subjects’ responses is perfectly predicted—can still cause estimation challenges.¹⁷ These situations, though not uncommon in biomedical datasets, are often overlooked. For example, current methods and packages that use likelihood estimation for dose–response or Emax models, such as the ClinDR package in R by Thomas¹⁸ and the Dosefinding package by Bornkamp et al.,¹⁹ do not account for the issue of separation. Moreover, it is well known that estimation of Emax model under finite small sample sizes may encounter unstable convergence due to its strictly increasing assumption, which further complicated the model inference with binary response.²⁰

In the framework of generalized linear models, Heinze and Schemper²¹ and Heinze²² demonstrated that Firth’s method,¹ initially designed to reduce the bias of maximum likelihood estimates (MLE), can also effectively address the problem of separation. However, their work did not account for the presence of missing data. Extending this approach, Maiti and Pradhan,²³ and Maity, Pradhan and Das²⁴ applied Firth’s method to reduce bias and adjust for separation under a nonignorable missing data mechanism within the logistic regression framework.

When fitting the binary Emax model, two distinct but related separation issues can arise. The first type of separation occurs in the binary outcome y, even in the absence of any missing data. This happens when a predictor or a linear combination of predictors perfectly separates the outcome classes, leading to estimation challenges. The strictly increasing assumption of the Emax model can particularly make it prone to this type of separation, especially with small sample sizes. The second type of separation is induced by the missingness mechanism itself. In the presence of nonignorable missing data, the missing data pattern may be perfectly predicted by model variables. This further complicates the estimation process, as the missingness mechanism now introduces an additional source of separation.

In the field of dose–response relationships, no prior research has addressed the issue of nonignorable missing data within the Emax model framework. Additionally, there is no existing literature on bias reduction and separation in the Emax model, even with complete data. The key contribution of this paper is the integration of the Firth-type bias correction with the weighted EM algorithm proposed by Ibrahim and Lipsitz to address both types of separation issues within the binary Emax model framework. While previous work has dealt with the issue of nonignorable missing data and Firth-type bias correction in the context of logistic regression, our paper is the first to extend this approach to the more complex nonlinear Emax model.

The remainder of the article is organized as follows: Section 2 outlines the proposed approach and details the derivation of the estimation algorithm. Section 3 explores the simulation settings and presents the results. Section 4 discusses a real-world application, providing insights into the practical implementation of the proposed method. Finally, Section 5 offers a discussion of the findings and suggests potential directions for future research. Additionally, we have developed an R package, ememax (https://github.com/Celaeno1017/ememax), that implements our methods.

2. Method

In this section, we present the proposed methodology in detail, beginning with the model formulation and assumptions that address these challenges.

2.1. Weighted EM procedure of Ibrahim and Lipsitz (IL)

Let $r$ be the missing indicator vector whose $i$ -th element is defined as

r_{i} = {\begin{cases} 1 & if y_{i} is missing \\ 0 & if y_{i} is observed \end{cases}

and is generated by

P (r_{i} = 1 | z_{i}, α) = p_{i} = \frac{\exp (z_{i}^{⊤} α)}{1 + \exp (z_{i}^{⊤} α)}, i = 1, 2, \dots, n

where

z_{i} = (x_{i}^{⊤}, {Dose}_{i}, y_{i})^{⊤}

is the covariate vector with

x_{i}

to be a

(p + 1) \times 1

covariate vector of interest including an intercept term, and

α = (α_{0}, \dots, α_{p}, α_{(p + 1)}, α_{(p + 2)})^{⊤}

is the corresponding parameter vector. If

α_{(p + 2)} = 0

, the missing data mechanism does not depend on

y_{i}

, and hence the missing mechanism is ignorable. However, if

α_{(p + 2)} \neq 0

, the missing mechanism depends on

y_{i}

and is therefore nonignorable. Note that when

α

is a null vector, the missing mechanism is MCAR. For fitting binary regression model with nonignorable missing values, Ibrahim and Lipsitz²⁵ proposed an EM algorithm to compute the estimate of the regression coefficient. Following Ibrahim and Lipsitz combined with the binary Emax model, the joint likelihood can be written as,

f (r, y | α, θ, Dose, x) = \prod_{i = 1}^{n} f (y_{i} | θ, {Dose}_{i}) f (r_{i} | z_{i}, α),

(5)

where

f (r_{i} | z_{i}, α) = {P (r_{i} = 1 | z_{i}, α)}^{r_{i}} {1 - P (r_{i} = 1 | z_{i}, α)}^{1 - r_{i}} .

(6)

The maximum likelihood of

(α, θ)

of (5) can be obtained via the EM algorithm by maximizing the expected log-likelihood. The E-step, expectation w.r.t. missing outcome conditional on observed data, of the

i

-th individual’s contribution can be written as

\begin{aligned} E [l (θ, α | z_{i}, r_{i})] = {\begin{cases} \sum_{y_{i} = 0}^{1} l (θ, α | z_{i}, r_{i}) f (y_{i} | r_{i}, {Dose}_{i}, x_{i}, θ, α) & if y_{i} is missing \\ l (θ, α | z_{i}, r_{i}) & if y_{i} is observed \end{cases} \end{aligned}

(7)

where

f (y_{i} | r_{i}, {Dose}_{i}, θ, α)

is the conditional distribution of the missing outcome given the observed data, and

l (θ, α | z_{i}, r_{i})

is the log-likelihood of the

i

-th individual given in (5). In the above expectation, let

w_{i y_{i}} = f (y_{i} | r_{i}, {Dose}_{i}, x_{i}, θ, α)

, which can be considered as weight. This can be written further using Bayes theorem as

\begin{aligned} f (y_{i} | r_{i}, {Dose}_{i}, x_{i}, θ, α) = w_{i y_{i}} & = \frac{f (r_{i}, y_{i} | α, θ, x_{i}, {Dose}_{i})}{\sum_{y_{i} = 0}^{1} f (r_{i}, y_{i} | α, θ, x_{i}, {Dose}_{i})} \\ = \frac{f (y_{i} | θ, {Dose}_{i}) f (r_{i} | z_{i}, α)}{\sum_{y_{i} = 0}^{1} f (y_{i} | θ, {Dose}_{i}) f (r_{i} | z_{i}, α)} \end{aligned}

(8)

Therefore, for all data with

i = 1, \dots, n

, the

(t + 1)

-th iteration of the E-step can be expressed as

\begin{aligned} Q (θ, α | θ^{(t)}, α^{(t)}) & = \sum_{i = 1}^{n} \sum_{y_{i} = 0}^{1} w_{i y_{i}}^{(t)} l (θ, α | z_{i}, r_{i}, {Dose}_{i}, y_{i}) \\ = \sum_{i = 1}^{n} \sum_{y_{i} = 0}^{1} w_{i y_{i}}^{(t)} {l (θ | {Dose}_{i}, y_{i}) + l (α | z_{i}, r_{i})} \end{aligned}

(9)

where the

t

-th stage weights of E-step are defined as

\begin{aligned} w_{i y_{i}}^{(t)} & = {\begin{cases} \frac{f (y_{i} | {Dose}_{i}, θ^{(t)}) f (r_{i} | z_{i}, α^{(t)})}{\sum_{y_{i} = 0}^{1} f (y_{i} | {Dose}_{i}, θ^{(t)}) f (r_{i} | z_{i}, α^{(t)})} & if y_{i} is missing \\ 1 & if y_{i} is observed. \end{cases} \end{aligned}

(10)

Since (9) is the sum of two equations involving parameter

α

for the logistic model and parameter

θ

for the Emax model, and both parameters

α

and

θ

are independent by assumption, for M-step, we can maximize these separately. Details about the derivation of the M-step can be found in supplementary Section 1. Finally, by repeating the E-step and M-step until convergence, we can get the estimates of

θ

and

α

2.2. Firth-type bias reduction on Ibrahim and Lipsitz (FIL)

Firth¹ introduced a modification to the score function of the likelihood to reduce the bias of the maximum likelihood estimator (MLE) in small sample settings. Firth demonstrated that when the parameter in question is the canonical parameter of a full exponential family—such as in the logistic regression model for the missing indicator model with $r$ —the modification of the score function is equivalent to applying a Jeffreys’ invariant prior to the likelihood function. In the spirit of Firth, to obtain explicit Firth-type bias-reduced estimates of the parameter $θ$ for the Emax model with missing outcome, the joint likelihood corresponding to (5) can be modified as follows:

f^{*} (r, y ∣ α, θ, x, Dose) = \prod_{i = 1}^{n} f^{*} (y_{i} ∣ θ, {Dose}_{i}) f^{*} (r_{i} ∣ z_{i}, α),

where the Firth-type bias reduction is achieved by penalizing each likelihood component on the right-hand side by multiplying it with the Jeffreys’ invariant prior as the penalty term. The penalized log-likelihood function for modeling nonignorable missingness is:

l^{*} (α ∣ r, z) = l (α ∣ r, z) + \frac{1}{2} | I (α) |,

(11)

where

| I (α) |

is the determinant of the observed information matrix, and

I (α) = Z^{T} V Z

, with Z as the design matrix and

V = diag (p_{i} (1 - p_{i}))

. The penalized log-likelihood for the Emax model part is given by:

l^{*} (θ ∣ y, Dose) = l (θ ∣ y, Dose) + \frac{1}{2} | I (θ) |,

(12)

where

I (θ) = \sum_{i = 1}^{n} (1 - π_{i}) π_{i} \nabla η {({Dose}_{i}, θ)}^{⊤} \nabla η ({Dose}_{i}, θ) .

The joint penalized log-likelihood can be obtained by combining (11) and (12), utilizing the independence of

θ

and

α

. The bias-reduced MLE can then be obtained by maximizing this joint penalized log-likelihood.

Subsequently, we apply IL method with penalized joint log-likelihood, and the E-step of the EM becomes

Q^{*} (θ, α | θ^{(t)}, α^{(t)}) = \sum_{i = 1}^{n} \sum_{y_{i} = 0}^{1} w_{i y_{i}}^{(t)} {l^{*} (θ | y_{i}, {Dose}_{i}) + l^{*} (α | r_{i}, z_{i}))} .

(13)

The modification in the E-step as shown above leads to a change in the M-step to obtain the maximizer as well. We present the detail derivation in the supplementary Section 1.

Finally, the estimator can be obtained via any iterative procedures such as Newton–Raphson, Gauss–Newton, and Fisher scoring. Details about computational considerations, including initial value selection and computational time, are discussed in the Supplementary material.

We summarize all the steps of the IL/FIL algorithm as Algorithm 1.

2.3. Variance estimator and confidence interval

Let $γ = (θ, α)$ denote the parameter for estimation via EM, and $\hat{γ}$ denote the final estimator obtained from EM. The approximate variance–covariance matrix of $\hat{γ}$ for both IL and FIL methods can be estimated via observed information matrix $I (γ)$ . In the EM setting, $I (γ)$ can be estimated using Louis²⁶ as below:

\begin{aligned} I (γ) = & - H (γ) - E [S (γ | y, r, Dose, z) S (γ | y, r, Dose, z)^{⊤}] \\ + E [S (γ | y, r, Dose, z)] E {[S (γ | y, r, Dose, z)]}^{⊤} \end{aligned}

where

S (γ | y, r, Dose, z)

is the complete-data score vector, and all the expectations are with respect to the conditional distribution of missing outcome given the observed data. Incorporating with the IL or FIL setting, we define

\begin{aligned} S_{i} (γ) & = \frac{\partial}{\partial γ} l (γ | y, r, Dose, z) \\ {\dot{q}}_{i} (γ) & = \sum_{y_{i} = 0}^{1} w_{i y_{i}} S_{i} (γ) . \end{aligned}

Thus, the estimated observed information matrix of

γ

giving observed data is

I (\hat{γ}) = - H (\hat{γ}) - \sum_{i = 1}^{n} \sum_{y_{i} = 0}^{1} {\hat{w}}_{i y_{i}} S_{i} (\hat{γ}) S_{i} (\hat{γ})^{⊤} + \sum_{i = 1}^{n} {\dot{q}}_{i} (\hat{γ}) {\dot{q}}_{i} (\hat{γ})^{⊤}

(14)

where

{\hat{w}}_{i y_{i}}

is the estimate at convergence of EM. Note that all the quantities in (14) can be obtained easily from the M-step as byproducts. We can get a consistent estimator of variance–covariance matrix as the inverse of

I (\hat{γ})

With the variance–covariance matrix estimator $I (\hat{γ})^{- 1}$ , confidence intervals for the $γ$ can be constructed. The standard errors of parameters can be obtained by taking the square root of the diagonal elements of the variance–covariance matrix, and the 100(1- $α$ )% confidence interval based on the normal approximation for parameter $γ_{i}$ is:

({\hat{γ}}_{i} - z_{α / 2} {\hat{s}}_{γ_{i}}, {\hat{γ}}_{i} + z_{α / 2} {\hat{s}}_{γ_{i}})

(15)

where

{\hat{γ}}_{i}

is the estimator of

γ_{i}

z_{α / 2}

is the

1 - α / 2

quantile of the standard normal distribution, and

{\hat{s}}_{γ_{i}}

is the estimated standard error of

{\hat{γ}}_{i}

3. Simulation study

To evaluate the performance of the proposed methodology, we conducted a thorough simulation study based on a general Phase-II dose–response clinical trial setting. We investigated the effect of sample size on the estimation process by considering sample sizes of n = 150, 250, 350, and 450. The total sample size was evenly distributed across five different treatment dose arms: Dose = (0, 7.5, 22.5, 75, 225), ensuring equal sample sizes in each treatment arm. The success rate for response in the placebo arm $({Dose}_{1} = 0)$ was set at 10%, the maximum success rate with an infinite dose at 80%, and the dose achieving a half-maximal effect at 7.5. This corresponds to $E_{0} = logit (0.1)$ , $E_{\max} = logit (0.8) - logit (0.1)$ , and ${ED}_{50} = 7.5$ .

The response variable $y_{i}$ was generated from a Bernoulli distribution with success probability $π_{i}$ as defined in (2), for $i = 1, \dots, n$ . After generating $y_{i}$ , the nonignorable missing data indicator $r_{i}$ was generated using a Bernoulli distribution with a missing probability $p_{i}$ , where $p_{i} = \exp (z_{i}^{⊤} α) / (1 + \exp (z_{i}^{⊤} α))$ . The covariate vector was defined as $z_{i} = (1, x_{i 1}, x_{i 2}, {Dose}_{i}, y_{i})^{⊤}$ , and $α = (α_{0}, α_{1}, α_{2}, α_{3}, α_{4})^{⊤}$ . We sampled $x_{i 1}$ and $x_{i 2}$ from two independent standard normal distributions. If $r_{i}$ was generated as 1, the corresponding $y_{i}$ value was masked as ‘NA’, indicating a missing response. The nonignorable nature of the missing mechanism was ensured by setting $α_{4} \neq 0$ . The overall rate of missing responses was controlled by selecting appropriate values for $α$ . For instance, to achieve an overall missing rate of approximately 15%, we set $α = (0.5, - 3, 0, - 0.05, 1)^{⊤}$ . A noteworthy aspect is that $α_{2}$ was set to 0 to introduce a degree of model mis-specification when using the full $z_{i}$ vector to predict missingness.

Table 1 summarizes the simulation results comparing five different estimation methods: Complete Case (CC), NRI, MI, Ibrahim–Lipsitz (IL), and Firth-Adjusted IL (FIL), based on $N = 1000$ replications. For each of the simulated data with missing values, MI was performed with $m = 100$ imputations using the mice package, employing predictive mean matching and lasso-logistic regression imputation methods along with the Namard–Rubin pooling rule. The CC and NRI analyses were conducted using the ClinDR package. For point estimation, we report the average estimated value of $θ_{i}$ , $\hat{θ_{i}} = (1 / s) \sum {\hat{θ_{i}}}^{(k)}$ for $i = 1, 2, 3$ corresponding to the three-parameter Emax model, where ${\hat{θ_{i}}}^{(k)}$ is the estimate of $θ_{i}$ in the $k$ -th replication, and $s$ is the number of valid estimates out of $N$ replications. Additionally, we report the mean bias error (MBE) as $MBE = (1 / s) \sum ({\hat{θ_{i}}}^{(k)} - θ_{i})$ , the root mean squared error (RMSE) as $RMSE = \sqrt{(1 / s) \sum ({\hat{θ_{i}}}^{(k)} - θ_{i})^{2}}$ , and the mean estimated standard error for ${\hat{θ}}_{i}$ , ${\hat{s}}_{θ_{i}} = (1 / s) \sum {\hat{s}}_{θ_{i}}^{(k)}$ . For confidence interval estimation, we report the coverage probability CP at 95% confidence level using $CP = (1 / s) \sum I ({\hat{θ_{i}}}^{(k)})$ , where $I (\cdot)$ is the indicator function for whether $θ_{i}$ falls within the estimated confidence interval. We also report the mean estimated interval length (Est. length) as $Est. length = (1 / s) \sum λ ({\hat{θ_{i}}}^{(k)})$ , where $λ ({\hat{θ_{i}}}^{(k)})$ represents the length of the confidence interval for ${\hat{θ_{i}}}^{(k)}$ .

Table 1.
Estimates, mean bias error, root mean squared error, estimated standard errors, coverage probabilities, and 95% Wald confidence intervals based on 1000 simulations with missing rate $\approx$ 15%. The best values for each metric are marked as bold.

Sample Size(N) Parameter Type Estimate MBE RMSE Est.SE CP Est.Length

150 log( ${ED}_{50}$ ) CC 2.104 0.089 0.679 0.698 0.979 2.736

NRI 2.564 0.549 0.858 0.666 0.896 2.610

MI 2.203 0.189 0.661 1.054 0.954 4.132

IL 2.026 0.011 0.762 0.738 0.982 2.895

FIL 2.099 0.084 0.581 0.597 0.973 2.341

$E_{\max}$ CC 3.900 0.316 0.742 0.862 0.989 3.377

NRI 4.185 0.602 0.940 0.877 0.974 3.438

MI 3.841 0.258 0.736 0.803 0.974 3.147

IL 3.704 0.121 0.752 0.851 0.964 3.338

FIL 3.626 0.043 0.604 0.766 0.976 3.004

$E_{0}$ CC $- 2.396$ $- 0.199$ 0.628 0.781 0.973 3.061

NRI $- 2.618$ $- 0.421$ 0.716 0.740 0.991 2.901

MI $- 2.339$ $- 0.141$ 0.586 0.663 0.966 2.600

IL $- 2.199$ $-$ 0.002 0.642 0.749 0.928 2.936

FIL $- 2.177$ 0.002 0.512 0.686 0.954 2.691

250 log( ${ED}_{50}$ ) CC 2.058 0.043 0.528 0.521 0.968 2.042

NRI 2.486 0.471 0.698 0.501 0.849 1.964

MI 2.181 0.166 0.525 0.473 0.942 1.853

IL 1.964 -0.051 0.772 0.535 0.970 2.097

FIL 2.032 0.017 0.471 0.479 0.966 1.878

$E_{\max}$ CC 3.898 0.315 0.707 0.667 0.979 2.614

NRI 4.137 0.554 0.799 0.627 0.947 2.456

MI 3.858 0.275 0.691 0.587 0.948 2.301

IL 3.740 0.156 0.667 0.660 0.967 2.587

FIL 3.678 0.095 0.580 0.619 0.971 2.428

$E_{0}$ CC $- 2.451$ $- 0.254$ 0.653 0.623 0.980 2.443

NRI $- 2.668$ $- 0.470$ 0.738 0.591 0.996 2.319

MI $- 2.386$ $- 0.189$ 0.612 0.530 0.958 2.078

IL $- 2.303$ $- 0.106$ 0.630 0.614 0.955 2.406

FIL $- 2.260$ $-$ 0.062 0.531 0.569 0.966 2.232

350 log( ${ED}_{50}$ ) CC 2.087 0.072 0.448 0.436 0.969 1.708

NRI 2.514 0.499 0.673 0.423 0.783 1.657

MI 2.162 0.147 0.447 0.394 0.925 1.543

IL 2.005 $- 0.010$ 0.459 0.449 0.967 1.759

FIL 2.041 0.026 0.412 0.414 0.962 1.624

$E_{\max}$ CC 3.863 0.280 0.649 0.553 0.962 2.169

NRI 4.106 0.523 0.743 0.519 0.896 2.034

MI 3.816 0.232 0.603 0.485 0.937 1.901

IL 3.712 0.128 0.626 0.550 0.944 2.155

FIL 3.670 0.086 0.561 0.526 0.949 2.060

$E_{0}$ CC $- 2.436$ $- 0.238$ 0.603 0.520 0.984 2.037

NRI $- 2.653$ $- 0.455$ 0.696 0.493 0.958 1.933

MI $- 2.365$ $- 0.168$ 0.544 0.440 0.957 1.725

IL $- 2.295$ $- 0.098$ 0.583 0.514 0.942 2.017

FIL $- 2.265$ $- 0.067$ 0.509 0.487 0.955 1.909

450 log( ${ED}_{50}$ ) CC 2.077 0.062 0.397 0.378 0.947 1.482

NRI 2.499 0.484 0.628 0.367 0.750 1.438

MI 2.156 0.141 0.397 0.340 0.916 1.333

IL 1.998 -0.017 0.401 0.389 0.953 1.524

FIL 2.025 0.010 0.370 0.366 0.950 1.435

$E_{\max}$ CC 3.863 0.280 0.565 0.481 0.962 1.884

NRI 4.099 0.516 0.680 0.448 0.869 1.757

MI 3.794 0.210 0.528 0.421 0.916 1.650

IL 3.712 0.128 0.529 0.478 0.966 1.875

FIL 3.682 0.098 0.482 0.463 0.961 1.813

$E_{0}$ CC $- 2.428$ $- 0.231$ 0.532 0.452 0.974 1.772

NRI $- 2.644$ $- 0.446$ 0.634 0.428 0.925 1.680

MI $- 2.358$ $- 0.161$ 0.491 0.381 0.932 1.495

IL $- 2.287$ $- 0.090$ 0.504 0.448 0.956 1.756

FIL $- 2.266$ $- 0.068$ 0.453 0.430 0.953 1.685

Sample Size(N)	Parameter	Type	Estimate	MBE	RMSE	Est.SE	CP	Est.Length
150	log( ${ED}_{50}$ )	CC	2.104	0.089	0.679	0.698	0.979	2.736
		NRI	2.564	0.549	0.858	0.666	0.896	2.610
		MI	2.203	0.189	0.661	1.054	0.954	4.132
		IL	2.026	0.011	0.762	0.738	0.982	2.895
		FIL	2.099	0.084	0.581	0.597	0.973	2.341
	$E_{\max}$	CC	3.900	0.316	0.742	0.862	0.989	3.377
		NRI	4.185	0.602	0.940	0.877	0.974	3.438
		MI	3.841	0.258	0.736	0.803	0.974	3.147
		IL	3.704	0.121	0.752	0.851	0.964	3.338
		FIL	3.626	0.043	0.604	0.766	0.976	3.004
	$E_{0}$	CC	$- 2.396$	$- 0.199$	0.628	0.781	0.973	3.061
		NRI	$- 2.618$	$- 0.421$	0.716	0.740	0.991	2.901
		MI	$- 2.339$	$- 0.141$	0.586	0.663	0.966	2.600
		IL	$- 2.199$	$-$ 0.002	0.642	0.749	0.928	2.936
		FIL	$- 2.177$	0.002	0.512	0.686	0.954	2.691
250	log( ${ED}_{50}$ )	CC	2.058	0.043	0.528	0.521	0.968	2.042
		NRI	2.486	0.471	0.698	0.501	0.849	1.964
		MI	2.181	0.166	0.525	0.473	0.942	1.853
		IL	1.964	-0.051	0.772	0.535	0.970	2.097
		FIL	2.032	0.017	0.471	0.479	0.966	1.878
	$E_{\max}$	CC	3.898	0.315	0.707	0.667	0.979	2.614
		NRI	4.137	0.554	0.799	0.627	0.947	2.456
		MI	3.858	0.275	0.691	0.587	0.948	2.301
		IL	3.740	0.156	0.667	0.660	0.967	2.587
		FIL	3.678	0.095	0.580	0.619	0.971	2.428
	$E_{0}$	CC	$- 2.451$	$- 0.254$	0.653	0.623	0.980	2.443
		NRI	$- 2.668$	$- 0.470$	0.738	0.591	0.996	2.319
		MI	$- 2.386$	$- 0.189$	0.612	0.530	0.958	2.078
		IL	$- 2.303$	$- 0.106$	0.630	0.614	0.955	2.406
		FIL	$- 2.260$	$-$ 0.062	0.531	0.569	0.966	2.232
350	log( ${ED}_{50}$ )	CC	2.087	0.072	0.448	0.436	0.969	1.708
		NRI	2.514	0.499	0.673	0.423	0.783	1.657
		MI	2.162	0.147	0.447	0.394	0.925	1.543
		IL	2.005	$- 0.010$	0.459	0.449	0.967	1.759
		FIL	2.041	0.026	0.412	0.414	0.962	1.624
	$E_{\max}$	CC	3.863	0.280	0.649	0.553	0.962	2.169
		NRI	4.106	0.523	0.743	0.519	0.896	2.034
		MI	3.816	0.232	0.603	0.485	0.937	1.901
		IL	3.712	0.128	0.626	0.550	0.944	2.155
		FIL	3.670	0.086	0.561	0.526	0.949	2.060
	$E_{0}$	CC	$- 2.436$	$- 0.238$	0.603	0.520	0.984	2.037
		NRI	$- 2.653$	$- 0.455$	0.696	0.493	0.958	1.933
		MI	$- 2.365$	$- 0.168$	0.544	0.440	0.957	1.725
		IL	$- 2.295$	$- 0.098$	0.583	0.514	0.942	2.017
		FIL	$- 2.265$	$- 0.067$	0.509	0.487	0.955	1.909
450	log( ${ED}_{50}$ )	CC	2.077	0.062	0.397	0.378	0.947	1.482
		NRI	2.499	0.484	0.628	0.367	0.750	1.438
		MI	2.156	0.141	0.397	0.340	0.916	1.333
		IL	1.998	-0.017	0.401	0.389	0.953	1.524
		FIL	2.025	0.010	0.370	0.366	0.950	1.435
	$E_{\max}$	CC	3.863	0.280	0.565	0.481	0.962	1.884
		NRI	4.099	0.516	0.680	0.448	0.869	1.757
		MI	3.794	0.210	0.528	0.421	0.916	1.650
		IL	3.712	0.128	0.529	0.478	0.966	1.875
		FIL	3.682	0.098	0.482	0.463	0.961	1.813
	$E_{0}$	CC	$- 2.428$	$- 0.231$	0.532	0.452	0.974	1.772
		NRI	$- 2.644$	$- 0.446$	0.634	0.428	0.925	1.680
		MI	$- 2.358$	$- 0.161$	0.491	0.381	0.932	1.495
		IL	$- 2.287$	$- 0.090$	0.504	0.448	0.956	1.756
		FIL	$- 2.266$	$- 0.068$	0.453	0.430	0.953	1.685

When the sample size was n=150, failures or unstable estimates occurred in 1.8% (CC), 1.0% (NRI), 1.0% (MI), and 0.8% (IL) of replicates, whereas FIL always converged and yielded stable estimates. For larger sample sizes, all methods converged without issues. As shown in Table 1, NRI estimators exhibit poor MBE and RMSE compared to other methods across all scenarios, leading to lower CP despite having lower mean estimated standard errors. For CC estimators, although the RMSE occasionally performs better than IL, and they are nearly unbiased for ${ED}_{50}$ , the MBEs for $E_{\max}$ and $E_{0}$ are consistently large, resulting in a highly biased estimator. MI estimators show some reduction in bias for point estimates across all scenarios compared to CC and NRI; however, the bias remains substantial. Furthermore, the standard error for MI is underestimated, leading to a reduced nominal coverage for the confidence intervals. FIL always outperforms IL, especially in estimating $E_{0}$ and $E_{\max}$ in terms of MBE. Concerning RMSE and mean estimated standard errors, FIL estimators always achieve the lowest values across all scenarios. Due to the reduced standard errors, FIL also produces the narrowest confidence intervals. Moreover, the estimated standard errors for both IL and FIL tend to be slightly overestimated compared to the true simulated standard error of the parameters, with this overestimation being more pronounced in smaller sample sizes. This leads to conservative confidence intervals and an overestimation of nominal coverage. Alternatively, a profile likelihood-based paralleled parametric bootstrap confidence interval can be considered with relatively computational cost when sample sizes are small.

The FIL method for bias correction is particularly valuable when the sample size is small and separation is encountered in the missingness pattern. Additional simulations were performed where missingness was more correlated with the dose, leading to severe separation in the placebo treatment arm (results reported in Table S1). In this scenario, all methods perform worse due to the systematic bias introduced by the loss of extreme scorers (success cases) in the placebo group. However, FIL and IL still outperform CC and NRI, with FIL constantly achieving better results. Notably, IL significantly underestimates the nominal coverage, whereas FIL maintains a reasonable coverage probability. Additionally, the standard errors estimated by FIL are smaller than those from IL, demonstrating the effectiveness of the Firth-type method in mitigating the effects of separation.

Figure 1 illustrates the distribution of point estimates for different methods, along with estimates based on the full data without missingness, across varying sample sizes. The true value of $θ$ is indicated by a vertical line on the boxplot, and the mean estimates are shown as blue dots. The plot highlights that FIL effectively reduces bias due to small sample sizes and separation, particularly when $n = 150$ . Even with the full dataset, separation due to small sample size results in unstable MLE estimates, characterized by large variations and outliers. However, the application of Jeffreys’ prior modification in FIL, which introduces strong convexity in the likelihood function with respect to the parameter of interest, mitigates the constancy of the likelihood due to separation, leading to more reliable estimates. Additionally, while both IL and FIL perform similarly in terms of median estimates, the box lengths (indicating variability) for FIL are narrower than those for IL.

Figure 1.

Boxplots comparing the distribution of point estimates with true parameter $E_{0} = - 2.197$ , $E_{\max} = 3.584$ , and ${ED}_{50} = 2.015$ , based on 1000 replications and missing rate approximately 15%.

Further simulations were conducted with varying missing rates and a fixed sample size of $n = 350$ . The missing rates were set at 10%, 15%, 25%, and 30%, with different $α$ combinations used to produce the varying rates. As in real randomized control trials, the datasets with much higher missing rates are not suitable for analysis and require new samples to be included, scenarios with missing rates higher than 30% will not be considered. As shown in Table S2, the results are similar to those in 1. FIL outperforms IL, and both methods provide better estimates than CC, NRI, and MI. Additionally, FIL consistently achieves the lowest mean estimated standard errors, with coverage probabilities close to the desired 95% nominal level.

4. Real data analysis

This section presents an example of fitting a dose–response model using data from the TURANDOT study,²⁷ a Phase II randomized, double-blind, placebo-controlled clinical trial for ulcerative colitis in patients with moderate to severe disease. In this study, 357 patients were randomly assigned to either a placebo group or one of four active dose groups: 7.5 mg, 22.5 mg, 75 mg, and 225 mg. As reported by Vermeire et al., a non-monotone dose–response profile was observed, with lower efficacy in the highest dose group (225 mg). Among the 357 patients, 73 received placebo, 71 received 7.5 mg, 72 received 22.5 mg, 71 received 75 mg, and 70 received 225 mg. The primary endpoint was clinical remission at Week 12, which included several missing values believed to be nonignorable, as summarized in Table 2.

Table 2.
Sample size, the number of missing response cases, the number of remission cases, and baseline statistics for each dosage group in TURANDOT study.

Dose (mg) Placebo 7.5 22.5 75 225

Sample size 73 71 72 71 70

Missing response 6(8.2%) 8(11.2%) 1(1.4%) 3(4.2%) 6(8.6%)

Remission (Yes) 2(2.7%) 8(11.3%) 12(16.7%) 11(15.5%) 4(5.7%)

Previous TNF therapy (Yes) 42(57.5%) 41(57.7%) 41(57.7%) 41(57.5%) 40(57.1%)

Sex (Male) 44(60.3%) 39(54.9%) 46(64.3%) 37(52.1%) 42(60.0%)

IS (Yes) 15(20.5%) 23(32.4%) 23(32.9%) 21(26.0%) 20(30.0%)

SD (Yes) 31(42.5%) 38(53.5%) 38(52.8%) 36(50.7%) 36(51.4%)

ASA (Yes) 48(65.8%) 37(52.1%) 36(51.4%) 44(60.3%) 36(51.4%)

MCSBASE (Mean(SD)) 8.4(1.7) 8.7(1.7) 8.1(1.6) 8.4(1.9) 8.7(1.6)

CRPBASE (Mean(SD)) 1.106(0.5) 1.097(0.5) 1.149(0.5) 0.979(0.4) 0.892(0.4)

Age (Mean(SD)) 38.6(12.7) 41.3(14.7) 42.1(14.7) 37.5(12.4) 41.3(13.2)

Dose (mg)	Placebo	7.5	22.5	75	225
Sample size	73	71	72	71	70
Missing response	6(8.2%)	8(11.2%)	1(1.4%)	3(4.2%)	6(8.6%)
Remission (Yes)	2(2.7%)	8(11.3%)	12(16.7%)	11(15.5%)	4(5.7%)
Previous TNF therapy (Yes)	42(57.5%)	41(57.7%)	41(57.7%)	41(57.5%)	40(57.1%)
Sex (Male)	44(60.3%)	39(54.9%)	46(64.3%)	37(52.1%)	42(60.0%)
IS (Yes)	15(20.5%)	23(32.4%)	23(32.9%)	21(26.0%)	20(30.0%)
SD (Yes)	31(42.5%)	38(53.5%)	38(52.8%)	36(50.7%)	36(51.4%)
ASA (Yes)	48(65.8%)	37(52.1%)	36(51.4%)	44(60.3%)	36(51.4%)
MCSBASE (Mean(SD))	8.4(1.7)	8.7(1.7)	8.1(1.6)	8.4(1.9)	8.7(1.6)
CRPBASE (Mean(SD))	1.106(0.5)	1.097(0.5)	1.149(0.5)	0.979(0.4)	0.892(0.4)
Age (Mean(SD))	38.6(12.7)	41.3(14.7)	42.1(14.7)	37.5(12.4)	41.3(13.2)

We fitted the Emax model, assuming all missing values in remission were nonignorable. For modeling the missingness indicator, we performed model selection based on the Akaike Information Criterion (AIC) with a lot of covariates and their interactions, and the final model included the following covariates: remission response (y), dose, Mayo score at baseline (MCSBASE), C Reactive Protein score at baseline (CRPBASE), age, sex, immune suppressant status (IS), steroid use history (SD), and Acetylsalicylic acid use history (ASA). We compared our proposed method to complete case (CC) analysis, NRI, and MI. Note that for MI, we used the imputation model identical to the one employed in IL and FIL, excluding the response y. The results of the dose–response relationship using different methods are shown in Table 3.

Table 3.

Analysis result of TURANDOT data with different missing data handling methods and the proposed method.

Parameter	Method	Estimate	StdErr	95% CI
log( ${ED}_{50}$ )	CC	0.480	1.856	$(- 3.159, 4.119)$
	NRI	0.756	1.484	$(- 2.153, 3.664)$
	MI	0.462	4.611	$(- 8.576, 9.500)$
	IL	$- 1.775$	15.065	$(- 31.303, 27.752)$
	FIL	1.030	0.907	$(- 0.747, 2.808)$
$E_{\max}$	CC	1.938	0.788	(0.394, 3.481)
	NRI	2.017	0.788	(0.472, 3.561)
	MI	1.851	0.741	(0.400, 3.303)
	IL	1.563	0.748	(0.098, 3.029)
	FIL	1.836	0.721	(0.423, 3.249)
$E_{0}$	CC	$- 3.484$	0.718	$(- 4.890, - 2.077)$
	NRI	$- 3.576$	0.716	$(- 4.980, - 2.172)$
	MI	$- 3.400$	0.667	$(- 4.707, - 2.092)$
	IL	$- 3.080$	0.685	$(- 4.423, - 1.738)$
	FIL	$- 3.285$	0.642	$(- 4.543, - 2.027)$

We observe that the IL method yields an unstable estimate for log( ${ED}_{50}$ ), with a very wide 95% confidence interval. This instability is likely due to separation issues in the logistic regression for the missingness indicator. Similarly, the MI method also produces a high standard error for estimating log( ${ED}_{50}$ ). In contrast, the FIL method provides a more stable estimate due to the use of penalized maximum likelihood estimation. Notably, the standard errors of the estimated parameters using FIL are consistently smaller, except for the IL method, which is affected by separation. As expected, the NRI method estimates the smallest value for $E_{0}$ , resulting in the smallest maximum achievable effect ( $E_{\max} + E_{0}$ ) among all methods.

To further evaluate the performance of these methods, we estimated the probabilities of remission using the coefficient estimates of the fitted Emax model. To quantify uncertainty of probabilities, we use a parametric bootstrap with $B = 5000$ replications based on the asymptotic normality of coefficient estimates: we sample the parameter vector $θ$ from its estimated joint asymptotic distribution, transform each draw through the Emax model to obtain estimated probability for dose $d$ : ${\hat{p}}^{(b)} (d)$ , and use the 2.5th and 97.5th percentiles across bootstrap draws at each dose $d$ as boundaries of 95% confidence interval. The results are as shown in Figure 2. We observe that NRI provides the lowest estimated probabilities across all dose groups due to its imputation strategy, but with relatively narrow confidence intervals. FIL, on the other hand, offers stable estimations across all dose groups with the smallest variances. Due to the unstable estimation of ${ED}_{50}$ , the IL method estimates the probability that nearly reaches the maximum effect at 7.5 mg, accompanied by comparably large and asymmetric confidence intervals across all dose groups. It is also noteworthy that the MI method exhibits the largest variance in estimations, despite its mean estimates being similar to those of FIL. This may be attributed to a violation of the missing mechanism assumption, which leads to biased estimation of parameter covariance under the non-linear model.

Figure 2.

Estimated dose response remission probabilities based on different methods with their bootstrapped 95% confidence intervals.

Table 4 presents the parameter estimates from the logistic regression of the missingness indicator, including their standard errors, Z-values, and p-values. The response variable clinical remission, coded as y, is significant at the 5% significance level while fitting the model in both the Ibrahim-Lipsitz (IL) and Firth-adjusted IL (FIL) methods, indicating that the probability of a response being missing depends on the response itself. This result provides supportive evidence that the missing data mechanism in this dataset may be nonignorable, while it is impossible to distinguish from MAR with the observed data only.²⁸

Table 4.

Estimates of missing data model for TURANDOT study.

Method	Variable	Estimate	StdErr	Z-value	p-value
IL	intercept	$- 4.006$	2.233	$- 1.794$	0.073
	y	2.532	0.595	4.254	<0.001
	Dose	$- 0.027$	0.016	$- 1.726$	0.084
	MCSBASE	0.222	0.202	1.098	0.272
	CRPBASE	0.260	0.127	2.042	0.041
	AGE	-0.042	0.027	$- 1.550$	0.121
	Sex (Male)	1.083	0.680	1.592	0.113
	IS (Yes)	$- 0.855$	0.840	$- 1.018$	0.309
	SD (Yes)	0.481	0.671	$- 0.915$	0.360
	ASA (Yes)	$- 0.614$	0.676	0.710	0.478
FIL	intercept	-2.510	1.824	$- 1.375$	0.169
	y	1.063	0.493	2.156	0.031
	Dose	-0.015	0.011	$- 1.396$	0.163
	MCSBASE	0.123	0.169	0.726	0.467
	CRPBASE	0.223	0.113	1.975	0.048
	AGE	$- 0.038$	0.022	-1.735	0.083
	Sex (Male)	0.967	0.551	1.756	0.079
	IS (Yes)	$- 0.663$	0.657	$- 1.008$	0.313
	SD (Yes)	0.250	0.550	0.455	0.650
	ASA (Yes)	$- 0.642$	0.550	$- 1.168$	0.243

5. Conclusion and discussion

In this article, we addressed the challenge of estimating the coefficients for the binary Emax model in the presence of missing responses under a nonignorable missing data mechanism. Our simulation studies demonstrate that the proposed Firth-type corrected weighted EM procedure of Ibrahim and Lipsitz (FIL) outperforms commonly used missing data handling strategies such as NRI and MI. Additionally, when fitting the binary Emax model with small or medium sample sizes, maximum likelihood-based approaches often face convergence issues due to complete separation or produce unstable estimates with significant bias and variation due to quasi-separation. In such scenarios, the FIL method offers a robust and reliable solution, effectively addressing the complexities and challenges posed by the data.

For both the IL and FIL methods, it is crucial to select an appropriate model for describing the likelihood of the missingness indicator, $f (r_{i} ∣ z_{i}, α)$ . Although it is known that identifiability of parameters is problematic without making additional untestable assumptions, a practical way of model selection in the range of potential assumed models still deserves discucssion. Standard model selection techniques, such as backward selection using AIC or likelihood ratio tests (as discussed by Ibrahim and Lipsitz²⁵), could be used in conjunction with considerations of scientific or clinical relevance. These techniques should be applied with the joint likelihood $f (θ, α ∣ z_{i}, r_{i})$ while keeping $f (y_{i} ∣ {Dose}_{i}, θ)$ fixed. Additionally, commonly used penalized likelihood variable selection methods, such as Lasso and Elastic-net regression, could be considered. In cases where the assumed missingness mechanism is well separated by regression models, a Firth-corrected regression model, as implemented in the FIL method, may be appropriate. Our simulations indicate that the estimation of $α$ can significantly impact the estimation of $θ$ under the Emax model, particularly when separation occurs in predicting the missingness indicator. Therefore, it is advisable to examine the estimation results of $α$ and consider adjustments if instability is observed.

The rationale for choosing an EM approach rather than multiple imputation (MI) to address missingness may warrant further explanation. As discussed in the introduction, MI is primarily designed for the MAR mechanism, which does not apply when the missing data are nonignorable. While recent literature, such as Im & Kim²⁹ and Galimard et al.,¹³ has proposed methods to handle nonignorable missingness, the selection of an appropriate imputation model remains a challenge and is often determined through sensitivity analysis. In contrast, with the IL and FIL methods, model selection can be performed using standard techniques, making these methods more straightforward to implement and interpret.

In dose–response analysis, the Multiple Comparison Procedure-Modeling (MCP-Mod) approach, developed by Bretz, Pinheiro, and Branson,^30,31 combines hypothesis testing and modeling with Type I error control. MCP-Mod uses AIC to select the best model for fitting the dose–response relationship, making it a natural extension to apply our proposed methods within the MCP-Mod framework, given its likelihood-based foundation. Indeed, Diniz et al. have recently developed a Firth-type MCP-Mod for Weibull regression with time-to-event data in small sample sizes.³² Further research could explore the implementation of these methods for binary or count dose–response models. In fact, MCP-mod contains a four-parameter Emax model as potential, thus extending to MCP-mod can capture situations when the three-parameter Emax model is misspecified. Additionally, penalization methods could be applied to the potential models within MCP-Mod to control the risk of separation, offering another avenue for future investigation.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802251403356 - Supplemental material for Robust Emax model fitting: Addressing nonignorable missing binary outcome in dose–response analysis

Supplemental material, sj-pdf-1-smm-10.1177_09622802251403356 for Robust Emax model fitting: Addressing nonignorable missing binary outcome in dose–response analysis by Jiangshan Zhang, Vivek Pradhan and Yuxi Zhao in Statistical Methods in Medical Research

Footnotes

Acknowledgements

The authors are grateful to an Associate Editor and referees, whose insightful comments have helped improve the work and presentation.

ORCID iDs

Jiangshan Zhang

Vivek Pradhan

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Supplemental material

Supplemental material for this article is available online.

References

Firth

. Bias reduction of maximum likelihood estimates. Biometrika 1993; 80: 27–38.

Maiti

Pradhan

. A comparative study of the bias corrected estimates in logistic regression. Stat Methods Med Res 2008; 17: 621–634.

Little

Rubin

. Statistical analysis with missing data. Nashville, TN: John Wiley & Sons, 2025.

Little

RJA

. Models for nonresponse in sample surveys. J Am Stat Assoc 1982; 77: 237–250.

Thomas

Roy

. Analysis of clinical dose—response in small-molecule drug development: 2009–2014. Stat Biopharm Res 2017; 9: 137–146.

Kirby

Brain

Jones

. Fitting Emax models to clinical trial dose–response data. Pharm Stat 2011; 10: 143–149.

Banerjee

Jin

, et al. Clinical dose–response for a broad set of biological products: a model-based meta-analysis. Stat Methods Med Res 2017; 27: 2694–2721.

O’Neill

Temple

. The prevention and treatment of missing data in clinical trials: an FDA perspective on the importance of dealing with it. Clin Pharmacol Therapeut 2012; 91: 550–554.

Papp

Fonjallaz

Casset-Semanaz

, et al. Analytical approaches to reporting long-term clinical trial data. Curr Med Res Opin 2008; 24: 2001–2008.

10.

Rubin

. Multiple imputation for nonresponse in surveys. Nashville, TN: John Wiley & Sons, 1987.

11.

Giusti

Little

RJA

. An analysis of nonignorable nonresponse to income in a survey with a rotating panel design. J Off Stat 2011; 27: 211–229.

12.

Galimard

Chevret

Protopopescu

, et al. A multiple imputation approach for MNAR mechanisms compatible with heckman’s model. Stat Med 2016; 35: 2907–2920.

13.

Galimard

Chevret

Curis

, et al. Heckman imputation models for binary or continuous MNAR outcomes and MAR predictors. BMC Med Res Methodol 2018; 18: 90.

14.

van Buuren

Boshuizen

Knook

. Multiple imputation of missing blood pressure covariates in survival analysis. Stat Med 1999; 18: 681–694.

15.

Carpenter

Kenward

White

. Sensitivity analysis after multiple imputation under missing at random: a weighting approach. Stat Methods Med Res 2007; 16: 259–275.

16.

Albert

Anderson

. On the existence of maximum likelihood estimates in logistic regression models. Biometrika 1984; 71: 1–10.

17.

Altman

Gill

McDonald

. Numerical issues in statistical computing for the social scientist. Nashville, TN: John Wiley & Sons, 2003.

18.

Thomas

Smith

. clindr: simulation and analysis tools for clinical dose response modeling. CRAN: contributed Packages 2017. DOI: 10.32614/cran.package.clindr.

19.

Bornkamp

Pinheiro

Bretz

, et al. Dosefinding: planning and analyzing dose finding experiments. CRAN: contributed Packages 2010. DOI: 10.32614/cran.package.dosefinding.

20.

Chen

Fries

Leonov

. Longitudinal model for a dose-finding study for a rare disease treatment. Stat Papers 2023; 64: 1343–1360.

21.

Heinze

Schemper

. A solution to the problem of separation in logistic regression. Stat Med 2002; 21: 2409–2419.

22.

Heinze

. A comparative investigation of methods for logistic regression with separated or nearly separated data. Stat Med 2006; 25: 4216–4226.

23.

Maiti

Pradhan

. Bias reduction and a solution for separation of logistic regression with missing covariates. Biometrics 2009; 65: 1262–1269.

24.

Maity

Pradhan

Das

. Bias reduction in logistic regression with missing responses when the missing data mechanism is nonignorable. Am Stat 2019; 73: 340–349.

25.

Ibrahim

Lipsitz

. Parameter estimation from incomplete data in binomial regression when the missing data mechanism is nonignorable. Biometrics 1996; 52: 1071–1078.

26.

Louis

. Finding the observed information matrix when using the em algorithm. J R Stat Soc Ser B: Stat Methodol 1982; 44: 226–233.

27.

Vermeire

Sandborn

Danese

, et al. Anti-madcam antibody (pf-00547659) for ulcerative colitis (turandot): a phase 2, randomised, double-blind, placebo-controlled trial. The Lancet 2017; 390: 135–144.

28.

Molenberghs

Beunckens

Sotto

, et al. Every missingness not at random model has a missingness at random counterpart with equal fit. J R Stat Soc Ser B: Stat Methodol 2008; 70: 371–388.

29.

Kim

. Multiple imputation for nonignorable missing data. J Korean Stat Soc 2017; 46: 583–592.

30.

Bretz

Pinheiro

Branson

. Combining multiple comparisons and modeling techniques in dose–response studies. Biometrics 2005; 61: 738–748.

31.

Pinheiro

Bornkamp

Glimm

, et al. Model-based dose finding under model uncertainty using general parametric models. Stat Med 2014; 33: 1646–1661.

32.

Diniz

Gallardo

Magalhães

. Improved inference for mcp-mod approach using time-to-event endpoints with small sample sizes. Pharm Stat 2023; 22: 760–772.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.14 MB