Longitudinal covariate-adjusted response-adaptive randomized designs

Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

1

Contents lists available at SciVerse ScienceDirect

3

Journal of Statistical Planning and Inference

5

journal homepage: www.elsevier.com/locate/jspi

7 9 11

Longitudinal covariate-adjusted response-adaptive randomized designs

13 15 Q1

Tao Huang a,b,c, Feifang Hu a,b,c,n

17

a b c

19 21 23 25 27

School of Statistics and Management, Shanghai University of Finance and Economics, Shanghai, PR China Department of Statistics, University of Virginia, Charlottesville, VA, USA School of Statistics, Renmin University of China, Beijing, PR China

a r t i c l e i n f o

abstract

Article history: Received 12 May 2011 Received in revised form 9 April 2013 Accepted 12 April 2013

Clinical trials often involve longitudinal data set which has two important characteristics: repeated and correlated measurements and time-varying covariates. In this paper, we propose a general framework of longitudinal covariate-adjusted response-adaptive (LCARA) randomization procedures. We study their properties under widely satisfied conditions. This design skews the allocation probabilities which depend on both patients' first observed covariates and sequentially estimated parameters based on the accrued longitudinal responses and covariates. The asymptotic properties of estimators for the unknown parameters and allocation proportions are established. The special case of binary treatment and continuous responses is studied in detail. Simulation studies and an analysis of the National Cooperative Gallstone Study (NCGS) data are carried out to illustrate the advantages of the proposed LCARA randomization procedure. & 2013 Published by Elsevier B.V.

Keywords: Asymptotic properties Clinical trial Continuous response Covariate-adjusted response-adaptive design Treatment allocation

29 31 33 35

63

37

65

39

1. Introduction

41

Clinical trials in health care are usually conducted with multiple but competing objectives. Among these objectives the efficacy, ethics and economics of the trial are often of primary concerns (Jennison and Turnbull, 2000). It is desirable that a randomization procedure can accommodate these competing objectives simultaneously. Response-adaptive randomization (RAR) procedures use sequentially accrued data on patient responses to skew the probability of patient allocation. Various approaches (Wei and Durham, 1978; Eisele, 1994; Ivanova, 2003; Hu and Rosenberger, 2003; Hu and Zhang, 2004) and their advantages have been well described in Hu and Rosenberger (2006). Often, information on covariates is available in addition to the primary outcome. Covariate-adjusted response-adaptive (CARA) (Rosenberger et al., 2001; Zhang et al., 2007) randomization procedures incorporate covariate information into the randomization procedure, and adjust the probability of patient allocation not only on the basis of previous responses but also on the current and past values of covariates. CARA randomization procedures are different from RAR procedures in that the allocation probability is determined without the assistance of some optimality criteria (Rosenberger and Sverdlov, 2008). Moreover, there are situations in practice where patients are monitored over a period of time and repeated responses are observed. In such situations, the accumulated information on longitudinal responses from active patients can provide additional information to assign incoming patients to the better treatment arm (Biswas and Dewanji, 2004a, 2004b; Sutradhar and Jowaheer, 2006).

67

43 45 47 49 51 53 55

Q3

61

71 73 75 77 79 81 83

57 59

69

n

Corresponding author at: Department of Statistics, University of Virginia, Charlottesville, VA, USA. E-mail addresses: [email protected] (T. Huang), [email protected] (F. Hu).

0378-3758/$ - see front matter & 2013 Published by Elsevier B.V. http://dx.doi.org/10.1016/j.jspi.2013.04.004

Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

85 87

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

2

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53

In this paper, we study the situations where information on both covariates and responses can be collected repeatedly over a period of time. There are several advantages to incorporate longitudinal covariates and response information into the randomization procedure. First, adaptive randomization designs based solely on the patient's response to treatment are inadequate and inappropriate when the probability distribution of the response is heterogeneous. For instance, a treatment may have quite different effects on patients with different characteristics, such as gender, age and smoking status, and it would be unreasonable to use data on the previous patient's response to determine the next patient's allocation probability without adjusting for these important covariates. Second, the probability distribution of the response can be time heterogeneous where the response to treatment presents a time trend (Bai and Hu, 2005; Duan and Hu, 2009). For instance, Grill et al. (1998) described a psychophysical experiment where the patients exhibit a clear learning curve in their response patterns. Under such scenarios, longitudinal data is intuitively appealing because it has the advantage to distinguish changes over time within individuals from differences among people in their baseline levels. Third, clinical trials are prospective studies in which staggered entrance and delayed response are common in practice. For longitudinal data, information is not only accumulating with each new patient, but also with each new observation from current patients. When treatment effects cannot be assessed immediately or the trajectory of treatment effects cannot be fully observed, the cumulative longitudinal information provided by multiple observations from current patients can still be utilized to better optimize the allocation of future patients. Finally, a randomization procedure under the longitudinal setting can simultaneously evaluate the treatment and covariate effects at the end of a trial, and then estimate and predict future responses. In this paper, we propose a general framework for longitudinal covariate-adjusted response-adaptive (LCARA) randomization procedures for a wide range of longitudinal data. The defining characteristic of a longitudinal data set is that each individual is measured repeatedly over time, and the sequential measurements on one individual are correlated. Consequently, the correlation structure in the data must be taken into account to draw valid inferences. Because of these characteristics of longitudinal data, the proposed LCARA randomization procedures are different from CARA randomization procedures (Zhang et al., 2007) in a number of ways. In particular, the proposed LCARA randomization procedures need to deal with (i) repeated and correlated measurements and (ii) time-varying covariates. It is especially interesting to investigate how these issues affect adaptive randomization procedures. Generalized estimation equation (GEE) (Liang and Zeger, 1986; Zeger and Liang, 1986) approach is often used for inference in generalized linear models because the likelihood sometimes is intractable or involves many nuisance parameters. A natural question is whether GEE approach is still applicable for inference in generalized linear models for longitudinal data when the longitudinal data is embedded in some adaptive randomization procedures. The answer is positive, and we will demonstrate that in general the GEE approach is still applicable for our proposed LCARA randomization procedure (see Theorem 1). In addition, in comparison with other covariates-adjusted randomization procedures, such as stratification and Pocock and Simon (1975) design, the proposed LCARA randomization procedures have the advantage to easily handle both discrete and continuous covariates. There are a limited number of work on adaptive randomization procedures for longitudinal data. For example, Biswas and Dewanji (2004a) and Sutradhar et al. (2005) extended the randomized play-the-winner rule (RPW) of Wei and Durham (1978) to longitudinal binary data with and without covariates, respectively, and Sutradhar and Jowaheer (2006) extended the RPW for longitudinal count data. Biswas et al. (2010) proposed a CARA design for binary longitudinal responses using the log-odds ratio within the Bayesian framework. These works are confined to the urn model and only for discrete longitudinal data set. Although some of these methods use the previous covariate information to determine the composition of the urn, they fail to incorporate the covariate information of an incoming patient to adjust the allocation probability. Therefore these procedures are not really covariate-adjusted response-adaptive randomization procedures. Moreover, as pointed in Hu and Rosenberger (2003), the randomized play-the-winner rule has high variability and is not recommended in practice. Our proposed LCARA randomization procedures can handle a wide range of longitudinal data: both discrete and continuous longitudinal responses as well as both discrete and continuous covariates. Asymptotic properties for the regression parameters and allocation proportions of LCARA randomization procedures are established to provide a statistical basis for inferences on treatment efficacy and covariate effects. Moreover, we apply the proposed LCARA randomization to binary treatment and continuous responses, and then evaluate the performance of the proposed LCARA randomization through extensive simulation studies and an analysis of the National Cooperative Gallstone Study data (Schoenfield et al., Q2 1981; Lachin 2007). The paper is organized as follows. Section 2 describes the general LCARA randomization procedures and their asymptotic properties. Examples of applications are presented in Section 3. Numerical studies and a real application are reported in Section 4. Conclusion remarks are in Section 5. The proofs are relegated to the Appendix.

