Keywords

Adaptive design, Biomarker targeted design, Maximum likelihood estimate, Non-parametric Bayes model, Sample size, Two-stage clinical trial

Introduction

Biomarker targeted design is an important step towards precision medicine, which can improve efficiency of randomized clinical trial [1,2]. Zhou, et al. [3] and Lee, et al. [4] proposed Bayesian approach and Wang, et al. [5] applied this design in therapeutic trials, Freidlin, et al. [6] discussed issues with this design, Tang and Zhou [7] proposed a general framework. The two or three stage designs are commonly used in recent phase IIb-IIIa trials comparing treatment vs. control. Such sequential monitoring has become an integral part of clinical trial. It allows for early stopping of the trial for extreme results observed in interim stage(s) [8-11].

More recently, Gao, Roy, and Tan [12,13] proposed a two-stage adaptive design for biomarker targeted population where only biomarker positive subjects enter the trial study. However, the test is imperfect subject to false positive and false negative errors. Thus a mixture normal model is used for the biomarker positive subjects. The final sample size re-estimation is based on the positive predict value (PPV), the proportion of true positives among test positives, estimated from the first stage data, deviates slight from a more common adaptive designs [14-16]. Proschan [17] and Xiong, Tan and Boyett [18] discussed sample size re-estimation in clinical trials. Since for early-stage clinical trial, often the sample size is relatively small, and so incorporating valuable prior information, if available, and using a robust method would be desirable. For the first goal a Bayesian model is preferable, while for robustness the nonparametric method is more suitable. Thus, we propose a nonparametric Bayesian method to achieve both goals.

The rest of this article is organized as follows. Section 2 introduces the problem, then develops the proposed nonparametric Bayesian method, and compares with the frequentist parametric method, the two-stage sequential test, and sample size re-estimation. In Section 3, extensive simulation studies are conducted to evaluate the performance of the proposed method, compared it with the existing normal parametric model. We leave the technical details in the Appendix.

Background and the Proposed Method

Background and review of the existing methods

We first introduce the setting by briefly reviewing the related targeted design in Gao, Roy, and Tan [13]. Consider a two-stage clinical trial with continuous or response endpoint comparing two groups, control and treatment, to assess the effect of a new treatment. The observed data at stage k is $D_{n_k}=\{x_{Ci},x_{Tj}:i=1,\dots ,n_{Ck},j=1,\dots ,n_{Tk}\}$ from independent biomarker test positive individuals (k = 1, 2), x_Ci is response from the i-th individual in the control group, and x_Tj is that from the control group. By convention, the sample size n_C2 and n_T2 include n_C1 and n_T1; n₁ = n_C1 + n_T1 is the planned sample size for stage I; and n₂ = n_C2 + n_T2 is the total sample size at end of the trial, which is subject to updating based on parameters estimates from stage I data. Let μ_C be the mean response of the control group, μ_T be that for the treatment group, and θ = μ_T - μ_C be the treatment difference in the overall study population. The objective of the trial is to test whether there is treatment difference between the two arms, i.e. to test H₀ : θ = 0 versus H_A : θ ≠ 0, with pre-specified significance level α and power β, and determine the total sample size n₂.

Since only biomarker test positive patients enter the trial, and a proportion ω of them are truly positive. Hence the trial population consists of true biomarker enrichment group (E) and the non-enrichment (NE) group. Let μ_0C and μ_1C be the mean response for NE and E group in control arm respectively; μ_0T and μ_1T be those for NE and E portion in the treatment arm respectively, and μ_C and μ_T be those of the control and treatment groups. Then

$${\mu }_C=\left(1-\omega \right){\mu }_{0C}+\omega {\mu }_{1C}\text{ and }{\mu }_T=\left(1-\omega \right){\mu }_{0T}+\omega {\mu }_{1T}.$$

In Gao, Roy and Tan [13], the following normal mixture model are used

$$X_{Ci}\sim (1-B_i)N({\mu }_C,{\sigma }^2)+B_{i}N({\mu }_C,{\sigma }^2)\sim N({\mu }_C,{\sigma }^2),(i=1,\dots ,n_C),$$

$$X_{Ti}\sim (1-B_i)N({\mu }_C,{\sigma }^2)+B_{i}N({\mu }_T,{\sigma }^2),(i=1,\dots ,n_T),$$

