Improved inference for MCP-Mod approach using time-to-event endpoints with small sample sizes

Adequate designs for early-phase trials is essential to a successful clinical drug development. Traditionally, investigators evaluate safety in phase I trials, proof-of-concept (PoC) in phase IIa trials, and efficacy in phase IIb trials. When drugs are promising in these early stages, phase III trials are designed accruing a large number of patients to provide a definitive evidence of efficacy. Nonetheless, Jardim et al. [1] showed that 20% of 80 cancer drug programs submitted between 2009 to 2015 to the Food and Drug Administration (FDA) did not have any data for phase II trials such that 31% obtained FDA approval; 46% had positive PoC with of 76% FDA-approval; and 34% had negative PoC with only 15% of FDA-approval. Therefore, attaining proof-of-concept is a predictor to a successful phase III trial.

A classic PoC is designed with a highest dose allowable based on Phase I clinical trial results to compare with placebo in a two-arm design or a historical threshold in a one-arm design. Once proof of concept is shown, dose range studies to identify a minimum effective dose (MED) or target dose (TD) are conducted following two possible approaches: multiple comparisons (MCP) and modeling (Mod). The MCP approach corresponds to evaluate contrasts among doses while preserving the family-wise error rate (FWER), which is robust to distribution assumptions but restricted to the set of doses under investigation; the Mod approach assumes a dose-response relationship allowing to estimate the response for a dose even if that dose is not under investigation but it heavily depends on choosing the appropriate functional form.

Bretz et al. [2] proposed a framework named Multiple Comparison Procedures with Modeling Techniques (MCP-Mod) unifying phase IIa and IIb trials into a seamless design while taking advantage of both traditional approaches for normally distributed endpoints. Later, Pinheiro et al. [3] extended this methodology to general parametric models using generalized least squares estimation allowing statisticians to consider more complex designs and other types of endpoints such as binary and time to event. Regulatory agencies have stated their approval of the MCP-Mod framework as an adequate and efficient methodology for design and analysis of phase II dose-finding studies that will guide dose selection for phase III trials.

The MCP-Mod framework is implemented in two steps: (i) MCP-step consists of a trend test to assess the presence of a dose response signal among a set of pre-specified candidate models while preserving FWER; (ii) Mod-step corresponds to estimate dose-response curves in order to identify the optimal dose that achieves a desired level of response in comparison to placebo among the models that were selected in the MCP-step. Therefore, the properties of the trend test defined as a Wald statistic used in the MCP-step and the maximum likelihood (ML) estimators with their covariance matrix used in the Mod-step are essential to the successful implementation of MCP-Mod approach. However, both steps rely on the the asymptotic properties of the ML estimators, which are only valid for large sample sizes.

In cancer mouse studies, time-to-death is an endpoint that could be used to guide the identification of the minimum effective dose (MED) with large expected effect sizes and small sample sizes. Nonetheless, the application of the MCP-Mod framework is limited due the underlying asymptotic assumptions. In this context, parametric survival models such as the censored Weibull regression model (WRM) was showed to be useful given that provides clinical meaningful interpretation based on event time ratios or hazard ratios[4]. Moreover, the asymptotic properties of the ML estimators can be studied and refined for small sample sizes. Recently, Magalhães et al. [5] obtained the skewness coefficient of the distribution of the maximum likelihood estimators for WRM, and Magalhães et al. [6] derived improved test statistics for LR, score and gradient tests but not for the Wald statistic.

Our main goal is to derive improved ML estimators for the regression parameters and its second-order covariance matrix for WRM, then use them as input to the generalized least squares procedure proposed by Pinheiro et al. [3] yielding a type I error probability closer to the nominal value when testing proof-of-activity and more accurate MED estimates. Moreover, we particularize the results from Cox and Snell [7] and Magalhães et al. [8] that are very general to WRM and they can be used in a much broader context than the MCP-Mod framework.

The remaining paper is organized as follows: in Section 2, we revisit the Weibull distribution and its properties; in Section 3, we review the main concepts of the MCP-Mod framework; in Section 4, we introduce the improved estimators; in Section 5, we conducted a simulation study to evalute the use of improved estimators in the MCP-Mod framework; finally, in Section 6, we presented some concluding remarks.

The Weibull distribution [9] is commonly used to analysis of time-to-event or lifetime data and a continuous random variable $T$ is called Weibull, if its probability density function (pdf) is

