# Additive Hazard Regression for the Analysis of Clustered Survival Data from Case-Cohort Studies

Authors: June Liu; Yi Zhang
Journal: Journal of Mathematics
 (2020)

---

## Abstract

The case-cohort design is an effective and economical method in large cohort studies, especially when the disease rate is low. Case-cohort design in most of the existing literature is mainly used to analyze the univariate failure time data. But in practice, multivariate failure time data are commonly encountered in biomedical research. In this paper, we will propose methods based on estimating equation method for case-cohort designs for clustered survival data. By introducing the event failure rate, three different weight functions are constructed. Then, three estimating equations and parameter estimators are presented. Furthermore, consistency and asymptotic normality of the proposed estimators are established. Finally, the simulation results show that the proposed estimation procedure has reasonable finite sample behaviors.


---

## Body

## 1. Introduction

In many failure time studies, there is a natural clustering of study subjects such that failure times within the same cluster may be correlated. One example of clinical trial is that the failure times of patients within a famlily may be correlated because they have common genetic characteristics and environmental factors. Another example is the time to onset of blindness in the left and right eyes of the patients with diabetic retinopathy. So, in this case, the focus is on the potential clustering in the data. Moreover, the data structure is that the observations from different clusters are independent, while observations within a cluster may be correlated. Interested readers can refer to the papers [1–4] and the references therein.In epidemiological research, the collection of detailed follow-up information is costly and time consuming, especially when the incidence of disease is low. One may consider the case-cohort design, proposed by Prentice [5], which is widely used for the large cohort studies. It entails collecting covariate data for all subjects who experienced the event of interest in the full cohort and for a random sample from the entire cohort. It has the same goals in studying risk factor as collecting the full data. Therefore, many studies about case-cohort data have been yielded. For example, it has been applied to the Cox proportional risk model [6], additive risk regression model [7], semiparametric transfer model [8], accelerated failure time model [9], and additive multiple risk regression model [10, 11].However, there is little research on group failure time data under case-cohort design. For example, Lu and Shih [12] applied the proportional risk model to clustered failure time data under case cohort and proposed three design methods of case cohort to extract subcolumns. Under different sampling mechanisms, the estimation equations are established, and then the large sample nature of the obtained estimation is proved. Unfortunately, the risk set considered by Lu and Shih [12] only contains the individuals in the subcolumns. Furthermore, Zhang et al. [13] added the failure information out of the subcolumns to the proportional risk model and proposed three kinds of different estimation equations and parameter estimates. The numerical simulation results show that the proposed model is effective. Note that the above literature studied the clustered data from case cohort based on the marginal proportional hazard models. In some practical cases, the effect of covariates may be additive. Then, the proportional risk model is no longer applicable. Therefore, in this paper, an additive risk model will be applied to clustered data in case-cohort design. Moreover, we will construct the weighted estimation equation, propose the estimation of regression parameters, and prove the asymptotic properties of the estimators.The remainder of the paper is organized as follows. In Section2, we describe the proposed estimation procedures. Then, we establish the consistency and asymptotic normality of the resulting estimator in Section 3. Numerical evidence in support of the theory is presented in Section 4. Finally in Section 5, some conclusions are drawn.


## 2. Model and Estimation Procedures