where B_i ~ Bernoulli (ω).

Let ${\overline{x}}_C^{(1)}$ and ${\overline{x}}_T^{(1)}$ be the sample means of the control and treatment group at end of stage I,

$${\widehat{\mu }}_C={\overline{x}}_C^{(1)}/\omega ,{\widehat{\mu }}_T^{(1)}={\overline{x}}_T^{(1)}/\omega ,$$

$${\widehat{\sigma }}_C^2=\frac{1}{n_{C1}-1}{\displaystyle \sum {(x_{Ci}{}^{(1)}-{\overline{x}}_C^{(1)})}^2},{\widehat{\sigma }}_T^2={\widehat{\sigma }}_C^2+\omega (1-\omega ){({\widehat{\mu }}_T^{(1)}-{\widehat{\mu }}_C^{(1)})}^2.$$

Then the total sample size n₂ needed for the whole trial is estimated as

$${\widehat{n}}_2=\frac{{(Z_{\alpha /2}+Z_{\beta })}^2{\sigma }_{T-C}^2}{{\mu }_{T-C}^2},$$

$$\begin{array}{l}where{\mu }_{T-C}=(1-\omega )({\widehat{\mu }}_{0T}-{\widehat{\mu }}_{0C})+\omega ({\widehat{\mu }}_{1T}-{\widehat{\mu }}_{1C})=\omega ({\widehat{\mu }}_T-{\widehat{\mu }}_C),and{\sigma }_{T-C}^2=\\ [{\sigma }^2+\omega (1-\omega ){(}^{{\widehat{\mu }}_{1T}}]+[{\sigma }^2+\omega (1-\omega ){(}^{{\widehat{\mu }}_{1C}}]=2{\sigma }^2+\omega (1-\omega ){(}^{{\widehat{\mu }}_T}.\end{array}$$

Maximum likelihood estimation (MLE) for mixture model is more conveniently obtained by the EM algorithm, given in the Appendix.

The proposed method

For early-stage clinical trial often the sample size is relatively small, and in some cases there is prior knowledge about the data distribution. A subjectively specified parametric model however may not be able to describe the distribution, so a nonparametric prior is preferred. Thus we adopt a nonparametric Bayesian model for this problem.

Let F_C and F_T be the distribution function of the control and treatment arm respectively. Often there are prior information for them. Let π(F_C) and π(F_T) be their priors, we assume $\pi (F_C)\sim D(P_C(\cdot )),$ and $\pi (F_T)\sim D(P_T(\cdot )),$ where $D(P_C(\cdot ))$ is the Dirichlet process with parameter PC(.); similarly for $D(P_T(\cdot )).$ The distribution PC is the prior knowledge about F_C, and similarly for PT.

We adopt the following assumptions assumed in Gao, Roy and Tan [13].

A1) μ_0C = μ_1C = μ_0T := μ_C;

A2) Var(X0C) = Var(X1C) = Var(X0T) = Var(X1T) := 2.

The reason for A1) is that a predictive biomarker is associated with response or lack of response to a particular therapy. Ideally, a predictive biomarker positive patient receiving therapy is expected to show a substantially higher response than negatively-biomarker patients receiving the therapy as well as those in the control group regardless of the marker status. Therefore, A1) with the treatment potentially making μ_1T different from μ_0C, μ_1C and μ_0T. Thus, the treatment effect, if any, is assumed to be a result of differential response to the treatment in the positively-biomarker group. A2) is a reasonable assumption to reduce model complexity.

For notational brevity, we will just write n_C for n_C,1 and n_T for n_T,1 etc. At end of stage I, using the non-parametric Bayesian formula for mean [19], in our case we have for μ_C and μ_T

$${\stackrel{⌣}{x}}_C^{(1)}=\frac{1}{n_C}{\mu }_{PC}+\frac{n{}_C}{n{}_C+1}{\overline{x}}_C^{(1)},{\stackrel{⌣}{x}}_T^{(1)}=\frac{1}{n_T}({\omega }_0{\mu }_{PT}+(1-{\omega }_0){\mu }_{PC})+\frac{n_T}{n_T+1}{\overline{x}}_T^{(1)},$$