63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99 101 103 105 107 109 111 113 115 117

55 2. Longitudinal covariate-adjusted response-adaptive randomization

119

57 2.1. General framework

121

59 61

The following notation and definitions are introduced to describe the general framework of LCARA randomization design. Given a clinical trial with K treatments, let X 1 ; X 2 ; … be the sequence of random treatment assignments. For the mth patient, Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

1 3 5 7 9 11 13 15

3

X m ¼ ðX m;1 ; …; X m;K Þ represents the treatment assignment such that if the patient is allocated to treatment k, then all elements in Xm are 0 except for the kth component, Xm,k, which is 1. Let Nm,k be the number of patients assigned to treatment k in the first m assignments, and write N m ¼ ðNm;1 ; …; N m;K Þ ¼ ∑m i ¼ 1 X i . In addition, assume that longitudinal covariate information is available in the study. For the mth patient, let

Z m ¼ ðZ Tm;1 ; …; Z Tm;Jm ÞT

be repeated covariate vectors recorded at a

sequence of time points T m ¼ ðt m;1 ; …; t m;Jm Þ. Correspondingly, let Y m;k ¼ ðY m;1;k ; …; Y m;Jm ;k ÞT be the response vector of the mth patient to treatment k. Denote Y m ¼ ðY m;1 ; …; Y m;K Þ. In practice, only the column Y m;k with X m;k ¼ 1 is observed. Let X m ¼ sðX 1 ; …; X m Þ, Z m ¼ sðZ 1 ; …; Z m Þ and Y m ¼ sðY 1 ; …; Y m Þ be the corresponding sigma fields. In addition, let F m ¼ sðX m ; Y m ; Z m Þ be the sigma field of the history. A general longitudinal covariate-adjusted response-adaptive randomization design is characterized by ψ mþ1 ¼ EðX mþ1 ∣F m ; Z mþ1;1 Þ ¼ EðX mþ1 ∣X m ; Y m ; Z m ; Z mþ1;1 Þ;

ð1Þ

the conditional probabilities of assigning treatments 1; …; K to the (m+1)th patient, conditioning on the entire history including the information of all previous m assignments, longitudinal responses and longitudinal covariates, plus only the information of the first observed covariate vector of the (m+1)th patient.

63 65 67 69 71 73 75 77

23

Remark 1. Allocation (1) only provides a general framework of LCARA randomization design, and we rigorously characterize the LCARA randomization procedure in Sections 2.2 and 2.3. Heuristically, in (1), the longitudinal data information are used in two successive stages. On the estimation stage, the entire history of information of all previous m patients, ðX m ; Y m ; Z m Þ, is used to estimate the unknown parameters that relate the longitudinal covariates to the longitudinal responses (see (2)). On the allocation stage, only the information of the first observed covariate vector of the (m+1)th patient, Zm+1, 1 (the only observed covariate vector at that time), is used to adjust the allocation probability ψ mþ1 (see (3)). In many applications, the covariate vector is fixed (or partial fixed) for each patient, namely Z m;1 ¼ ⋯ ¼ Z m;Jm .

25

2.2. Longitudinal model

87

27

Following the above notation, let Ym,t,k be the response to treatment k at time t of the mth patient with a covariate vector Zm,t at time t. We then assume for each treatment separately

89

17 19 21

29 31 33 35 37 39 41

ð2Þ

where pk ð; Þ are known functions, θk are unknown parameter vectors, and Θk is the parameter space of θk . Denote θ ¼ ðθ1 ; …; θK Þ and Θ ¼ Θ1  ⋯  ΘK . (ii) A parametric model for the conditional variance of Ym,t,k given Zm,t, namely Var½Y m;t;k ∣Z m;t  ¼ ϕbk ðpk ðθk ; Z m;t ÞÞ;

47 49

83 85

93 95 97 99

where bk ðÞ are known functions and ϕ is a possible unknown scale parameter. (iii) A parametric correlation structure on the repeated responses Y m;k ¼ ðY m;1;k ; …; Y m;Jm ;k ÞT , namely,

101

CorrðY m;k ÞÞ ¼ Σ k ðαk Þ; where Σ k ðÞ are known structures and αk are unknown parameter vectors.

103 105

43 45

81

91

(i) A parametric model for the conditional expectation of Y m;t;k given Z m;t , namely, E½Y m;t;k ∣Z m;t  ¼ pk ðθk ; Z m;t Þ;

79

Remark 2. The model (2) assumes that all treatment effects are different across the K treatment groups, and in the definition of LCARA allocation (Section 2.3), it is assumed that treatment effects for group k are estimated using data only from group k. This means that all covariates are predictive (contribute treatment–covariate interaction effects). In some applications, some covariates may be prognostic and contribute additive effects to the model (Simon, 2010). In this paper, we did not consider the case of common parameters. The assumption (ii) is satisfied for generalized linear model based on exponential family (see McCullagh and Nelder, 1989, pp. 29, 30). Assumption (iii) is a widely used condition in longitudinal model.

107 109 111 113

51 2.3. Allocation

115

53 55 57 59 61

Now we rigorously define the LCARA randomization procedure. To start, assign m0 patients to each treatment by using a restricted randomization (Rosenberger and Lachin, 2002; Efron, 1971). Assume that m ðm 4Km0 Þ patients have been assigned to treatments. Their longitudinal covariates fZ i ; i ¼ 1; …; mg and the corresponding responses fY i ; i ¼ 1; …; mg are observed. Let θbm ¼ ðθbm;1 ; …; b θ m;K Þ be an estimate of θ ¼ ðθ1 ; …; θK Þ. Here θbm;k ; k ¼ 1; …; K, is the estimator of θk based on the observed Nm,k samples (Zi, Yi,k) for which X i;k ¼ 1; i ¼ 1; …; m. When the (m+1)th patient is ready for randomization and the first covariate vector Zm+1, 1 is recorded, we assign this patient to treatment k with a probability of ψ mþ1;k ¼ PðX mþ1;k ¼ 1∣F m ; Z mþ1;1 Þ ¼ π k ðθbm ; Z mþ1;1 Þ;

ð3Þ

Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

117 119 121 123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

4

1 3 5 7 9 11 13 15 17 19 21 23

where π k ð; Þ, k ¼ 1; …; K, are some given functions. We call the function πð; Þ ¼ ðπ 1 ð; Þ; …; π K ð; ÞÞ the allocation function. It satisfies π 1 þ ⋯ þ π K ¼ 1. Remark 3. Note that, in (3), only the first observed covariate vector is used to adjust the allocation probability ψ mþ1;k , and all follow-up observations are unavailable in this allocation stage. Moreover, the proposed LCARA randomization procedures (3) are based on sequential estimation. In order to simplify our theoretical study, we assume that the entire longitudinal history of all previous m patients, namely F m , is used to estimate the unknown parameters θ. However, in practice, we may not fully observe the entire trajectory of treatment effects of each of these m patients, and for some patients, only a few observations are available (possible only one). If so, the partial longitudinal information provided by current patients can still be utilized to estimate the unknown parameters θ. In fact, the estimator of θ can and should be updated whenever a new observation becomes available. One can show that this will not affect the theoretical properties of the design as in Hu et al. (2008) under some additional conditions on how many observations are delayed for each patient. In this paper, our primary objective is to develop a general framework of LCARA randomization procedures, and we assume that the allocation function π is fixed. For a trial with some specific objectives, one can tailor the allocation function π based on some optimal criteria, for example, the hypothesis testing criterion as Tymofyeyev et al. (2007) or the D-optimal criterion as Atkinson (1982). For instance, assuming that Z and θk ; k ¼ 1; …; K, have the same dimensions (otherwise, slight modifications are necessary), we can take π k ðθ; ZÞ ¼ Rk ðθ1 Z T ; …; θk Z T Þ ¼ Rk ðu1 ; …; uK Þ ¼ Rk ðuÞ; k ¼ 1; …; K and u ¼ ðu1 ; …; uK Þ. Here, 0≤Rk ðuÞ≤1, k ¼ 1; …; K, are real functions that are defined in RK with K