In lifetime data, there is the censoring restriction, i.e, if $T_{1},\ldots,T_{n}$ are a random sample from (1), instead of $T_{i}$, we observe, under right censoring, $t_{i}=\min(T_{i},L_{i})$, where $L_{i}$ is the censoring time, independent of $T_{i}$ and $\delta_{i}=1$, if $T_{i}\leq L_{i}$ or $\delta_{i}=0$, otherwise, $i=1,\ldots,n$. Here, we consider an hybrid censoring scheme, where the study is finalized when a pre-fixed number $r\leq n$ out of $n$ observations have failed, as well as when a prefixed time, say $L_{1}=\cdots=L_{n}=L$, has been reached. The type I censoring is a particular case for $r=n$ and the type II censoring appears when $L_{1},\ldots,L_{n}=+\infty$. Additionally, we add the non-informative censoring assumption, i.e., the random variables $L_{i}$ does not depend on $\lambda$. Usually, the regression modeling considers the distribution of $Y_{i}=\log(T_{i})$ instead of $T_{i}$, which is an accelerated lifetime model form, see Kalbfleisch and Prentice[10]. The distribution of $Y_{i}$ is of the extreme value form with pdf given by

The total score function and the total Fisher information matrix for $\beta$ are, respectively, ${\mbox{\boldmath{$U$}}}_{{\bm{\beta}}}=\sigma^{-1}{\mbox{\boldmath{$X$}}}^{\top}{\mbox{\boldmath{$W$}}}^{1/2}{\mbox{\boldmath{$v$}}}$ and ${\bm{K}}_{{\bm{\beta}}{\bm{\beta}}}=\sigma^{-2}{\mbox{\boldmath{$X$}}}^{\top}{\mbox{\boldmath{$W$}}}{\mbox{\boldmath{$X$}}}$, where ${\mbox{\boldmath{$X$}}}=({\mbox{\boldmath{$x$}}}_{1},\ldots,{\mbox{\boldmath{$x$}}}_{n})^{\top}$, the model matrix, assuming rank$({\mbox{\boldmath{$X$}}})=p$, ${\mbox{\boldmath{$W$}}}=$ diag$(w_{1},\ldots,w_{n})$, $w_{i}=\mathds{E}\left[\exp\left(\frac{y_{i}-\mu_{i}}{\sigma}\right)\right]$ and ${\mbox{\boldmath{$v$}}}=(v_{1},\ldots,v_{n})^{\top}$, $v_{i}=\left\{-\delta_{i}+\exp\left(\frac{y_{i}-\mu_{i}}{\sigma}\right)\right\}w_{i}^{-1/2}$. It can observed that the value of $w_{i}$ depends on the mechanism of censoring. That means $w_{i}=q\times\left\{1-\exp\left[-L_{i}^{1/\sigma}\exp(-\mu_{i}/\sigma)\right]\right\}+(1-q)\times\left(r/n\right)$, where $q=\mathds{P}\left(W_{(r)}\leq\log L_{i}\right)$ and $W_{(r)}$ denotes the $r$th order statistic from $W_{1},\ldots,W_{n}$. Note that $q=1$ and $q=0$ for types I and II censoring, respectively, as showed in Magalhães et al. [5]. The maximum likelihood estimator (MLE) of $\beta$, $\widehat{{\mbox{\boldmath{$\beta$}}}}$, is the solution of ${\mbox{\boldmath{$U$}}}_{{\bm{\beta}}}={\bf 0}$. The $\widehat{{\mbox{\boldmath{$\beta$}}}}$ can not be expressed in closed-form. It is typically obtained by numerically maximizing the log-likelihood function using a Newton or quasi-Newton nonlinear optimization algorithm. Under mild regularity conditions and in large samples,

We briefly summarize the two-stage procedure discussed by Pinheiro et al. [3] following the same notation. We consider that time-to-event responses $t_{ij}$ for doses $x_{i}$ given to $j$th subject for $i=0,\ldots,p$ and $j=1,\ldots,n_{i}$ can be described by the Weibull distribution with scale parameters $\lambda_{i}$ and shape parameter $\sigma$ defined in (1) and (3). Then, we consider the parameter $\mu$ as our response for the dose-response model such that it could be defined as the median survival time (2) for $\rho=0.5$ or, alternatively, $\log(\lambda)$.

Initially, a set of candidate dose-response models $\mu_{i}=f(x_{i},\mbox{\boldmath{$\theta$}})$ is considered such that each model can be rewritten as function of standardized model as below

In our collaborative work, investigators wanted to establish a dose-response relationship between a new inhibitor agent for pancreatic cancer in combination with a given dose of gemcitabine in mouse models. Based on preliminary data, a survival median time of 4 months was estimated in control-treated KPC mice model such that previous studies showed no survival benefit with only gemcitabine [13]. Assuming a Weibull distribution with $\sigma=0.5$, we calculated the placebo effect equal to 1.57 and the maximum effect of 2.26 considering a hazard ratio of 4. Investigators were interested in the minimum effective dose yielding a minimum hazard ratio of 2 corresponding to $\Delta=0.693$.