Suppose the full cohort consists ofn independent clusters, and the i-th cluster i=1,…,n has mi correlated subjects. We assume that subjects within the same cluster are exchangeable. In advance of follow-up, a random sample of the entire cohort, called the subcohort, is selected. Covariate data are then collected from individuals in the subcohort as well as those observed to fail in the entire cohort. Now, we consider the following three designs to obtain the subcohort [13]:Design A: randomly sample individuals from each cluster with Bernoulli sampling. In other words, each individual in each cluster has an independent fixed probability of being selected to the subcohort.Design B: randomly sample clusters from the full cohort with Bernoulli sampling.Design C: randomly sample clusters from the full cohort with Bernoulli sampling and then randomly sample subjects with Bernoulli sampling from the selected clusters.Note that Design A and Design B are special cases of Design C.LetTij and Cij denote the potential failure time and the potent censoring time for the j-th subject in the i-th cluster, respectively. Let Zijt represent a possibly time-dependent p×1 vector of covariates, restricted to be external. We assume that Tij and Cij are independent conditional on the given Zij⋅. The observed time is Xij=minYij,Cij. Let δij=ITij≤Cij be the indicator of failure, Nijt=IXij≤t,δij=1, and Yijt=IXij≥t.LetHi denote the indicator for the i-th subject being selected into the subcohort and Hij denote the indicator for the jth subject in the i-th cluster being selected into the sample. Both Hi and Hij are the independent Bernoulli variables, where ℰHi=γ,ℰHij=θ. Under Design A, Hi=1, for i=1,…,n. Under Design B, Hij=1, for i=1,…,n,j=1,…,mi.Suppose that the marginal hazard functionλijt is associated with Zijt as follows:(1)λijtZijt=λ0t+β0TZijt,where λ0⋅ is an unspecified baseline hazard function and β0 is a p×1 vector of regression parameters.If the full cohort data are available, the estimate of the true regression parameterβ0 in (1) could be obtained by solving the following estimating function [14]:(2)UHβ=∑i=1n∑j=1mi∫0τZijt−Z¯HtdNijt−YijtβTZijtdt,where Z¯Ht=∑i=1n∑j=1miYijtZijt/∑i=1n∑j=1miYijt.Then, the estimator ofβ0 in (1) can be estimated by β^H, which is the solution to the estimating equations UHβ=0. That is,(3)β^H=∑i=1n∑j=1mi∫0τYijtZijt−Z¯t⊗2−1∑i=1n∑j=1mi∫0τZijt−Z¯tdNijt,where τ<∞, a⊗2=aaT.For clustered failure time data from case-cohort studies, we consider the observed-event probability,Prδij=1, which we denote by p0 [15]. We propose three procedures to estimate β0, which are the same designs proposed by Zhang et al. [13] except that we consider the additive hazard regression model. We then develop a weighted estimating equation as follows:(4)Uβ,p=∑i=1n∑j=1miρij∫0τZijt−Z¯tdNijt−YijtβTZijtdt,where Z¯t=∑i=1n∑j=1miρijYijtZijt/∑i=1n∑j=1miρijYijt, ρij=p/N1/Nδij+1−p/n0/N1−δijHiHij with N=∑i=1nmi, N1=∑i=1n∑j=1miδij, and n0=∑i=1n∑j=1mi1−δHiHij.By solving the estimation equationUβ,p0=0, we could obtain β^ to estimate β0:(5)β^=∑i=1n∑j=1miρij∫0τZijt−Z¯t⊗2dt−1∑i=1n∑j=1miρij∫0τZijt−Z¯tdNijt.However,p0 is usually not known in most situations. In case where the study cohort is well defined, we can use p^w, the full cohort case proportion, to estimate p. When the study cohort is less well defined, p^s, the subcohort case proportion, is a suitable alternative. The corresponding estimator of β0 is written as β^w or β^s.The cumulative baseline hazard function,Λ0t=∫0tλ0udu, can be consistently estimated by(6)Λ^0t;β^,p^=∫0t∑i=1n∑j=1miρijdNiju−YijuβTZijudu∑i=1n∑j=1miρijYiju.In (4) and (6), either p^s, p^w, or p0 could be used to replace the parameter p^.


## 3. Asymptotic Properties

To derive the asymptotic properties of the proposed estimates, we impose the following assumptions:(C1)Ti,Ci,Zi⋅:i=1,…,n are independently and identically distributed, where Ti=Ti1,…,TimiT,Ci=Ci1,…,CimiT, and Zi=Zi1⋅,…,Zimi⋅T.(C2)PrYijτ>0>0 and Λ0τ<∞, for i=1,…,n,j=1,…,mi.(C3) The matrixℰW1β0,p0⊗2 is positive definite where(7)Wiβ,p=∑j=1mi∫0τZijt−etdMijt+1−δijγθHiHij−γθ∫0τZijt−etdMijt+D1βG1ip−G2ipp+D2βG1ip−G3ipγθ.(C4) There exists a constantCz such that Zij0+∫0τdZiju<Cz<∞ holds almost everywhere.Theorem 1.
Under the regularity conditions (C1)–(C4),n−1/2Uβ0,p0 converges to a mean zero normal distribution with covariance ∑β0,p0=ℰW1β0,p0⊗2, with(8)Wiβ,p=∑j=1mi∫0τZijt−etdMijt+1−δijγθHiHij−γθ∫0τZijt−etdMijt+D1βG1ip−G1ipp+D2βG1ip−G3ipγθ,D1β=ℰ∑j=1m1δ1jμ∫0τZ1jt−etdM1jt,D2β=ℰ∑j=1m11μγθ1−δ1jH1H1j∫0τZ1jt−etdM1jt,G1ip=n−1∑i=1nmi−μ,G2ip=n−1∑j=1miδij−μp,G3ip=n−1∑j=1mi1−δijHiHij−μγθ1−p,et=EY11tZ11tEZ11t.Theorem 2.
Under conditions (C1)–(C4),β^t converges in probability to β0, and each of n1/2β^t−β0 converges in distribution to a mean zero normal distribution with covariance matrix.Theorem 3.
Under conditions (C1)–(C4), bothβ^s and β^w converge in probability to β0, and each of n1/2β^s−β0 and n1/2β^w−β0 converges in distribution to a zero mean normal with covariance matrices Aβ0−1Ωsβ0,p0Aβ0−1 and Aβ0−1Ωwβ0,p0Aβ0−1, respectively, where for a=s or w, Ωaβ=ℰψ1β,p0⊗2, ψiβ,p=Wiβ,p+BβQiap, Ωaβ=ℰψ1β,p0⊗2, ψiβ,p=Wiβ,p+BβQiap, Qisp=1/μγθ∑j=1miHiHijδij−p,Qiwp=1/μ∑j=1miδij−p, Bβ=∫0τηijYijtZijt−Z¯t⊗2dt, and ηij=1/N1/Nδij−1/n0/N1−δijHiHij.Theorem 4.
Under the regularity conditions (C1)–(C4), for eachk=1,…,n, Λ^0β^,p^,t converges in probability to Λ0t uniformly in t∈0,τ. In addition, n1/2Λ^0β^,p^,t−Λ0t converges weakly to a Gaussian process with zero mean and covariance function at s,t given by ℰϕ1β0,p0,sϕ1β0,p0,t.