∑ Rk ðuÞ ¼ 1;

k¼1

Ri ðuÞ ¼ Rj ðuÞ whenever ui ¼ uj :

Gðuk Þ ; Gðu1 Þ þ ⋯ þ GðuK Þ

65 67 69 71 73 75 77 79 81 83

For instance, functions Rk can be defined as Rk ðuÞ ¼

25

and

63

85

k ¼ 1; …; K; 87

33

where GðÞ is a smooth real function that is defined in R and satisfies 0≤GðÞ o∞. In the two-treatment case, we can let R1 ðu1 ; u2 Þ ¼ Gðu1 −u2 Þ and R2 ðu1 ; u2 Þ ¼ Gðu2 −u1 Þ, where GðÞ is a real function defined on R satisfying G(0)¼1/2, G(−u)¼G(u) and 0 o GðuÞ o 1 for all u. For instance, GðÞ can be any symmetric univariate cumulative distribution function. For the logistic regression model, Rosenberger et al. (2001) suggested using the estimated covariate-adjusted odds ratio, Rk ðu1 ; u2 Þ ¼ expðuk Þ=ðexpðu1 Þ þ expðu2 ÞÞ; k ¼ 1; 2, to allocate patients, which is equivalent to defining GðuÞ ¼ ð1 þ expð−uÞÞ−1 : When πðθ; ZÞ does not depend on Z, Bandyopadhyay and Biswas (2001) suggested choosing GðÞ ¼ Φð=sÞ for the normal linear regression model, where ΦðÞ is the cumulative normal distribution and s is a tuning parameter.

35

2.4. Asymptotic properties

97

37

To emphasize that the allocation function only depends on the first observed covariate vector, let Z ;1 denote the generic observed covariate vector at time one. Let νk ¼ g k ðθÞ ¼ E½π k ðθ; Z ;1 Þ; k ¼ 1; …; K, and denote gðθÞ ¼ ðg 1 ðθÞ; …; g K ðθÞÞ, ν ¼ ðν1 ; …; νK Þ. It follows from (3) that

99

27 29 31

39 41

PðX mþ1;k ¼ 1∣F m Þ ¼ g k ðθbm Þ;

k ¼ 1; …; K:

89 91 93 95

101 103

In addition, we assume 43

Assumption 1. Θ is a bounded domain in RKd and the true value θ is an interior point of Θ. Here d ¼ dimðΘk Þ for k ¼ 1; …; K.

105

45

Assumption 2. Z 1;1 ; …; Z n;1 are i.i.d. random vectors.

107

47

Assumption 3. For each fixed Z ;1 ¼ z;1 , π k ðθ; z;1 Þ 4 0 is a continuous function of θ; k ¼ 1; …; K.

109

49

Assumption 4. For k ¼ 1; …; K, π k ðθ; z;1 Þ is differentiable with respect to θ under the expectation, and there is a δ 40 such that T  ∂g k jθn g k ðθÞ ¼ g k ðθn Þ þ ðθ−θn Þ þ oð∥θ−θn ∥1þδ Þ; ∂θ

111

51 53 55 57 59 61

where ∂g k =∂θ ¼ ð∂g k =∂θ11 ; …; ∂g k =∂θKd Þ.

115 117

Theorem 1. Suppose that for k ¼ 1; …; K, 1 n ∑ X h ðY ; Z m Þð1 þ oð1ÞÞ þ oðn−1=2 Þ θbnk −θk ¼ n m ¼ 1 m;k k m;k

113

ð4Þ

119

where hk are K functions with E½hk ðY k ; ZÞ∣Z ¼ 0; and E∥hk ðY k ; ZÞ∥2 o∞. Then under Assumptions 1–4, we have, for k ¼ 1; …; K,

121

PðX n;k ¼ 1Þ-νk ;

a:s:;

PðX n;k ¼ 1jF n−1 ; Z n;1 ¼ z;1 Þ-π k ðθ; z;1 Þ

a:s:

ð5Þ

Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

1

5

63

and rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! Nn log log n −ν ¼ O ; n n

3

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! log log n b θ n −θ ¼ O : n

ð6Þ

65 67

5 Further, assume that, for k ¼ 1; …; K, 7

69

n

n−1 ∑ Efπ k ðθ; Z m;1 Þðhk ðY m;k ; Z m ÞÞT hk ðY m;k ; Z m Þg-V k : m¼1

9

71

13 15 17 19 21 23 25

35

N n;kjz;1 -π k ðθ; z;1 Þ a:s:; Nn ðz;1 Þ

53 55

79 81 83 85 87

k ¼ 1; …; K;

93 95 97

ð9Þ

99 101

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N D njz;1 −πðθ; z;1 Þ -Nð0; Σ jz;1 Þ; Nn ðz;1 Þ N n ðz;1 Þ

ð10Þ

103 105

where

Σ jz;1

0 1T ∂πðθ; z;1 Þ @ ∂πðθ; z;1 Þ A ¼ diagfπðθ; z;1 Þg−πðθ; z;1 Þ πðθ; z;1 Þ þ 2PðZ ;1 ¼ z;1 Þ ∑ Vk : ∂θk ∂θk k¼1 T

107

K

109 111

49 51

77

91

and

45 47

75

89

Theorem 2. Fix z;1 such that PðZ ;1 ¼ z;1 Þ 4 0. Under Assumptions 1–4 and Condition (4), we have

41 43

ð8Þ

where, for k≤K, N n;kjz;1 is the number of patients with first observed covariate z;1 who are randomized to treatment k by the first n random allocations, and Nn ðz;1 Þ is the total number of patients with covariate z;1 . Write N njz;1 ¼ ðN n;1jz;1 ; …; Nn;Kjz;1 Þ. The following theorem establishes some asymptotic properties of these proportions.

37 39

73

N n;kjz;1 ∑nm ¼ 1 X m;k IfZ m;1 ¼ z;1 g ¼ ; Nn ðz;1 Þ ∑nm ¼ 1 IfZ m;1 ¼ z;1 g

29

33

and Σ ¼ Σ 1 þ 2Σ 2 . Then,

In Condition (4), note that we assume that the entire longitudinal history of all previous m patients is used to estimate the unknown parameters. Condition (4) is assumed to facilitate the proof of Theorem 1, and it depends on different estimation methods. For generalized linear models, it is easy to check that the GEE estimator satisfies Condition (4) provided that the marginal regression models (2) are correctly specified. The correlation in the data is reflected in the expression of functions hk. Specifically, in many problems, hk ðY k ; ZÞ ¼ ðY k −pk ðθk ; ZÞÞT CovðY k Þ−1 ð∂pk ðθk ; ZÞ=∂θk ÞT (see Diggle et al., 2002, p. 139). Theorem 1 provides general results on the asymptotic properties of the allocation proportions Nn;k =n; k ¼ 1; …; K. Equally important, it provides general results on the asymptotic properties of the unknown parameters allowing to make inferences about the treatment and covariate effects at the end of clinical trials. Note that the asymptotic covariance matrices Σ and V in (8) depend on functions hk. In Appendix, via a special case, we show the expression of functions hk and also provide a simple and naive estimate for the covariance matrix V. Sometimes, one may be interested in the proportions for a particular stratum, which often corresponds to a discrete covariate. The proportion of patients from stratum z;1 who are assigned to treatment k is