Five scenarios (constant, emax, exponential, logistic and beta model) presented in Table 1 were studied. True parameters were defined based on the aforementioned placebo and maximum effects, and the percent of maximum effect that is achieved at given dose as discussed in Bornkamp et al.[14]. For each scenario, the doses 0, 5, 25, 50 and 100 (mg/kg) were considered such that true MED was calculated as continuous dose given in Table 1. Five candidate models (linear, emax, exponential, logistic and beta model) were considered with guess estimates defined as the true parameters to calculate the optimal contrasts in (8). For each scenario, the null hypothesis of non-signal of the new inhibitor agent was tested at 5% significance level, and we used Akaike Information Criteria (AIC) to choose the model to estimate MED when more than one model rejected the non-signal hypothesis. Furthermore, we assumed that the target dose is estimated as a continuous dose for any value within the range of the dose grid. In case, the target dose is estimated outside of the dose grid, then the closest dose is selected.

For both steps of the MCP-Mod framework, the five strategies (MLE, MLE2, BCE, BCE2 and Firth) were evaluated based on a Monte Carlo simulation study with 100,000 replicates for sample sizes ranging from 5 to 100 mice per dose and censoring rate of 10%, 25% and 50%. The following operating characteristics were studied: (i) convergence rate of GLS algorithm in the MCP-Mod General framework; (ii) Probability of incorrectly detecting a dose-response signal under the scenario with a constant dose response curve, i.e., type I error; (iii) Probability of correctly detecting a dose-response signal under scenarios with a non-constant dose-response curve (Emax, Exponential, Logistic, Beta), i.e., power; (iv) Probability of selecting the true model given that there is a signal; (v) Bias of MED estimate and (vi) RMSE of MED estimate.

We have derived improved inferences based on the Wald statistic for WRM particularizing general results from Cox and Snell [7] and Magalhães et al. [8], which complements previous results for improved inference based on likelihood ratio, Rao score and gradient statistics discussed in Magalhães and Gallardo [6]. Few authors have presented improved inference for survival models with small sample sizes under the classical approach: Cordeiro and Colosimo [15, 16] and Medeiros [17] derived those statistics, respectively, for censored exponential regression models (ERM), a particular case of censored WRM. Also for ERM, Lemonte [18] presented the second-order covariance matrix of the MLE.

We have also proposed strategies based on bias-corrected (BCE, BCE2 and Firth) and second-order covariance matrices (MLE2, BCE2) as input for the general MCP-Mod framework introduced by Pinheiros et al.[3], which addresses the issue of relying on the asymptotic properties of MLE, which might not be valid for small sample sizes. To the best of our knowledge, this work is the first attempt to apply refined estimators for small sample sizes in the general MCP-Mod framework. Two simulations studies were performed to study the properties of refined estimators in an usual context of regression models and relevant operating characteristics in the MCP-Mod framework.

In the simulation study presented for general censored Weibull regression model in the supplementary material, we showed numerical evidences that BCE and Firth estimator have lower bias and RMSE than MLE; second-order covariance matrices evaluated at MLE and BCE are closer to their respective empirical covariance matrices in comparison to the first-order covariance matrices; and Wald statistics derived from the combination between bias-corrected estimators and second-order covariance matrices yielded type I error probability closer to the nominal value than the standard Wald statistic with no loss of power.

In the simulation study for the MCP-Mod framework, we have found that refined estimators and second-order covariance matrices approximate type I error probability to its nominal value in the MCP-step, while there are negligible differences between MLE and refined estimators in the probability of correctly detecting the dose-response signal, probability correctly selecting the dose-response model, bias and RMSE when censoring rates are up to 25%. For censoring rate of 50%, we found convergence issues for the corrected estimators and second-order covariance matrices and poorer performance in the assessed operating characteristics.

In conclusion, we recommend the use of Firth as a strategy using refined estimators for small sample sizes in the MCP-Mod framework. In the context of basic science with limited sample sizes and large effect sizes, we do not expect large censoring rates in mouse experiments. In human trials, smaller effect sizes are pursued such that larger sample sizes are required. In this case, proposed strategies are not needed because all estimators present comparable performance for large sample sizes. Nonetheless, we showed that type I error probability is still inflated even for sample sizes of 25 subjects per dose, therefore, the use of proposed strategies would avoid to dedicate further efforts on non-promising drugs.

We hope that refined estimators allow statisticians to implement the MCP-Mod framework with small sample sizes accelerating the pre-clinical and clinical drug development process. An R-package, MCPModBC [19], is available on CRAN to fit the Weibull model with refined estimators and perform simulations for power considerations. Similar ideas can be applied to other distributions such as Binomial, Negative Binomial and Poisson, and they are currently under investigation. Furthermore, the impact of model misspecification requires further study.

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.