## 4. Numerical Studies

Simulation studies are conducted to examine the finite sample properties of the proposed estimators. In the study, we generate the clustered failure time data fromn=100 clusters. mi are simulated from a binomial (50, 0.8) distribution for i=1,…,n, with u=ℰmi=40. The total number in the whole queue is about 4000. The covariate Zij takes values 1 and 0, with probabilities 0.5 and 0.5, respectively. The failure time Ti1,Ti2,…,Timi in the i-th cluster is simulated by the joint survival function (for details, see [16]):(9)Sti1,…,timizi1,…,zimi=∑j=1miexp∫0tij2t+βTzijdtκ−mi−1−κ,where κ>0 is a parameter representing the degree of dependence among variables. Smaller value of κ represents stronger correlation.In our experiments, we chooseκ = 0.5 or 2 and β0=log0.5 or 0. The censoring times Cij are constant and equal to 1. The observed-event probabilities of p are equal to 0.45 or 0.63.For each data generation, for Design A, individuals within each cluster are selected into the subcohort by Bernoulli sampling with equal probability 0.2 or 0.15. For Design B, we select clusters by Bernoulli sampling with probability 0.2 or 0.15. For Design C, we first sample clusters by Bernoulli sampling with probability 0.4 or 0.3 and then sample individuals from those selected clusters by Bernoulli sampling with probability 0.5. Therefore, for each design, we would expect approximately 800 or 600 individuals in the subcohort.To assess the performance of the proposed estimator, we calculate the estimated standard error based on the bootstrap method. Because the failure clusters inside and outside the subcohort have different structures, we use the following method to conduct bootstrap sampling. Suppose the failure clusters in the subcohort contain M failures and Q nonfailures, whereas the failure clusters outside the subcohort contain D failures only. For each bootstrap sampling, we propose to obtain a bootstrap sample by sampling D failures from the original D failures outside the subcohort with replacement and then sampling M failures and Q nonfailures from the original M failures and Q nonfailures in the subcohort, respectively. The simulation results are based on 1000 replications.Tables1 and 2 report the simulation results of our proposed estimators. For simplicity, the notation Bias denotes the difference of the average of the estimates and the true values, SE denotes the average of standard errors, SEE denotes the average asymptotic standard error, BSE denotes the estimated standard error based on the bootstrap method, and CP denotes the coverage probability of the nominal 95% confidence interval.Table 1
Simulation results for the estimation ofβ0=log0.5.
DesignκMethodn=100,ns=800n=100,ns=600BiasSESEEBSECPBiasSESEEBSECPA0.5SC0.00100.05680.05740.05000.948−0.00210.06770.06350.06780.944WC−0.00260.05260.05760.05500.969−0.00230.06030.06460.06300.966T−0.00210.05330.05250.04960.953−0.00150.05940.05970.06420.9562SC0.00150.05730.05690.05060.948−0.00200.06680.06560.06720.942WC−0.00310.05310.05710.05600.953−0.00310.06250.06480.06020.953T−0.00160.05290.05360.05620.942−0.00100.05840.06030.06180.947B0.5SC−0.00350.05780.05690.04600.937−0.00120.06830.06440.07660.920WC0.00230.05930.06120.06090.9430.00040.06420.06400.06350.938T−0.00300.05340.05110.06050.934−0.00390.06290.06130.06140.9362SC−0.00120.05680.05730.05890.941−0.00090.06790.06540.06850.920WC0.00270.05780.05810.05130.9390.00150.06530.06370.05910.938T−0.00190.05340.05260.04800.940−0.00260.06290.06130.05570.936C0.5SC−0.00090.05760.05700.05510.948−0.00190.06850.06420.06510.929WC−0.00320.05430.05710.05710.958−0.00390.06050.06500.06320.963T−0.00150.05400.05180.04810.941−0.00470.06050.05910.06790.9482SC−0.00210.05760.05700.05150.948−0.00240.06790.06530.06140.934WC−0.00170.05470.05620.05650.942−0.00260.06140.06470.06440.956T−0.00100.05490.05230.05100.937−0.00380.05870.06040.06130.943Estimate ofβ0 from five methods: SC = estimating p0 using the subcohort p^s; WC = estimating p0 using the whole cohort, p^w; T = using true value, p0.Table 2
Simulation results for the estimation ofβ0=0.
DesignκMethodn=100,ns=800n=100,ns=600BiasSESEEBSECPBiasSESEEBSECPA0.5SC−0.00230.07300.07230.07480.945−0.00180.08390.08230.08410.950WC−0.00120.07390.07220.07170.9400.00010.08410.08230.07910.943T0.00020.07160.07230.07450.7460.00130.08250.08200.08680.9492SC−0.00100.07250.07360.06640.947−0.00220.08250.08310.08080.947WC−0.00170.07320.07290.06830.9320.00120.08370.08300.07090.937T0.00210.07210.07300.07100.9440.00310.08190.08270.08060.948B0.5SC−0.00130.07410.07140.07080.938−0.00570.08570.08000.06840.924WC−0.00020.07360.07160.07650.9310.00530.08550.07990.07830.925T0.00210.07340.07090.06190.927−0.00300.08420.08060.07560.9332SC−0.00230.07380.07250.06250.938−0.00420.08490.08130.08780.928WC−0.00150.07240.07210.06970.9340.00190.08580.08240.07740.930T0.00320.07410.07140.07250.931−0.00280.08360.08130.08330.940C0.5SC0.00050.07320.07150.07080.9490.00070.08430.08110.07460.937WC0.00020.07370.07180.07660.948−0.00220.08140.08580.08780.937T−0.00070.07070.07210.07440.948−0.00350.08090.08130.08360.9412SC0.00140.07390.07230.07260.9410.00140.08380.08200.08060.941WC0.00100.07410.07300.06850.938−0.00220.08230.08000.07630.935T−0.00220.07120.07310.07090.948−0.00150.08210.08170.07460.944Estimate ofβ0 from four methods: SC = estimating p0 using the subcohort p^s; WC = estimating p0 using the whole cohort, p^w; T = using true value, p0.It can be seen from Tables1 and 2 that the estimates of both regression parameters β0=log0.5 and β0=0 are unbiased, and both the estimated standard deviation and the empirical standard deviation are also close. The coverage probability of the nominal 95% confidence intervals seems to be very reasonable. Furthermore, when ns=600 for Design B in Tables 1 and 2, we can obtain that the values of SEE and CP are slightly underestimated. It is due to the small number of clusters in the subcohort. In addition, with the increase of the capacity of the subcolumns, i.e., when ns is increased to 800, the effect of the estimation is obviously improved. Among the above three designs, Design A and Design C are more effective than Design B. The conclusion is consistent with that in [13]. By comparing SE, SEE, and BSE, SE and SEE agree with the BSE. Furthermore, the bootstrap standard errors and the simulated standard errors are similar. The proposed semiparametric estimator works well. Therefore, one may use the bootstrap variance estimate for statistical inference in practice.


## 5. Conclusion

This paper proposed methods of fitting additive hazard regression models for clustered survival data from case-cohort studies. Risk difference can provide information value to medical research. Specifically, risk differences can provide information regarding the reduction in the number of cases developing a certain disease due to a decrease in a particular exposure. An advantage of the additive hazard models is that risk difference between different exposure groups can be readily derived from the coefficients in the additive hazard models. With respect to differing observed-event probability, we propose three procedures for parameter estimator. Consistency and asymptotic normality of the proposed estimators are established. The simulation results show that Design A results are more efficient estimators than Design C, and Design C has greater efficiency than Design B when subcohort size is approximately equal. It can be attributed to differences in the number of sampled clusters in the subcohort.



---

*Source: 9587870-2020-07-16.xml*