27

31

Σ 2 ¼ ∑Kk ¼ 1 ð∂g=∂θk ÞV k ð∂g=∂θk ÞT ,

Let V ¼ diagfV 1 ; …; V K g, Σ 1 ¼ diagfν1 ; …; νk g−ν ν,   pffiffiffi Nn pffiffiffi D D n −ν -Nð0; ΣÞ and nðθbn −θÞ-Nð0; VÞ: n T

11

ð7Þ

Q4 Remark 4. Theorems 1 and 2 are different from Theorems 2.1 and 2.2 in Zhang et al. (2007) in the following aspects: (i) the

covariates Z 1 ; …; Z n are independent, but not identically distributed in our Theorems 1 and 2, while the covariates ξ1 ; …; ξn are i.i.d. in Zhang et al., 2007; (ii) the hk ðY m;k ; Z m Þ depends on both the correlated structure of the repeated measurement and the time-varying covariate Zm in our Theorems 1 and 2, while the hk ðY 1;k ; ξ1 Þ; …; hk ðY m;k ; ξm Þ are i.i.d. in Zhang et al. (2007); (iii) the allocation X m;k in (4) depends on Zm only through Zm,1 in our Theorems 1 and 2, while the Xm,k depends on ξm in Zhang et al. (2007). These differences will make the proof different and more difficult. The details are in the Appendix.

113 115 117 119

57 3. Comparing two treatments with continuous responses

121

59 61

In this section, the general results of Section 2 are applied to the special case of binary treatment and continuous responses. Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

6

1

3.1. Parameter estimation

63

3

Similar to Section 2.2, let Ym,t,k be the response variable of the mth patient with a covariate vector Zm,t to treatment k at time t. We then assume

65 67

5 (i) A linear model for the conditional expectation of Ym,t,k given Zm,t, namely, 7 9 11 13 15 17

ð11Þ

69

where ðμ1 ; μ2 Þ are the treatment effects, and ðβ1 ; β2 Þ are the covariate effects. Following the notation in Section 2, θ ¼ ðθ1 ; θ2 Þ ¼ ðμ1 ; β1 ; μ2 ; β2 Þ. Moreover for simplicity we assume that Z m ; m ¼ 1; 2; …; are not only independent, but also identically distributed. (ii) The random error process εm;t;k follows a normal distribution and

71

Y m;t;k ¼ μk þ βTk Z m;t þ εm;t;k ;

Eðεm;t;k Þ ¼ 0;

k ¼ 1; 2;

75

Covðεm;ti ;k ; εm;tj ;k Þ ¼ s2k Σðαk Þij ; k ¼ 1; 2;

where the correlation function Σðαk Þ has a parametric form, αk are unknown parameter vectors and Σðαk Þij is the corresponding (i,j) element of matrix Σðαk Þ. Common choices for the correlation function Σðαk Þ include compound symmetric (see the example in Section 4) and autoregressive (AR) or autoregressive moving average (ARMA) structures (see Diggle et al., 2002, Chapter 5 for details).

23 25 27 29 31 33 35 37 39 41

In practice, the maximum likelihood method can be used to estimate the parameters θ and η ¼ ðη1 ; η2 Þ ¼ ðs21 ; α1 ; s22 ; α2 Þ. In particular, one can simplify the computation using an iterative algorithm to estimate the parameters θ and η in a two-step fashion until the result converges. Step 1: For treatment k, k¼ 1, 2, given an estimate b η k for the parameter ηk , define a weight matrix Wðb η k Þ ¼ fb s 2k Σðb α k Þg−1 , we b then can obtain an estimate θ k for θk by using a weighted least squares method based on the observed samples (Zi,Yi,k) for which Xi,k ¼1. Step 2: For treatment k, k¼1, 2, given an estimate b θ k for the parameter θk , we can update the estimate for the parameter ηk based on the residuals derived from the observed samples (Zi, Yi,k) for which Xi,k ¼1.

The allocation scheme is as follows. Initially, enough patients are assigned to the two treatments by using a restricted randomization scheme so that the parameters in (11) can be accurately estimated. Assuming that m patients have been randomized, and their longitudinal responses and longitudinal covariates information are recorded, we then fit the b2 Þ. As discussed in Section 2.3, we regression model (11) and obtain the maximum likelihood estimator b θ ¼ ðb μ 1 ; βb1 ; μ b2 ; β specify the allocation function π for the two-treatment case as follows: π 1 ðθ; ZÞ ¼ 1−π 2 ðθ; ZÞ ¼ R1 ðθT1 Z; θT2 ZÞ ¼ GðθT1 Z−θT2 ZÞ ¼ Φðððθ1 −θ2 ÞT ZÞ=sÞ;

49 51 53 55 57 59 61

83 85 87 89

93 95 97 99

where ΦðÞ is the cumulative standard normal distribution and s is a tuning parameter. Henceforth an incoming (m+1)th patient with a covariate vector Zm+1,1 observed at time one is assigned to the treatment one with probability ! b 1 −β b2 ÞT Z mþ1;1 ðb μ 1 −b μ 2 Þ þ ðβ : ð12Þ pmþ1;1 ðZ mþ1;1 Þ ¼ Φ s

101 103 105

−1

Corollary 1. Suppose that Assumptions 1–4 are satisfied, Z 1 ; Z 2 ; … are i.i.d., and given a weight matrix WðηÞ ¼ fs ΣðαÞg , the matrix E½Z T WðηÞZ is nonsingular. We then have (5), (6), and (8) with V k ¼ fE½π k ðθ; Z ;1 ÞZ T Wðηk ÞZg−1 ; k ¼ 1; 2. Moreover, if PðZ ;1 ¼ z;1 Þ 40 for a given z;1 , then (9) and (10) hold. 2

47

79

91

3.2. LCARA allocation

43 45

77

81

19 21

73

Remark 5. From Corollary 1, it follows that qffiffiffiffiffiffiffiffiffi D Nn;k ðθbn;k −θk Þ-Nð0; νk V k Þ; k ¼ 1; 2:

107 109 111

ð13Þ

It should be noted that the asymptotic variance in (13) is different from that of linear models with a fixed allocation or noncovariate-adjusted randomization procedure. For the latter, we have qffiffiffiffiffiffiffiffiffi D ð14Þ Nn;k ðθbn;k −θk Þ-Nð0; fE½Z T Wðηk ÞZg−1 Þ; k ¼ 1; 2: Under a fixed allocation or non-covariate-adjusted randomization procedure, the asymptotic covariance matrix of estimators is fE½Z T Wðηk ÞZg−1 , because the allocation functions π k ðθ; Z ;1 Þ; k ¼ 1; 2, do not depend on the covariate Z ;1 , and π k ðθ; Z ;1 Þ ¼ g k ðθÞ ¼ νk . Remark 6. For Corollary 1, we assume that the weight matrix is the inverse of the true covariance matrix of random errors. It can be shown that, given an arbitrary weight matrix, say W, the asymptotic result (13) continues to hold, but with a Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

113 115 117 119 121 123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

1

5 7 9 11

νk fE½π k ðθ; Z ;1 ÞZ WZg−1 E½π k ðθ; Z ;1 ÞZ Ws2k Σðαk ÞWZfE½π k ðθ; Z ;1 ÞZ WZg−1 : T

T

65

It can also be shown that the choice of weight matrix W does not affect the consistency of estimation of θ, therefore the allocation function (12) is robust to the selection of the weight matrix. A simple and naive estimate of the asymptotic covariance matrix in (13) is, for k ¼ 1; 2, 8 9−1 > > < Nn;k = b 2k α k ÞZ m ; ∑ ωm Z Tm Σ −1 ðb s > > :m¼1 ; where

13 15

63

different asymptotic covariance matrix, T