where μ_PC and μ_PT are the prior means of P_C and P_T. We have the estimates for means of the two distributions

$${\stackrel{⌣}{\mu }}_C={\stackrel{⌣}{x}}_C^{(1)},{\stackrel{⌣}{\mu }}_T=\frac{{\stackrel{⌣}{x}}_T-(1-{\omega }_0){\stackrel{⌣}{x}}_C}{{\omega }_0}.$$

Also, using the non-parametric Bayesian formula for variance [19], we have ${\stackrel{⌣}{\sigma }}_C^2$ of Var(XC) and ${\stackrel{⌣}{\sigma }}_T^2$ of Var(XT ) in our case as

$${\stackrel{⌣}{\sigma }}_C^2=\frac{{\mu }_{2,P_C}+{\sum }_{i=1}^{n_C}{x}_{Ci}^2}{n_C+2}-\frac{{({\mu }_{P_C}+{\sum }_{i=1}^{n_C}{x}_{Ci})}^2}{(n_C+1)(n_C+2)},$$

$${\stackrel{⌣}{\sigma }}_T^2=\frac{{\mu }_{2,P_T}+n_T({\stackrel{⌣}{\sigma }}_C^2+{\omega }_0{\stackrel{⌣}{\mu }}_T+(1-{\omega }_0){\stackrel{⌣}{\mu }}_C)}{n_T+2}-\frac{{({\mu }_{P_T}+n_T({\omega }_0{\mu }_{P_T}+(1-{\omega }_0){\mu }_{P_C}))}^2}{(n_T+1)(n_T+2)},$$

where ${\mu }_{2,P_C}={\displaystyle \int x^{2}d{P}_C(x)}$ is the prior second moment, similarly for ${\mu }_{2,P_T}.$

To update the estimation of PPV ω, let $S^{(1)}=n_1{\overline{x}}_C^{(1)}+n_1{\overline{x}}_T^{(1)}$ be the total responses at end of stage I. Since $E[S^{(1)}]=\omega n_T{\mu }_T+[n_C+n_T(1-\omega )]{\mu }_C,$ we substitute means and total responses to estimate ω as

$$\stackrel{⌣}{\omega }=\frac{S^{(1)}-n_1{\stackrel{⌣}{\mu }}_C}{n_T({\stackrel{⌣}{\mu }}_T-{\stackrel{⌣}{\mu }}_C)}.$$

The two-stage sequential test

For given significance level α and power β, the decision boundaries are determined to satisfy type I error no greater than α and with power at least β. Consider test statistics T_j at stage j (j = 1, ..., k), the flexible class of boundaries proposed by Wang and Tsiatis [20] are, for some (c, γ) to be determined,

$$b_j=c{j}^{(\gamma -0.5)}.$$

Under H₀, (T₁, ..., T_k) is multivariate normal distributed with zero mean vector. The sequential test will reject the null hypothesis at stage j if |T_j| ≥ b_j. Then we have equation

$$P_{H_0}\left\{\underset{j=1}{\overset{k}{\cap }}\right|T_j\left|\text{ < }c{j}^{\gamma -0.5}\right\}=1-\alpha$$

Here we adopt O'Brien-Fleming boundaries having shape parameter γ = 0. In our case, k = 2, α = 0.05, c = 2.7967. The corresponding threshold of p-value at stage I is approximately 0.0054.

Sample size re-estimation

If the null hypothesis is not rejected at stage I, the sample size for the next stage need be determined by considering the required global power. Recall that n₁ and n₂ are the sample sizes at the end of stage I and II respectively. For simplicity we assume n_T = n_C = n₁/2. Then under H₀,

$$T_j=\sqrt{n_j}\frac{{\overline{x}}_T^{(j)}-{\overline{x}}_C^{(j)}}{2\sigma }\sim N(0,1),(j=1,2).$$

Under H₁,

$${(T_1,T_2)}^T\sim N(\mu ,\Omega ),$$

$$\mu =\left[\begin{array}{c}\frac{\sqrt{n_1}\omega ({\mu }_T-{\mu }_C)}{2\sigma }\\ \frac{\sqrt{n_2}\omega ({\mu }_T-{\mu }_C)}{2\sigma }\end{array}\right],\Omega =\left(1+\frac{\omega (1-\omega ){({\mu }_T-{\mu }_C)}^2}{2{\sigma }^2}\right)\left[\begin{array}{cc}1& \sqrt{\frac{n_1}{n_2}}\\ \sqrt{\frac{n_1}{n_2}}& 1\end{array}\right],$$

