Targeted design for adaptive clinical trials via semiparametric model

1 Introduction

Precision medicine, in which the treatment regimen is assigned to each patient according to that patient’s personal profile, is revolutionizing health care. The targeted phase II clinical trial design is the initial assessment of efficacy or screening for promising treatments. Such clinical trials use the patients’ personal profiles, such as biomarkers, to identify and help patients likely to benefit significantly from the treatments. In contrast to conventional clinical trials, targeted designs are aimed at identifying sub-population of patients with specific profiles for which the treatments are particularly effective. Such a design can improve the efficiency of a randomized clinical trial [1], [2]. Multi-stage design, also called group sequential design, is commonly used in phase II clinical trials to monitor the trial process, and the trial can be stopped early if extreme response effect or the lack of it thereof is observed at the intermediate stage [3], [4], [5]. Zhou et al. [6] and Lee et al. [7] proposed Bayesian methods for such design; Tang and Zhou [8] considered using biomarkers and optimal allocation for such design; and Gao [9] proposed multi-stage adaptive biomarker-directed clinical design.

Here, the goal is to find possible favorable subgroups for each treatment based on the subjects’ covariate profiles, and the two subgroups for treatments 1 and 2 may not be the same. This problem has some connection but differs from subgroup analysis, in which the goal is to identify a subgroup on which one treatment has significantly better effect than the other [10], [11], [12]. Shen and He [13] proposed a logistic-normal mixture model for subgroup analysis.

Since the patients’ personal profiles are used along with their responses in the assessment of treatment outcomes, the parametric linear regression mixture model will be used. Existing methods often assume some known model, for example, a Gaussian linear model, to estimate the parameters. In practice, the specified model is subjective and misspecified, so that the estimated parameters may not be consistent or may have an inflated error. Various methods have been considered to relax this assumption. Some authors considered different or more general parametric error distributions. Hanson and Johnson [14] considered the error from a mixture of Polya trees. Atsedeweyn and Rao [15] considered the error from a generalized new symmetric distribution. The most robust method is to make no assumption on the error distribution. Yuan and De Gooijer [16] considered such a scenario, used a kernel estimator to estimate the error density and construct the estimated likelihood, and then estimated the regression parameters by the MLE under the estimated likelihood. Yao and Zhao [17] considered a two-step method: they first estimate the errors by the difference between response and the covariates at the LSE value, estimate the error density by kernel estimator, and then estimate the regression parameters by the MLE under the estimated likelihood. However, the kernel or other density estimation methods, such as spline, basis method or penalized smoothing method, still have subjective inputs, such as the choices of kernel, bandwidth, basis, number of bases, penalty term, and smoothing parameter. Determining these quantities in many cases is not easy. The nonparametric maximum likelihood estimate (NPMLE) is appealing because it is completely objective, but it does not exist in the general case. Under certain shape constraint(s), such as log-concave, the NPMLE is unique and well-studied [18], [19], and the asymptotic distributions of the NPMLE and its mode are studied by [20]. Application of this method can be found in [21] and others.

The mixture regression model is increasingly utilized in practice, for example, in clustering, classification, personalized precision medicine, subgroup analysis, and survival analysis. To address robustness in a phase II clinical trial, we propose and study a semiparametric regression mixture model in which the mixing proportions are specified according to the subjects’ profiles, and for robustness, each sub-group distribution is assumed to be unimodal with mode at zero. The nonparametric maximum likelihood estimate (NPMLE), in this case as a form of isotonic regression, is used to estimate the error densities along with the MLE of the regression parameters. The asymptotic properties of the estimators are investigated.

The rest of this article is organized as follows. In Section 2, we introduce our methodology. Section 3 presents the asymptotic properties of the proposed method. A real data set is used for illustration in Section 4, followed by simulation studies in Section 5. Concluding remarks are given in Section 6.

2 The method