3

7

ωm ¼

p^ m N

∑mn;k¼ 1 p^ m

;

p^ m ¼ Φ

67 69

ð15Þ

71 73

! b 1 −β b2 ÞT Z m;1 ðb μ 1 −b μ 2 Þ þ ðβ ; s

75 77

17

b1 ; μ b2 Þ and b2 ; β and b θ ¼ ðb μ1; β

19

4. Numerical study and applications

81

21

In this section, simulations are first used to evaluate the performance of our proposed LCARA randomization procedure. Then we demonstrate the proposed LCARA procedure through redesigning and analyzing a subset of data from the National Cooperative Gallstone Study (NCGS). For the simulation purpose, we consider the following linear model, given treatment k,

83

23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61

b b1 ; s b 22 ; α b2 Þ η ¼ ðb s 21 ; α

are MLE for θ and η at the end of randomization.

Y m;t;k ¼ θk;1 þ θk;2 T m;t þ θk;3 Z m;1 þ θk;4 Z m;2 þ εm;t;k ;

k ¼ 1; 2;

79

ð16Þ

where

87 89

θ1 ¼ ðθ1;1 ; θ1;2 ; θ1;3 ; θ1;4 Þ ¼ ð103; 1:3; 2:5; 1Þ; θ2 ¼ ðθ2;1 ; θ2;2 ; θ2;3 ; θ2;4 Þ ¼ ð100; 1; −1:5; 0:9Þ: Further, we assume that each patient has an initial observation at time zero once assigned to a treatment arm, and subsequently three repeated observations from time zero to time ten, specifically, Tm,1 ¼0, and T m;t ∼Uniform ½0; 10; t ¼ 2; 3; 4. Let Zm,1 be a binary indicator, say gender with one being male and zero being female, and P(Zm,1 ¼1) ¼ 0.5. Let Zm,2 be a continuous predictor uniformly distributed over [5, 7]. In addition, we assume that the random error process εmtk follows a normal distribution, and Eðεm;t;k Þ ¼ 0, Varðεm;t;k Þ ¼ s2k with ðs1 ; s2 Þ ¼ ð4:5; 3Þ. For both treatments, the correlation function has a compound symmetric structure, namely, Corrðεm;ti ;k ; εm;t j ;k Þ ¼ ΣðαÞij

85

where ΣðαÞ ¼ αJ þ ð1−αÞI; k ¼ 1; 2; and α ¼ 0:8:

Here the matrices I and J are the 4  4 identity and unit matrices respectively. We assume that there are n ¼200 patients to be allocated randomly. First, we randomly allocated 20 patients evenly into two arms by using a restricted randomization scheme. To show the improvement on the efficiency of the estimation of all parameters, we compare our proposed LCARA randomization scheme with CARA randomization scheme. Moreover, for the LCARA randomization scheme, we consider two weight matrices, one being the inverse of the true covariance matrix and the other being the identity matrix which ignores the correlation of repeated responses within a patient. We denote these two LCARA randomization schemes as LCARA-T and LCARA-I, respectively. We run 2000 iterations, with the tuning parameter s¼22, which results in ν ¼ ðν1 ; ν2 Þ ¼ ð0:60; 0:40Þ. Table 1 demonstrates the estimation of parameters at the end of randomization. In the table, “Bias” represents the sample average over 2000 iterations subtracting the true value, and “SD” represents the sample standard deviation over these 2000 estimates and can be thought as the true standard deviation of estimators. By comparing with the true values, the biases are very small for all parameters. It is observed that all three randomization schemes give accurate estimates for the parameters. But among them, LCARA-T has the smallest standard deviation and CARA has the largest standard deviation. This shows that our proposed LCARA randomization schemes (LCARA-T and LCARA-I) are more efficient than the CARA randomization scheme for estimating all parameters. We also use the proposed formula (15) to estimate the standard deviation of estimators of parameters. In the table, “Est Std” in the parentheses represents the average of 2000 estimates for standard deviation by using (15). For LCARA-T, it is observed that our proposed estimator (15) is quite accurate as the estimated standard deviation (Est Std) is very close to the true standard deviation (Std). For LCARA-I, the estimated standard deviation (Est Std) could be quite different from the true standard deviation (Std). This may due to the fact that the LCARA-I is based on an incorrect covariance matrix. Furthermore, by choosing the tuning parameter s¼44, 22, 14.4, 10.5, the target allocations are ν1 ¼ 0:55; 0:60; 0:65; 0:70, respectively. Table 2 shows that all three randomization procedures target the true allocations. However, among these three randomization schemes, CARA has the largest variability, which means that CARA randomization procedures approach Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

91 93 95 97 99 101 103 105 107 109 111 113 115 117 119 121 123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

8

1

63

Table 1

Q5 Parameter estimation.

3

LCARA-T

5

LCARA-I

Std (Est Std)

67

0.112

4.532 (4.318)

69

0.003

0.155 (0.153)

0.768 (0.403)

−0.001

0.821 (0.843)

0.381 (0.368)

−0.027

0.680 (0.707)

Bias

Std (Est Std)

Bias

Std (Est Std)

θb11 θb12

−0.117

1.234 (1.248)

−0.274

2.256 (2.130)

0.001

0.031 (0.030)

−0.002

0.053 (0.063)

9

θb13 θb14

0.011

0.783 (0.769)

0.010

0.014

0.187 (0.188)

0.036

11

θb21 θb22

13

θb23 θb24

7

65

CARA Bias

0.021

0.937 (0.958)

−0.139

1.585 (1.672)

−0.054

−0.001

0.024 (0.024)

−0.002

0.040 (0.049)

0.009

−0.060

0.563 (0.611)

−0.070

0.612 (0.368)

−0.028

0.720 (0.681)

0.014

0.157 (0.159)

0.026

0.270 (0.284)

−0.001

0.589 (0.561)

3.356 (3.528) 0.124 (0.116)

19

25

TP, s

TA, ν1

45 22 14.4 10.5

0.55 0.60 0.65 0.70

LCARA-T

LCARA-I

81

CARA

Estimate

Std

Estimate

Std

Estimate

Std

83

0.546 0.593 0.641 0.690

0.034 0.035 0.035 0.038

0.547 0.593 0.640 0.690

0.035 0.036 0.036 0.039

0.545 0.593 0.641 0.685

0.039 0.045 0.063 0.076

85

31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61

87 89

27 29

75

79

Table 2 Target allocation (TP, tuning parameter; TA, target allocation).

21 23

73

77

15 17

71