with family-wise power function

$$\beta =P_{H_1}\left(\left|T_1\right|\ge b_1\right)+P_{H_1}\left(\left|T_2\right|\ge \frac{b_2}{\sqrt{2}}\left|\left|T_1\right|<b_1\right.\right)$$

$$=P_{H_1}\left(\left|T_1\right|\ge b_1\right)+\frac{P_{H_1}\left(\left|T_1\right|<b_1\right)-P_{H_1}\left(\left|T_2\right|<\frac{b_2}{\sqrt{2}},\left|T_1\right|<b_1\right)}{P_{H_1}\left(\left|T_1\right|<b_1\right)}.$$

Denote $p_1=P_{H_1}\left(\left|T_1\right|\ge b_1\right)$ and $p(n_2)=P_{H_1}\left(\left|T_2\right|<b_2/\sqrt{2},\left|T_1\right|<b_1\right),$ which is determined by n₂ given other parameters. Then n₂ is determined by the least integer satisfies

$$p_1+1-\frac{p(n_2)}{1-p_1}\ge 1-\beta ,\text{ or }(p_1+\beta )(1-p_1)\ge p(n_2).$$

Simulation Study

Extensive simulation studies are conducted to compare three methods, the method of Gao, Roy and Tan [13], the proposed nonparametric Bayesian method, and the parametric maximum likelihood (MLE) based method. We considered both truncated normal and skewed normal (with skewness parameter α). The reason for the first is that the treatment effects should be negative values; the latter represent departure from the normality assumption. For the sequential test, we assume n₁ = n₂/2, i.e., the sample sizes are the same for the two stages. If the re-estimated sample size is less than n₂, we will keep the original design. Moreover, the estimation of ${\stackrel{⌣}{\sigma }}_T^2$ can be obtained using the formula below

$${\stackrel{⌣}{\sigma }}_T^2={\stackrel{⌣}{\sigma }}_C^2+{\omega }_0(1-{\omega }_0){({\stackrel{⌣}{\mu }}_T-{\stackrel{⌣}{\mu }}_C)}^2$$

which is the formula we actually used in the simulation study.

The results are displayed in Table 1, Table 2, Table 3, Table 4 and Table 5, with different parameter settings and sample sizes, with estimated treatment effects μ_T, μ_C, their difference θ, test result, and estimated sample size (if the null hypothesis is not rejected in stage I), the probability of reject H⁰ (so not continue the trial to stage II). In all the tables, results from the method of Gao, Roy and Tan [13] is named `empirical'; the proposed method, named `NP Bayes', and the MLE with EM-algorithm, named `MLE(EM), are compared. The sample size estimate n₂ depends on the estimated PPV is not stable. So we display the mean, median and trimmed mean (with trim proportion 10% on both sides) from all three methods, with 500 repetitions. Also, the corresponding disease prevalence (prev), sensitivity (sen), specificity (spec), the true ω₀, and the effect size (EF, |mean| /s.d) used in the simulation are given at the top of each Table.

To reflect the effects on estimated sample size of the skewness of the skewed normal and of EF, Figure 1, Figure 2, Figure 3 and Figure 4 are shown for the mean and trimed mean estimation methods below. We see that the empirical Bayesian estimate seems more reasonable, the empirical methods tend to over-estimate the sample size, and the MLE tends to under estimate it, and that the trimed mean method is much more stable. For comparison, we altered settings for parameters and kept EF the same for each graph.

From the Tables, the empirical and non-parametric Bayesian estimates are close to each other when the sample size is relatively large, which is consistent with the Bayes-frequentist estimation theory. However, the non-parametric Bayesian estimators were more robust than empirical ones, especially in cases where the sample size is relatively small. The variance of the empirical estimators was hugely inflated, which would require large increase of patients in stage-II. Moreover, as expected the non-parametric Bayesian estimators were not very sensitive to the selection of the prior, likely due to the fact that essential information was captured by the prior already. MLE method by EM algorithm can give smallest increase in sample size, especially when the sample were precisely from normal distribution, with zero-skewness. However, the estimators were significantly biased, which may be caused by the exceedingly high power of the stage-I test than the specified level. Therefore, non-parametric Bayesian estimators performed the best in the proposed adaptive design since they are more robust and require weak assumption on prior. They utilize the prior information but are not too dependent on prior selection. In smaller scales, it could reduce the unexpected variance, which may lead to large sample size increase in the next stage. In addition, they asymptotically converged to empirical estimators, indicating unbiasedness with large sample sizes.