In the targeted clinical trial design, a patient population with a specific disease etymology is identified by some covariate profile, such as biomarker(s). Before the trial, patients are tested on the specific profile, such that often only profile-positive patients enter the trial. This design is also referred as enrichment design in the sense that a selected population is further studied in the trial. However, there is debate as to whether test-negative patients should also enter the trial and be investigated along with the test-positive patients. In this case, we just record the test result for each patient and register it as one component in the profile. In such phase II trial, each patient is randomized into several treatment arms; here, without loss of generality, we consider the case of two arms. In each treatment group, there may be a sub-population of patients whom respond to the treatment significantly better than the others. Below, we describe the proposed model for the problem and its estimation. Then, we describe the specifics of the multi-stage design trial.

3 Asymptotic results

Recall the definition of g r ( ɛ r | θ , f , x r ) given in (3). Let s = (y, x, t) and γ = P(t = 1). Omitting the density for x, the density-mass function for s is

Denote ℓ ( s | θ , f ) = log   g ( s | θ , f ) , ‖ f − f 0 ‖ = ∫ ​ | f ( t ) − f 0 ( t ) | d t , and let ‖ θ ˆ n − θ 0 ‖ be the Euclidean distance between θ ˆ n and θ 0 . We need the following regularity conditions:

1. Θ is bounded.

2. The support of x is bounded.

3. sup f ∈ ℱ f ( 0 ) < ∞ .

4. f 0 ∈ ℱ and f 0 ( ⋅ ) is continuous.

5. ℓ ( s | θ , f ) is Lipschitz in ( θ , f ) .

4 Application to real data

We analyze the real data AIDS Clinical Trials Group (ACTG) 175 with the proposed method. The drugs used in the trial were provided by three manufacturers Glaxo-Wellcome Company (AZT), Bristol-Myers Squibb Inc (ddI), and Hoffmann-La Roche Inc (ddC). The trial was conducted at ACTG’s 43 sites with both adult and pediatric patients and at nine sites of the National Hemophilia Foundation. Participants were enrolled in the study between December 1991 and October 1992 and received treatment through December 1994. Follow-up and final evaluations of participants took place between December 1994 and February 1995.

The purpose of the trial was to evaluate whether treatment of HIV infection with one drug (monotherapy) was the same, better than, or worse than treatment with two drugs (combination therapy) in patients with CD4+ T cells between 200 and 500/mm³. The results of original the ACTG 175 together with the results from earlier studies demonstrate that antiretroviral therapy is beneficial to HIV-infected people who have less than 500 CD4+ T cells/mm³.

We analyze these data using the proposed method on treatment 1 (with AZT) and treatment 4 (with ddI), assuming they construct the first stage of a multi-stage adaptive trial. There were 532 and 561 patients in treatments 1 and 4, respectively. In this analysis, the response variable is the CD4 counts after 20 weeks of the corresponding treatment, square root transformed, and the covariates are baseline age (re-scaled as 10 years incremental) and baseline CD4 counts, also square root transformed.

As a comparison, we also analyze the data using a parametric normal model (NM). For the semiparametric method, we fit the data assuming unimodality (at zero) and symmetric f ( ⋅ ) (SP1) or asymmetric (SP2). In the EM algorithm, we first fit a linear regression on the covariates to obtain the starting values β ( 0 ) . Initial membership (for the treatment-favorable subgroup) classification circles through three magnitudes — 1/2, 1/4 and 1/8 of the standard deviation of the residuals — of the choice of μ r ( 0 ) . Logistic regression is then applied using the candidate membership as a response on the set of covariates (again with the intercept suppressed) to obtain the starting value for α r ( 0 ) . This set of initials, θ ( 0 ) = ( μ r ( 0 ) , α r ( 0 ) , β ( 0 ) ) , are fed into the “NM” model and “SP” models for the EM algorithm to calculate θ ( k ) = ( μ r ( k ) , α r ( k ) , β ( k ) ) . The convergence of the algorithm is accessed by the criterion, ‖ θ ( k + 1 ) − θ ( k ) ‖ ‖ θ ( k ) ‖ ≤ ρ where ρ = 10 − 3 .