the target allocation more slowly than LCARA randomization procedures. This again shows that LCARA randomization procedures result in more efficient estimation of θ than CARA randomization procedures. We now redesign and analyze a subset of data from the National Cooperative Gallstone Study (NCGS) (Lachin, 2008) by using the proposed LCARA procedures. The ultimate objective of this study is to evaluate the safety and efficacy of the drug chenodiol for the treatment of cholesterol gallstones. In this study, each patient was randomly assigned to one of the three treatments, namely the high dose of 750 mg/day, the low dose of 375 mg/day, or the placebo, double-masked. In the NCGS it was suggested that chenodiol would dissolve gallstones by altering the metabolic pathway of cholesterol to desaturate cholesterol in gallbladder bile, but in doing so could theoretically increase level of serum cholesterol and predispose the patient to atherosclerosis. Each patient's serum cholesterol was thus monitored at baseline and roughly at 1, 2, 3, 6, 9, 12, 16, 20 and 24 months of follow-up, but substantial measurements were missing for various reasons. For illustration, we only consider those patients who were assigned to the placebo and the high dose treatments with at least three repeated longitudinal observations on serum cholesterol. In total, there are total 591 patients, 295 from the high dose group and 296 from the placebo group. On average, 8.943 observations per patient. The side effect, the increment of serum cholesterol, Ymtk, of using the drug chenodiol is our primary concern. We are interested in adjusting various covariates in assigning patients. In particular, five covariates, weight, age, gender, smoking status and drinking status, are considered. Of these 591 patients, 324 (54.82%) are female, 132 (22.34%) are smokers, 240 (40.61%) are drinkers. In addition, we assume that the treatment effects are linear, namely baseline treatment effects θ1 and treatment effects over time θ2 . The model can then be written as E½Y m;t;k ∣Z m;t  ¼ θk;1 þ θk;2 T þ θk;3 Weight þ θk;4 Age þ θk;5 Gender þ θk;6 Smoking þ θk;7 Drinking: The NCGS randomization sequence was implemented using a modified biased-coin design, however, retrospectively, we redesign the study by assuming that the NCGS randomization sequence was implemented using different procedures. In particular, we consider two randomization procedures, our proposed LCARA randomization procedure and the CARA randomization procedure. For LCARA randomization procedures, we consider two covariance matrices, a compound symmetric matrix with α ¼ 0:51 based on the real data and an identity matrix. For all procedures, we first assign 20 patients to each treatment using a restricted randomization scheme, then adaptively allocate another 551 patients. We target a sequence of allocation proportions, ν1 ¼ ð0:40; 0:42; …; 0:48Þ, for the high dose treatment, and perform 2000 replications. Similar to the simulation studies, we report the parameter estimation for ν1 ¼ 0:40 in Table 3 and the target allocation in Table 4. Both tables show that the proposed LCARA randomization procedures result in more efficient estimations for all parameters. Now we compare our proposed LCARA randomization procedure with Neyman randomization procedure based on the simulation model (16). In particular, we consider the null hypothesis, H 0 : Δ ¼ θ1;1 −θ2;1 −3 ¼ 0, and evaluate the power in a sequence of alternatives with Δ ¼ 0:5; 1:0; …; 4:5. In Table 5, for both allocation procedures, we report the estimated allocation proportions for treatment one, the power and the proportions male assigned to treatment one (MIO: Male in One) Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

91 93 95 97 99 101 103 105 107 109 111 113 115 117 119 121 123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

1

LCARA-T

5

9 11 13 15 17 19

63

Table 3 Parameter estimation.

3

7

9

LCARA-I

Bias

Std (Est Std)

65

CARA

Bias

Std (Est Std)

Bias

Std (Est Std)

67 69

θb11 θb12

1.139

15.628 (14.769)

1.302

16.098 (6.936)

0.372

20.986 (19.632)

0.000

0.068 (0.068)

0.000

0.079 (0.096)

0.047

0.288 (0.271)

θb13 θb14 θb15

−0.007

0.120 (0.121)

−0.008

0.121 (0.056)

−0.016

−0.105

3.831 (3.754)

−0.098

3.756 (1.803)

0.257

5.437 (4.983)

−0.010

0.173 (0.153)

−0.011

0.171 (0.080)

−0.009

0.218 (0.208)

θb16 θb17

−0.463

4.093 (4.124)

−0.604

4.309 (2.089)

−1.103

6.843 (6.530)

−0.402

3.420 (3.367)

−0.420

3.613 (1.590)

−0.110

4.462 (4.433)

θb21 θb22

−1.113

13.563 (13.278)

−1.201

13.583 (6.452)

−1.389

18.067 (17.661)

0.004

0.078 (0.083)

0.003

0.236 (0.247)

θb23 θb24

0.004

0.101 (0.097)

0.006

0.110 (0.049)

0.000

0.128 (0.121)

0.103

3.622 (3.350)

0.125

3.427 (1.512)

0.250

4.487 (4.142)

0.004

0.067 (0.063)

0.173 (0.161)

θb25 θb26

0.006

0.149 (0.147)

0.005

0.157 (0.073)

0.014

0.201 (0.199)

−0.031

3.319 (3.401)

−0.083

3.219 (1.501)

−0.005

5.023 (4.380)

θb27

−0.012

3.102 (2.909)

−0.102

3.234 (1.420)

0.170

4.106 (3.983)

25

Table 4 Target allocation (TP, tuning parameter; TA, target allocation). TP, s

TA, ν1

27 29 31

37

26 33 44 66 133

0.40 0.42 0.44 0.46 0.48

LCARA-T

LCARA-I

CARA

43 45 47 49 51 53 55

77 79 81

87

Estimate

Std

Estimate

Std

Estimate

Std

0.412 0.426 0.442 0.462 0.481

0.048 0.041 0.032 0.025 0.020

0.411 0.427 0.442 0.461 0.480

0.048 0.041 0.033 0.024 0.020

0.407 0.422 0.448 0.464 0.482

0.086 0.072 0.058 0.036 0.027

89 91 93 95

Table 5 LCARA and Neyman allocations. Δ

97

Allocation

Power

Covariate-Adjust, MIO

LCARA

Neyman

LCARA

Neyman

LCARA

Neyman

0.597 0.594 0.596 0.595 0.594 0.597 0.594 0.596 0.597 0.592

0.601 0.602 0.602 0.601 0.602 0.605 0.602 0.604 0.603 0.601

0.040 0.059 0.118 0.147 0.231 0.368 0.480 0.606 0.727 0.820

0.044 0.053 0.117 0.143 0.232 0.345 0.481 0.605 0.711 0.816

0.529 0.529 0.527 0.523 0.519 0.522 0.513 0.515 0.517 0.517

0.503 0.501 0.500 0.501 0.502 0.504 0.496 0.503 0.494 0.497

99 101

39 41

75

85

33 35

73

83

21 23

71

0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5

based on 2000 replications. It is observed that both our proposed LCARA procedure and the Neyman procedure approach the true target allocation 0.60, and that they have similar power. But between the two allocation procedures, our proposed LCARA method allocates more male patients (larger than 50 percentage) to treatment one, which is consistent with θ13 ¼ 2:5 4θ23 ¼ −1:5. In other words, in comparison to fixed allocation procedures or non-covariate-adjusted random procedures when targeting the same allocation proportions, our proposed LCARA procedures can assign more patients to the better or more appropriate treatment according to the covariate information, still preserving power.

103 105 107 109 111 113 115 117

57

5. Discussion

119

59

In this paper, we have proposed a general framework of longitudinal covariate-adjusted response-adaptive randomization procedures, which apply to discrete and continuous responses as well as discrete and continuous covariates. The proposed LCARA designs incorporate the special features of longitudinal data, such as repeated and correlated

121

61

Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

10

1 3 5 7 9 11 13 15 17 19 21 23

measurements, and time-varying covariates. Asymptotic properties are established to provide a statistical basis for inferences on treatment efficacy and covariate effects under widely satisfied conditions. For parameter estimation from a trial using LCARA randomization, asymptotic normality continues to hold under these conditions. However, the asymptotic variance is different from that of fixed allocation or non-covariate-adjusted randomization. From simulation studies and an analysis of the NCGS data, we can see the advantages of using LCARA randomization procedures. In comparison to fixed allocation or non-covariate-adjusted randomization, LCARA randomization balances the allocation across all covariates, and assigns more patients to the better treatment arm for a given covariate. In comparison to CARA randomization, LCARA randomization significantly improves the estimation efficiency of all regression parameters. Thus, ethically, LCARA randomization procedures may reduce the number of participants needed for a given study. The theoretical and numerical properties provide a guideline to apply LCARA randomization procedures in practice. In the simulated study (with sample size n ¼200), we start the LCARA procedures after first 20 patients (m0 ¼10). In the study of the NCGS data (sample size n¼591), we start the LCARA procedures after first 40 patients (m0 ¼20). In both studies, we estimate the parameters by using GEE method. In practice, one may choose a suitable m0 based on the number of covariates, the correlation structure and the total sample size. Several issues not addressed here could be worthy topics of future research. One question is how to choose a tuning parameter s in (12) to target a specific allocation proportion. The proposed LCARA procedures are based on a given allocation function π k ðθ; Z ;1 Þ. It is unclear how to obtain optimal π k ðθ; Z ;1 Þ under longitudinal setting. It is also interesting and important to investigate the behavior of the power function when a LCARA randomization procedure is used in practice. In general, the theoretical formulation can be really complicated because it involves the target allocation νk , the allocation function π k ðθ; Z ;1 Þ, the covariate Z as well as the variance Vk. However, this paper provides a framework and a starting point to study these topics. Missing data is quite common in clinical studies. We have not considered subject dropout and missing in this paper. Both optimal target allocations (Biswas and Mandal, 2004; Zhang and Rosenberger, 2006) and optimal designs (Hu et al., 2006, 2009) have been well studied in the literature of response-adaptive randomization. We may apply these results to LCARA procedures as well.