Concluding Remarks

We have proposed a nonparametric Bayesian method for the two-stage adaptive biomarker-targeted clinical trial design. Compared with the existing parametric model, it has the advantage of incorporating prior information into the design, and being robust to model assumption. Extensive simulation studies are conducted to evaluate the performance of the proposed method, and compare with the commonly used methods for this problem. In our simulation studies, we considered the non-skewed and skewed, to reflect correct and incorrect model specifications. It was found that the skewness will influence the estimation accuracy in this design. The estimate is most accurate with left-skewed distributions, least accurate with right-skewed distributions, and modest with truncated normal distributions. Cases gave moderate results, and would be the hardest to estimate. Moreover, the more left-skewed, the more accurate for the estimations.

References

1. Simon R, Maitournam A (2004) Evaluating the efficiency of targeted designs for randomized clinical trials. Clin Cancer Res 10: 6759-6763.
2. Maitournam A, Simon R (2005) On the efficiency of targeted clinical trials. Statist Med 24: 329-339.
3. Zhou X, Liu S, Kim ES, Herbst RS, Lee JJ (2008) Bayesian adaptive design for targeted therapy development in lung cancer - a step toward personalized medicine. Clin Trials 5: 181-193.
4. Lee J, Gu X, Liu S (2010) Bayesian adaptive randomization designs for targeted agent development. Clin Trials 7: 584-596.
5. Wang SJ, Hung HMJ, O'Neill RT (2009) Adaptive patient enrichment designs in therapeutic trials. Biom J 51: 358-374.
6. Freidlin B, McShane LM, Korn EL (2010) Randomized clinical trials with biomarkers: Design issues. J Natl Cancer Inst 102: 152-160.
7. Tang L, Zhou XH (2013) A general framework of marker design with optimal allocation to assess clinical utility. Stat Med 32: 620-630.
8. Pocock SJ (1977) Group sequential methods in the design and analysis of clinical trials. Biometrika 64: 191-199.
9. O'Brien PC, Fleming TR (1979) A multiple testing procedure for clinical trials. Biometrics 35: 549-556.
10. Tan M, Xiong X, Kutner MH (1998) Clinical trial designs based on sequential conditional probability ratio tests and reverse stochastic curtailing, Biometrics 54: 682-695.
11. Jennison C, Turnbull B (2000) Group sequential methods with applications to clinical trials. Chapman and Hall, Boca Raton, Florida.
12. Gao Z, Roy A, Tan M (2015) Multistage adaptive biomarker-directed targeted design for randomized clinical trials. Contemp Clin Trials 42: 119-131.
13. Gao Z, Roy A, Tan M (2016) A two-stage adaptive targeted clinical trial design for biomarker performance-based sample size re-estimation. Statistics in Biosciences.
14. Proschan MA, Hunsberger SA (1995) Designed extension of studies based on conditional power. Biometrics 51: 1315-1324.
15. Lan KKG, Trost DC (1997) Estimation of parameters and sample size re-estimation. In: Proceedings of the Biopharmaceutical Section. American Statistical Association, 48-51.
16. Kieser M, Friede T (2000) Re-calculating the sample size in internal pilot study designs with control of the type I error rate. Stat Med 19: 901-911.
17. Proschan MA (2005) Two-stage sample size re-estimation based on a nuisance parameter: A review. J Biopharm Stat 15: 559-574.
18. Xiong X, Tan M, Boyett J (2003) Sequential conditional probability ratio tests for normalized test statistic on information time. Biometrics 59: 624-631.
19. Ferguson TS (1973) A Bayesian analysis of some nonparametric problems. Annals of Statistics 1: 209-230.
20. Wang SK, Tsiatis AA (1987) Approximately optimal one-parameter boundaries for group sequential trials. Biometrics 43: 193-199.

---