The NPMLE under shape constraint(s) is uniquely obtained as a formulation of the isotonic regression problem.

The results of the parameter estimations are presented in Table 1, in which “est” stands for estimate, and “[sd]” stands for estimated standard deviation. We observe that treatment-favorable subgroups are identified based on the profile of two baseline biomarkers, the parametric approach (NM) and the semiparametric approaches (SP1 and SP2). The estimates of β are comparable between the three models, but the NM approach produces different profiles, judging by the λ ‾ r , the mean values of estimated λ _r . Interestingly, between the semiparametric approaches, similar membership profiles are estimated, although the model with the symmetric assumption (SP1) produces slightly larger μ _r than the model where the asymmetric assumption is made (SP2). In general, the SP-based models produce larger μ _r values than the NM method. When a 20% threshold is used for λ ₀, both approaches warrant the hypothetical multi-stage trial proceeding. We tend to favor the results from the SP models because it does not make any distributional assumptions, except unimodality at zero. In the next section, we further study the methods with simulation.

5 Simulation study

Utilizing the data structure of the real trial analyzed above, we simulate n = 1000, with n ₁ = n ₂ = 500, 1-dimensional response y _i and covariates x r i = ( x r i , 1 , x r i , 2 ) . We first generate the covariates, sample the x _ri from the bivariate normal distribution with mean vector as ( 3.52 , 1.85 ) ′ and a given covariance matrix ( 0.75 − 0.01 − 0.01 0.09 ) . We generate the response data where we set the first n ₁ = 500 covariates as being assigned to treatment 1 and the rest to treatment 2. Given the covariates, we have that y r i = β ′ x r i + μ r + ϵ i with probability Δ r ( x r i | α r ) , and y r i = β ′ x r i + ϵ r i with probability 1 − Δ r ( x r i | α r ) . The values β , μ r , α r are taken from the real data analysis where Δ r ( x | α r ) = exp ( α r ′ x ) / ( 1 + exp ( α r ′ x ) ) . For the error term ɛ , we simulate two types of unimodal density: one is based on a Laplace distribution, which is symmetric, and one is based on a skewed normal distribution.

We also compare the performance of the estimates from the parametric model where a normal distribution is assumed for f ( ⋅ ) . We generate 200 simulation runs, and the results are summarized in Tables 2 and 3, where “cr” is the 95% coverage rate for the true parameter. Figure 1 provides a representative model fitting for the density functions during the simulation.

Figure 1:

Plot of estimated density function from the semiparametric approach (solid blue line) when the true distribution is Laplace (left panel) or skewed normal (right panel). The solid line represents the density of the generating distribution. The dashed line represents a normal curve from the parametric model.

Overall, the models produce the most accurate estimates for the regression coefficients β, followed by the estimates of μ r , and then the α r parameters. As expected, compared to the normal distribution-based model, the semiparametric method produces a much more accurate estimation with coverage rates very close to the nominal level.

6 Concluding remarks

We have proposed and studied a semiparametric model for targeted group sequential clinical trials in which the response-covariates relationship is specified by a linear regression model, and the error distribution is specified as nonparametric for model robustness, with the only assumption being unimodality. The semiparametric maximum likelihood estimate is used for model parameter estimation. Our simulation study shows that the proposed method has a clear advantage over the commonly used normal model, and then the method is used to analyze real clinical trial data. Here we focused on model parameter estimation. As studying the asymptotic distribution of the regression parameter estimation in the single index semiparametric model is still an open question, due to the fact that the shape restricted NPMLE of the nonparametric component is a non-smooth step function, it is challenging to obtain the asymptotic distribution of the corresponding Wald statistic. Thus currently the critical value or/and p-value of the test statistic can only be obtained via simulation. In addition it is known [13] that the likelihood ratio test does not exist for this problem. So the analytic form of the hypothesis test on existence of subgroups is a topic of our further study.