63 65 67 69 71 73 75 77 79 81 83 85 87

25 27

Acknowledgments

89

29

91

31

The authors thank Professor Lachin for providing the NCGS data. Special thanks go to anonymous referees, the associate editor and the editor for the constructive comments, which led to a much improved version of the paper. This research work is partially supported by the National Science Foundation Grants DMS-0906661, DMS-0907297 and DMS–1209164.

93

33

Appendix A. Technical proofs

95

35

Proofs of Theorems 1 and 2. To prove Theorem 1, we first prove the asymptotic properties of θ^ n . Let F m ¼ sðX m ; Y m ; Z m Þ be the s field of the history as defined in Section 2 and the symbol ∇ denote the differencing operator of a sequence fxn g, i.e., ∇xn ¼ xn −xn−1 . Further let ∇T m;k ¼ X m;k hk ðY m;k ; Z m Þ, k ¼ 1; …; K. Because E½hk ðY k ; ZÞ∣Z ¼ 0, we have

97

37 39 41 43 45 47 49 51 53 55

E½ðX m;k hk ðY m;k ; Z m Þ∣F m−1  ¼ E½fE½ðX m;k hk ðY m;k ; Z m Þ∣F m−1 ∣Z m  ¼ E½πðθ^ m−1 ; Z m;1 ÞE½hk ðY m;k ; Z m Þ∣Z m  ¼ 0:

101

Therefore Tn,k is a martingale sequence. Further because ∥X m;k hk ðY m;k ; Z m Þ∥≤∥hk ðY m;k ; Z m Þ∥, E∥hk ðY k ; ZÞ∥2 o ∞ and the law of the iterated logarithm for martingales, we have pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T n;k ¼ Oð n log log nÞ a:s: for k ¼ 1; …; K: From Eq. (4), we have rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! log log n ^θ n −θ ¼ O n

a:s:

99

ð17Þ

From the above derivation, it is easy to see that the Lindberg condition is satisfied for the martingale sequence T n;k ðk ¼ 1; …; KÞ. We just need to obtain the conditional variance–covariance matrices of the martingale difference f∇T m;k g satisfying E½ð∇T m;k ÞT ∇T m;k ∣F m−1  ¼ E½πðθ^ m−1 ; Z m;1 Þðhk ðY m;k ; Z m ÞÞT hk ðY m;k ; Z m Þ: Because of condition 2 of Assumption 1, and Eq. (17), we have

103 105 107 109 111 113 115 117

n

57 59

n−1 ∑ E½πðθ^ m−1 ; Z m;1 Þðhk ðY m;k ; Z m ÞÞT hk ðY m;k ; Z m Þ

119

m¼1

n

−n−1 ∑ E½πðθ; Z m;1 Þðhk ðY m;k ; Z m ÞÞT hk ðY m;k ; Z m Þ⟶0

121

m¼1

123

61 Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

1

11

63

in probability. Further, by condition (7), we have n

n−1 ∑ E½πðθ^ m−1 ; Z m;1 Þðhk ðY m;k ; Z m ÞÞT hk ðY m;k ; Z m Þ⟶V k

3

65

m¼1

5

in probability. Therefore we obtain the asymptotic normality of θ^ n . We may write

7 9

15 17 19

for k ¼ 1; …; K, where θ^ m is the estimator of θ based on F m . In general, we have

25

n−1

m¼1

m¼1

71 ð18Þ

73

Based on the results we obtained about θ^ m and the fact that Z 1;1 ; …; Z n;1 are i.i.d. random vectors, the rest of the proof follows from the proof of Theorem 2.1 in Zhang et al. (2007). The proof of Theorem 2 follows from the proof of Theorem 1 and the proof of Theorem 2.2 in Zhang et al. (2007). We omit the details here. □

75

Proof of Corollary 1. For the observations up to stage m, let θbm;k minimize the error sum of squares m

T

Sk ðθk Þ ¼ ∑ X i;k ðY i;k −θk Z i Þ Wðηk ÞðY i;k −θk Z i Þ i¼1

over θk ∈Θk ;

m

m

i¼1

i¼1

83

ð20Þ

m

m

i¼1

i¼1 m

¼ m−1 ∑ E½π k ðθ; ZÞZ T WðηÞZ∣

θ ¼b θ i−1 ;η ¼ b η i−1

i¼1

35

87 89 91 93

m−1 ∑ X i;k Z Ti Wðηk ÞZ i ¼ m−1 ∑ E½X i;k Z Ti Wðηk ÞZ i ∣F i−1  þ oð1Þ

33

85

T

Note that and fX i;k ðY i;k −θk Z i Þ Wðηk ÞZ i g are both sequences of martingale differences. It follows from the law of large numbers for martingales that

31

79

ð19Þ

ðθbm;k −θk Þm−1 ∑ X i;k Z Ti Wðηk ÞZ i ¼ m−1 ∑ X i;k ðY i;k −θk Z i ÞT Wðηk ÞZ i fX i;k Z Ti Wðηk ÞZ i −E½X i;k Z Ti Wðηk ÞZ i ∣F i−1 g

77

81

for k ¼ 1; …; K, where Wðηk Þ is a weight matrix. Note that Sk ðθk Þ is strictly convex function of θk . It follows that the weighted least squares estimator θbk for θk exists and is unique. Moreover, it is easy to show that θbm;k is the solution of the normal equation as

27 29

n

Nn;k −nνk ¼ ðE½X 1;k jF 0 −g k ðθÞÞ þ ∑ ðX m;k −E½X m;k jF m−1 Þ þ ∑ ½g k ðθ^ m Þ−g k ðθÞ:

21 23

69

X mþ1;k ¼ X mþ1;k −E½X mþ1;k jF m  þ g k ðθ^ m Þ;

11 13

67

þ oð1Þ

a:s:

ð21Þ

95 97

and 37

99

m

m−1 ∑ X i;k ðY i;k −θk Z i ÞT Wðηk ÞZ i -0

a:s:

i¼1

39

101

It follows that 41 ∥b θ m;k −θk ∥2

 min

∥b∥2 ¼ 1;θ∈Θk

43

"

 T bðE½π k ðθ; ZÞZ T WðηÞZÞb þ oð1Þ

m

103 105

#

45

η k ÞZ i ðθbm;k −θk ÞT ≤ðb θ m;k −θk Þm−1 ∑ X i;k Z Ti Wðb

47

¼ m−1 ∑ X i;k ðY i;k −θk Z i ÞT Wðη^ k ÞZ i ðb θ m;k −θk ÞT

109

¼ oð1Þ∥b θ m;k −θk ∥

111

107

i¼1

m

i¼1

49 51 53

a:s:

Hence

113 b θ m;k -θk

ð22Þ

a:s:

115

Now, by (21) and (22), 55 57 59 61

1 m ∑ X Z T Wðηk ÞZ i ¼ I Z;k þ oð1Þ m i ¼ 1 i;k i

117 a:s:; 119

where I Z;k ¼ E½π k ðθ; Z ;1 ÞZ T Wðηk ÞZ. Which, together with (20), implies that 1 m b θ m;k −θk ¼ ∑ X i;k ðY i;k −θk Z i ÞT Wðηk ÞZ i I −1 Z;k ð1 þ oð1ÞÞ m i

121

a:s:

Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

123

T. Huang, F. Hu / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]

12

1 3

Note that E½ðY i;k −θk Z i ÞT Wðηk ÞZ i I −1 Z;k ∣Z i  ¼ 0. Eq. (4) is satisfied with hk ¼ ðY k −θk ZÞ

T

Wðηk ÞZI −1 Z;k :

53 ð23Þ

55

Now based on Theorem 1 and (23), 5 7 ¼n 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51

m

57

i¼1 m

59

i¼1

61 63

V k ¼ n−1 ∑ Efπ k ðθ; Z i;1 Þðhk ðY i;k ; Z i ÞÞT hk ðY i;k ; Z i Þg þ oð1Þ −1

T −1 T ∑ Efπ k ðθ; Z i;1 ÞI −1 Z;k Z i Wðηk ÞðY i;k −θ k Z i ÞðY i;k −θ k Z i Þ Wðηk ÞZ i I Z;k g þ oð1Þ

T −1 −1 ¼ I −1 Z;k E½π k ðθ; Z ;1 ÞZ Wðηk ÞZI Z;k þ oð1Þ ¼ I Z;k þ oð1Þ;

where Wðηk Þ ¼ CovðY k Þ−1 . We proved Corollary 1.

□

References Atkinson, A.C., 1982. Optimal biased coin designs for sequential clinical trials with prognostic factors. Biometrika 69, 61–67. Bai, Z.D., Hu, F., 2005. Asymptotics in randomized urn models. Annals of Applied Probability 15, 914–940. Bandyopadhyay, U., Biswas, A., 2001. Adaptive designs for normal responses with prognostic factors. Biometrika 88, 409–419. Biswas, A., Dewanji, A., 2004a. A randomized longitudinal play-the-winner design for repeated binary data. Australian and New Zealand Journal of Statistics 46, 675–684. Biswas, A., Dewanji, A., 2004b. Sequential adaptive designs for clinical trials with longitudinal responses. Statistical Theory and Method Abstracts 173, 69–84. Biswas, A., Mandal, S., 2004. Optimal adaptive designs in phase III clinical trials for continuous responses with covariates. In: Bucchianico, A.D.i, Läuter, H., Wynn, H.P. (Eds.), m0Da 7-Advances in Model-Oriented Design and Analysis. Physica-Verlag, Heidelberg, pp. 51–59. Biswas, A., Park, E., Bandyopadhyay, R., 2010. Covariate-adjusted response-adaptive designs for longitudinal treatment responses: PEMF trial revisited. Statistical Methods in medical Research 21 (4), 379–392. Diggle, P.J., Heagerty, P., Liang, K.Y., Zeger, S.L., 2002. Analysis of Longitudinal Data. Oxford University Press Inc. Duan, L., Hu, F., 2009. Doubly adaptive biased coin designs with heterogeneous responses. Journal of Statistical Planning and Inference 139, 3220–3230. Efron, B., 1971. Forcing a sequential experiment to be balanced. Biometrika 58, 403–417. Eisele, J.R., 1994. The doubly adaptive biased coin design for sequential clinical trials. Journal of Statistical Planning and Inference 38, 249–261. Grill, S.E., Rosenberger, W.F., Boyle, K., Cannon, M., Hallet, M., 1998. Perception of timing of kinesthetic stimuli. Neuro Report 9, 4001–4005. Hu, F., Rosenberger, W.F., 2003. Optimality, variability, power: evaluating response-adaptive randomization procedures for treatment comparisons. Journal of the American Statistical Association 98, 671–678. Hu, F., Rosenberger, W.F., 2006. The Theory of Response-Adaptive Randomization in Clinical Trials. Wiley. Hu, F., Rosenberger, W.F., Zhang, L.X., 2006. Asymptotically best response-adaptive randomization procedures. Journal of Statistical Planning and Inference 136, 1911–1922. Hu, F., Zhang, L.X., 2004. Asymptotic properties of doubly adaptive biased coin designs for multi-treatment clinical trials. Annals of Statistics 30, 268–301. Hu, F., Zhang, L.X., Cheung, S.H., Chan, W.S., 2008. Doubly adaptive biased coin designs with delayed responses. Canadian Journal of Statistics 36, 541–559. Hu, F., Zhang, L.X., He, X., 2009. Efficient randomized adaptive designs. Annals of Statistics 37, 2543–2560. Ivanova, A., 2003. A play-the-winner type urn design with reduced variability. Metrika 58, 1–14. Jennison, C., Turnbull, B.W., 2000. Group Sequential Methods with Applications to Clinical Trials. Chapman and Hall, CRC. Lachin, J.M., 2008. National Cooperative Gallstone Study. In: D'Agostino, R., Sullivan, L., Massaro, J. (Eds.), Wiley Encyclopedia of Clinical Trials, pp. 1–6. Liang, K.Y., Zeger, S.L., 1986. Longitudinal data analysis using generalized linear models. Biometrika 73, 12–22. McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models. Chapman and Hall, CRC. Pocock, S.J., Simon, R., 1975. Sequential treatment assignment with balancing for prognostic factors in the controlled clinical trial. Biometrics 31, 103–115. Rosenberger, W.F., Lachin, J.M., 2002. Randomization in Clinical Trials: Theory and Practice. Wiley, New York. Rosenberger, W.F., Sverdlov, O., 2008. Handling covariates in the design of clinical trials. Statistical Science 23, 404–419. Rosenberger, W.F., Vidyashankar, A.N., Agarwal, D.K., 2001. Covariate-adjusted response-adaptive designs for binary response. Journal of Biopharmaceutical Statistics 11, 227–236. Schoenfield, L.J., Lachin, J.M., The Steering Committee, & the NCGS Group, 1981. Chenodiol (Chenodeoxycholic Acid) for dissolution of gallstones: the National Cooperative Gallstone Study. Annals of Internal Medicine 95, 257–282. Simon, R., 2010. Clinical trials for predictive medicine: new challenges and paradigms. Clinical Trials 7, 516–524. Sutradhar, B.C., Biswas, A., Bari, W., 2005. Marginal regression for binary longitudinal data in adaptive clinical trials. Scandinavian Journal of Statistics 32, 93–113. Sutradhar, B.C., Jowaheer, V., 2006. Analyzing longitudinal count data from adaptive clinical trials: a weighted generalized quasi-likelihood approach. Journal of Statistical Computation and Simulation 76, 1079–1093. Tymofyeyev, Y., Rosenberger, W.F., Hu, F., 2007. Implementing optimal allocation in sequential binary response experiments. Journal of the American Statistical Association 102, 224–234. Wei, L.J., Durham, S., 1978. The randomized play-the-winner rule in medical trials. Journal of the American Statistical Association 73, 840–843. Zeger, S.L., Liang, K.Y., 1986. Longitudinal data analysis for discrete and continuous outcomes. Biometrics 42, 121–130. Zhang, L.X., Hu, F., Cheung, S.H., Chan, W.S., 2007. Asymptotic properties of covariate-adjusted adaptive designs. Annals of Statistics 35, 1166–1182. Zhang, L., Rosenberger, W.F., 2006. Response-adaptive randomization for clinical trials with continuous outcomes. Biometrics 62, 562–569.

Please cite this article as: Huang, T., Hu, F., Longitudinal covariate-adjusted response-adaptive randomized designs. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.004i

65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99 101 103