A new approach for semi-parametric regression analysis of bivariate interval-censored outcomes from case-cohort studies

Abstract Interval-censored failure time data frequently occur in many areas and a great deal of literature on their analyses has been established. In this article, we discuss the situation where one faces bivariate interval-censored data arising from case-cohort studies, which are commonly used as a tool to save costs when disease incidence is low and covariates are difficult to obtain. For this problem, a class of copula-based semi-parametric models is presented and for estimation, a sieve weighted maximum likelihood estimation procedure is developed. The resulting estimators of regression parameters are shown to be strongly consistent and asymptotically normal. Furthermore, the proposed method is generalized to the situation of non rare diseases. A simulation study is conducted to assess the finite sample performance of the proposed method and suggests that it performs well in practice.


Introduction
Bivariate interval-censored data occur when there exist two correlated failure times of interest and the observations on both are interval-censored (Sun 2006).A more specific example is given by a diabetic macular oedema (DMO) study in which, among others, one variable of interest is the time of reaching visual acuity of 70 Early Treatment Diabetic Retinopathy Study (ETDRS) letters.Due to the nature of the study, only interval-censored data on the the time are available.Many methods have been developed for bivariate interval-censored failure outcomes and in general, they can be classified into two types.One is the marginal approach that relies on the working independence assumption and ignores the correlation between the two related failure times of interest (Tong, Chen, and Sun 2008;Chen, Tong, and Zhu 2013).It is easy to see that one disadvantage of these methods is that they may be less efficient and also they cannot provide an estimation of the association of the correlated failure times.Corresponding to these, the other type of method is to directly model the correlation by using either frailty variables or copula models (Chen, Tong, and Sun 2009;Zhou, Hu, and Sun 2017;Cui, Zhao, and Sun 2018;Sun and Ding 2021).
Case-cohort design was first proposed by Prentice (1986) and discussed by many as a tool to reduce costs for large epidemiological cohort studies (Cai andZeng 2004, 2007).It collects the covariate information only on the subjects whose failures are observable as well as a sub-cohort of the rest of the subjects rather than on all study subjects.One adventure of the case-cohort design is that the same sub-cohort can be used for a variety of diseases so as to reduce costs further and to allow simultaneous evaluation of multivariate diseases.Many authors have discussed the analysis of interval-censored failure time data arising from casecohort studies but most of the available methods are for univariate situations (Gilbert et al. 2005;Li, Gilbert, and Nan 2008;Li and Nan 2011;Zhou, Zhou, and Cai 2017).
Several authors have discussed regression analysis of case-cohort multivariate failure time data but most of them only focused on right-censored data (Kang and Cai 2009;Kim, Cai, and Lu 2013) except Zhou, Cai, and Zhou (2021), who developed a marginal approach under the proportional hazards model.As mentioned before, the marginal approach may be less efficient and does not provide an estimation of the dependence, which may be required sometimes.To address this, in the following, we will propose a copula model-based approach.In comparison to marginal or frailty approaches, copula model-based methods directly link the two marginal distributions through a copula function to generate a joint distribution, whose copula parameter determines the dependence.This unique feature enables the modeling of the margins to be separated from the computation of the copula function, which is highly beneficial from a modeling perspective and an interpreting perspective.More specifically, we will present a class of semiparametric copula-based models for bivariate interval-censored data and the approach allows one to simultaneously estimate the association between the two event times and the covariate effects on their marginal distributions.
The remainder of this article is organized as follows.In Section 2, we begin by introducing some notation and assumptions that will be used throughout the article in addition to the structure of the observed data.We will then present a weighted likelihood estimation method in Section 3, which will employ the inverse probability weighting (IPW) technique to account for sampling bias.The asymptotic properties of the resulting estimators will be discussed in Section 4. Section 5 presents some results obtained from a simulation study conducted to examine the empirical performance of the proposed method and they suggest that it works well for practical situations.In Section 6, the proposed approach is illustrated by using a set of real data from the DMO study mentioned above, and Section 7 contains some concluding remarks.

Assumptions, background, and data structure
Consider a failure time study that involves two related failure times of interest denoted by T 1 and T 2 and suppose that there exists a vector of covariates Z = (Z 1 , Z 2 ) , where Z 1 and Z 2 may affect the two failure times T 1 and T 2 , respectively.Here Z 1 and Z 2 could be the same, completely different, or share some common components and they will be written as Z i1 and Z i2 for individual i.We denote the marginal survival function for T j from subject i as S j (t ij |Z ij ) = P(T ij > t ij |Z ij ), j = 1, 2, and the joint survival function for subject i as By the Sklar's theorem (Sklar 1959), so long as marginal survival functions S j are continuous, there exists a unique copula function Here the parameter α measures the dependence between the two margins.
As mentioned above, a significant feature of the copula model is that it allows for the dependence to be modeled separately from the marginal distributions (Nelsen 2006).One popular copula family for bivariate censored data is the Archimedean copula family, which usually has an explicit formula.In the following, we will focus on the one-parameter Archimedean copula model given by C α (u, v; α) = ψ −1 {ψ(u; α) + ψ(v; α); α} Here for a fixed α, ψ(•; α) : [0, 1] → [0, ∞] denotes a continuous, strictly decreasing convex function with ψ(1; α) = 0, and ψ −1 (•; α) its inverse function.Two special cases are the Clayton and Gumbel copula models, which account for the lower or upper tail dependence between two margins using a single parameter.Examples of various one-parameter Archimedean copula models are discussed in Nelsen (2006).
Note that although the inference for the association between the two event times can be made through the inference about the association parameter α, in practice, the Kendall rank correlation coefficient (Kendall's τ ) is usually used.This is because the range of α may depend on the copula function used, while Kendall's τ is a scale-invariant measure of the correlation and directly related to α through To model the covariate effect, we will assume that given Z j , the failure time T j marginally follows the proportional hazards model with the cumulative hazard function given by j (t|Z j ) = j (t) exp(β j Z j ).
( 1 ) In the above, j (t) denotes an unknown disease-specific increasing function, j = 1, 2, and β = (β 1 , β 2 ) is a p-dimensional vector of regression parameters.Note that here we assume that T 1 and T 2 may have different baseline cumulative hazard functions.
The case-cohort design can be considered as a two-phase sampling scheme.Suppose that there are n subjects in the study cohort.In the first phase, we observe interval-censored data on the two diseases of interest for all n cohort members.Suppose that the occurrence of the failure is not exactly observed but only determined at a sequence of examination times, denoted as 0 = U j0 < U j1 < • • • < U jM j ≤ ζ j , j = 1, 2. The failure time T j is then known to lie in the interval (L j , R j ], where L j = max{U jm : U jm < T j , m = 0, . . ., M j } and R j = min{U jm : U jm ≥ T j , m = 0, . . ., M j + 1} with U j,M j +1 = ∞.In particular, if the failure T j occurs before the first examination time U j1 , then (L j , R j ] = (0, U j1 ]; if the failure T j has not occurred at the last examination time, then (L j , R j ] = (U jM j , ∞].Denote the observed data for all cohort members from the first phase as In the following, we will assume that the joint distribution of the L ij 's, R ij 's and Z ij 's does not involve unknown parameters in model (1).Note that I(R ij < ∞) = 1 suggests that the ith subject has developed the jth disease.
In the second phase, we first select a simple random sample (SRS), called sub-cohort, by Bernoulli sampling with probability p s from the study cohort.Let ξ i indicate whether the ith subject is selected into the sub-cohort, i = 1, . . ., n.As part of the second phase, we measure the expensive covariates for all subjects in the sub-cohort (i.e., subjects with ξ i = 1) and for each member of the cohort with the disease (i.e., subjects with I(R ij < ∞) = 1 for j).Then the observed data can then be expressed as First-phase data: {O ij , ξ i : j = 1, 2, i = 1, . . ., n}; Second-phase data: Note that sometimes although a disease is non rare or not-so-rare, measuring the expensive covariates for every diseased individual in the cohort may still not be possible or necessary.In this case, the generalized case-cohort design can be used in which only a fraction of the diseased subjects is assessed for the expensive covariates (Kim, Cai, and Lu 2013).More specifically, it will select an SRS by Bernoulli sampling with probability p cj among the subjects who have the jth disease but are not included in the subcohort of the first phase.Let η ij denote whether the ith subject is selected into this sample, j = 1, 2, i = 1, . . ., n.Then the covariate measurements are gathered from the subjects with ξ i = 1 or η ij = 1, and the observed data can be summarized as First-phase data: {O ij , ξ i , η ij : j = 1, 2, i = 1, . . ., n}; Second-phase data: To see graphically, Figure 1 illustrates the data structures under the case-cohort design and the generalized case-cohort design, respectively.

Estimation and inference procedures
Now we discuss the estimation of or inference about the regression parameter β under both the case-cohort design and the generalized case-cohort design.Note that for both designs, the unknown covariates can be regarded as missing at random, and in order to account for missing information, a common way is to employ the inverse probability weighting (IPW) approach.This motivates us to consider the following inverse probability weighted loglikelihood function (2) In the above, the w i 's are some weights to be discussed below and the inverse of the sampling probability will be used as the basis for weighing observations if it is known.Also for a given subject i, if T ij is right-censored, then any term involving R ij becomes 0 (since R ij is set to be ∞) and the joint survival function for subject i reduces to either only one term (if both T i1 and T i2 are right-censored) or two terms (if one T ij is right-censored).
For the selection of the weights w i 's, we will use the information from both diseases and construct them in accordance with the proposed sampling schemes as shown in Figure 1.Specifically, the weights are defined as for the case-cohort design, and defined as for the generalized case-cohort design.As a result of the case-cohort design, all cases are assigned a weight of 1 and non cases in the subcohort are assigned a weight of 1/p s , where p s is the sampling probability.For the generalized case-cohort design, all cases selected at the second phase will be weighed by ), the cases in the subcohort will be given the weigh of 1, and the non cases will be weighed by 1/p s , where p cj is the sampling probability of cases in the sub-cohort for the disease j, j = 1, 2. As for subjects with no observed covariable, we have w i = 0 in both designs.
As one can see, the unknown parameters in the weighted log-likelihood function (2) include both regression coefficients and the cumulative baseline hazard functions.To address this, by following Zhou, Hu, and Sun (2017) and others, we employ a sieve method with Bernstein polynomials to handle the unknown cumulative baseline functions { 1 , 2 }.In particular, let where p denotes the dimension of β, M a positive constant, ||v|| denotes the Euclidean norm for a vector v, and M j is the collection of all bounded and continuous non decreasing, non negative functions over the interval [c j , u j ] with 0 ≤ c j < u j < ∞, j = 1, 2. We define the sieve parameter space as Here with B k (t, m, c j , u j ) being Bernstein basis polynomials defined as with the degree m = o(n v ) for some v ∈ (0, 1), j = 1, 2. In practice, [c j , u j ] is usually taken as the range of the finite values of the observed L ij 's and R ij 's.Then it is natural to estimate θ by the sieve maximum likelihood estimator θn = ( βn , αn , ˆ 1n , ˆ 2n ) of θ defined to be the value of θ that maximizes the weighted log-likelihood function over n .
For the implementation of the estimation approach proposed above, it should be noted that there are some restrictions on the parameters due to non negativity and monotonicity of the functions 1 and 2 .On the other hand, they can be easily removed through reparameterization.An ideal method is to re-parameterize the parameters {φ j0 , . . ., φ jm }, j = 1, 2 as the cumulative sums of {exp(φ * j0 ), . . ., exp(φ * jm )}.For the determination of θn , an iterative algorithm is clearly needed.For the selection of initial estimates, one can first use the marginal method (Zhou, Cai, and Zhou 2021) to get the initial values of (β, 1 , 2 ), denoted as ( β(1) , ˆ (1)  1 , ˆ (1) 2 ).Then the initial values of α, denoted as α(1) , can be chosen by maximizing the weighted log-likelihood (2) with substituting ( β(1) , ˆ (1)  1 , 1 , ˆ (1) 2 |Z).One can simultaneously maximize (2) to get final estimates θn = arg max θ (θ|Z) with initial values ( β(1) , α(1) , ˆ (1)  1 , ˆ (1) 2 ).A further issue is the choice of the degree of Bernstein polynomials m, which dictates the smoothness of the approximation.For this, we suggest considering several different values of m and then applying the AIC that selects the value of m that minimizes In the numerical studies below, we employ the Quasi-Newton algorithm built in fminunc in Matlab.After getting the estimator θn , one can obtain a natural estimator of the marginal survival Ŝj (t|Z j ) = exp{− ˆ j (t) exp( β j Z j )}, j = 1, 2 by substituting j and β j to the survival function of (1).Furthermore, an estimator for the joint survival function

Asymptotic properties
In this section, we establish the large sample properties of θn .Note that all weights are bounded and do not depend on θ, and they satisfy E{w i |O ij , j = 1, 2} = 1.Let ||f || ∞ = sup t |f (t)|, the supremum norm for a function f .Furthermore, for any where v ∈ (0, 1) such that m = o(n v ) and d is defined in Condition 3. Theorem 4.3.(Asymptotic normality).Assume that Conditions 1 -5 given in Appendix hold.Then we have that in distribution, where is given in Appendix.
For inference about (β , α) , it is obvious that we need to estimate the covariance matrix of their estimators.For this, one approach is to derive a consistent estimator of but it can be seen from Appendix that this would not be straightforward.To provide a more stable variance estimator, we suggest to apply the following simple weighted bootstrap procedure by following others (Ma and Kosorok 2005), which is easy to implement and works well in the numerical study.In particular, the weighted bootstrap procedure is suitable for the missing information and the IPW likelihood by reweighting each individual to achieve the effect of resampling.More specifically, let {u 1 , . . ., u n } denote n independent realizations of a bounded positive random variable u satisfying E(u) = 1 and Var(u) = 0 < ∞.Define the new weights

A simulation study
To examine the empirical performance of the sieve maximum likelihood estimation procedure proposed in the previous sections, a simulation study was conducted.In this study, we considered two diseases and generated their failure times T 1 and T 2 based on the joint survival function S(t 1 , t 2 |Z) = C α (S 1 (t 1 |Z), S 2 (t 2 |Z); α) and model (1) with 1 (t) = log(1 + t/10) and 2 (t) = log(1 + t/8).It was assumed that both covariates Z 1 and Z 2 are one-dimensional and follow the Bernoulli distribution with the success probability 0.5 and the standard normal distribution N(0, 1), respectively.For the copula model, we considered That is, they follow either Clayton or Gumbel copula model.For the two models, we have that α = 2τ/(1 − τ ) or α = 1/(1 − τ ), for the relationship between α and τ .
For the generation of the censoring intervals or observed data, for each subject, we first generated a sequence of observation time U j1 < • • • < U jM j , j = 1, 2 over the study period [0, ζ ] such that we keep sampling the uniform random variables with the mean 0.05I(|Z| > 1) + 0.1I(|Z| ≤ 1) until the cumulative sum of these independent and identically distributed random variables is greater than ζ .Note that here the observation process was assumed to depend on the covariate to be more general.Then as above, we defined L j = max{U jl : U jl < T j , l = 0, . . ., M j } and R j = min{U jl : U jl > T j , l = 1, . . ., M j + 1}, where U j0 = 0 and U j,M j +1 = ∞.For ζ , we considered two values, 1.25 and 2.25, yielding the proportions (p d1 , p d2 ) of the cases with respect to T 1 and T 2 in the cohort being close to (0.08, 0.14) or (0.14, 0.22), respectively.
For the generation of the sub-cohort under the case-cohort design, we used the independent Bernoulli sampling with the success rate p s = 0.1.Under the generalized case-cohort design, we further generated a sub-cohort of all cases by using the independent Bernoulli sampling with the success rate (p c1 , p c2 ) = (0.5, 0.5).The results given below are based on the cohort size n = 2000 or 3000, the degree of Bernstein polynomials being m = 3, and the size of the weighted bootstrap sample being 50 with 1000 replications.Table 1 presents the results on estimation of β and τ given by the proposed estimation procedures with the true value of β being (−0.5, 0.5) and ζ = 1.25 for the case-cohort design or ζ = 2.25 for the generalized case-cohort design.Here for Kendall's τ , we considered two values τ = 0.3 and 0.6, representing the weak and strong associations between T 1 and T 2 , respectively.The results include the empirical bias (Bias) given by the average of the proposed estimates minus the true value, the sample standard error of the estimates (SSE), the mean of the estimated standard errors (ESE), and the 95% empirical coverage probability (CP).
One can see from Table 1 that the proposed estimators seem to be unbiased and the bootstrap variance estimation procedure appears to work well.In addition, the results on the empirical coverage probabilities suggested that the normal approximation to the distributions of the estimators appears to be appropriate.As expected, the results became better when the sample size increased.To see the performance of the proposed method on estimation of the cumulative hazard function, Figures 2 and 3 show the averages of ˆ 1n and ˆ 2n for the casecohort design situations corresponding to Table 1 with τ = 0.6 and n = 2000 under Clayton copula and Gumbel copula, respectively.They indicate that the Bernstein polynomial-based method seems to work well.We also considered other copulas such as the FGM and Frank copulas and obtained similar conclusions.
We also carried out a simulation study by considering both Z 1 and Z 2 being twodimensional and assuming that the marginal cumulative hazard function of T j follows the proportional hazards model   1 under the casecohort design with Gumbel copula, τ = 0.6 and n = 2000 (left: 1 (t); right: 2 (t)): the true function (solid), the estimated function (dashed), and the 95% confidence bands (dash-dot).j = 1, 2. In the above, 1 (t) and 2 (t) are defined as above and it was assumed that Z1 ∼ N(0, 1), Z12 ∼ Ber(0.5) and Z22 ∼ Ber(0.5) with Z12 and Z22 being independent.We simulated interval-censored failure time data in the same way as before.Also the uniform random variable was used to generate the interval-censored failure time data, dependent on the covariates ( Z1 , Zj2 ) by setting its mean equal to 0.05I(| Z1 | > 1 and Zj2 = 1) + 0.1I(| Z1 | ≤ 1 or Zj2 = 0).Note that here in addition to having different covariates, the two failure times T 1 and T 2 also have one shared covariate Z1 .Table 2 presents the results on estimation of β and τ given by the proposed estimation procedures with n = 2000 and the true values of the regression parameters are (β 11 , β 12 ) = (−0.5,0.3) and (β 21 , β 22 ) = (0.5, 0.7).Here we set ζ = 0.65 for the case-cohort design and ζ = 1.20 and (p c1 , p c2 ) = (0.5, 0.5) for the generalized case-cohort design.They correspond to the proportions of cases in the cohort (p d1 , p d2 ) being close (0.07, 0.10) or (0.12, 0.20) for T 1 and T 2 , respectively.The proposed method yielded quite similar performance as above.
As discussed above, one advantage of the proposed copula model-based method over a marginal approach is that the former can be more efficient than the latter.To investigate or verify this, we applied the marginal estimation procedure that maximizes the following inverse probability weighted log-likelihood function obtained under the working independence assumption between T 1 and T 2 , to the simulated data giving Table 1 with n = 2000.For the variance estimation, as above, a similar weighted bootstrap procedure was employed, which is consistent with the argument in Zhou, Cai, and Zhou (2021).The obtained results on estimation of β and τ are given in Table 3 along with the results given by the proposed method.It is apparent that the proposed estimator is clearly more efficient than the marginal estimator.In particular, as expected, the relative efficiency gain of the proposed estimator over the marginal method increases as Kendall's τ increases.In the proposed estimation procedures, it has been assumed that the copula function is known and it is apparent that this may not be true in reality.Thus, it will be interesting to investigate how this may affect the proposed estimators.To see this, we repeated the study above given in Table 1 except that it was assumed that 1  4. Note that in the table, by Misspecified model, we mean that for the data generated under the Clayton model, the Gumbel model was assumed to be the true one, while for the data generated under the Gumbel model, the Clayton model was assumed to be the true one.It seems that the results obtained under the misspecified model are similar to those obtained under the correct one.In other words, the proposed estimation approach appears to be robust with respect to the copula model.

An illustration
In this section, we illustrate the sieve maximum likelihood estimation approach proposed in the previous sections by using a set of real data arising from the DMO study described above (Kern et al. 2021).It is a longitudinal cohort study containing electronic medical records for all patients undergoing intravitreal injections in a tertiary referral center between March 2013 and October 2018.Treatment response in terms of visual acuity (VA) outcomes was reported for all eyes over the observation period.Among others, one objective is to evaluate the VA outcomes of intravitreal anti-vascular endothelial growth factor (VEGF) in DMO, and one particular VA outcome of interest is the Early Treatment Diabetic Retinopathy Study (ETDRS) letters score.The cohort includes 2614 DMO eyes of 1964 patients over 48 months.In the following, we will focus on the 650 subjects who gave full records of both left and right eyes.The failure time of interest here is the time to absolute VA attaining 70 ETDRS letters or above.Due to the nature of the study, the occurrence time is only known to be between the last observation time when the VA is below 70 and the first observation time when it had already attained.In other words, only interval-censored data are available for both eyes.
In the previous analysis, Kern et al. (2021) simplified the observed data into right-censored data and concluded that being male and the VA at baseline were positively associated with VA ≥ 70 and the age at baseline was correlated inversely with VA ≥ 70.In addition, they ignored the association between both eyes from one person.For the analysis below, by following Kern et al. (2021), we will focus on three time-independent covariates for left (L) and right (R) eyes: gender, the age at baseline, and the VA at baseline.In other words, we set Z 1 and Z 2 to the same and both are three-dimensional covariate vectors.To apply the proposed estimation procedures, let T 1 and T 2 denote the first time that the left eye and right eye attained 70 ETDRS, respectively.Under the case-cohort design, a sub-cohort was selected by using the independent Bernoulli sampling with a success rate of 0.2.With the generalized case-cohort design, we further sampled some cases with both eyes outside of the sub-cohort above using the independent Bernoulli sampling with a success rate of 0.5.
Table 5 presents the results given by the proposed methods under either the Clayton or Gumbel model assumption with m = 5 and 100 bootstrap samples for variance estimation.They include the estimated Kendall's τ and covariate effects (Est) along with the estimated error (SE) and the p-value for testing the regression parameter being zero.For comparison, we also obtained the results based on the full cohort, assuming that the covariate information is available for all subjects in the proposed method, and by applying the marginal approach given in Zhou, Cai, and Zhou (2021).It can be seen that under either Clayton or Gumble model, the full cohort method, case-cohort method, and generalized case-cohort method gave quite similar results.In particular, they indicated that there existed a significant dependence between left eye and right eye.
Although the results obtained under the two different copula models differ, they yielded similar conclusions.In particular, they suggested that the times of attaining 70 ETDRS letters for both eyes were significantly positively correlated with the VA at baseline and negatively correlated with the age at baseline.In addition, gender did not seem to have any effect on both conversions, and the estimated effects of the age at baseline and the VA at baseline are similar to the conclusions given in Kern et al. (2021).Also as seen in the simulation study, the marginal method is less efficient on most of the estimators than the proposed copula modelbased methods.We also tried different values for m and B and obtained similar results.

Concluding remarks
In the preceding sections, we have discussed regression analysis of bivariate interval-censored data arising from either case-cohort or generalized case-cohort studies.For the problem, copula model-based approaches were proposed and in the methods, sieve maximum likelihood estimation procedures based on Bernstein polynomials were developed.The use of the sieve approach makes the proposed method easily implemented and computationally fast.The proposed estimators of regression parameters were shown to be strongly consistent and asymptotically normally distributed.Compared to the marginal approach, the proposed method directly models the dependency between the correlated failure times and can yield more efficient estimators.In addition, a weighted bootstrap procedure was given to estimate the covariance matrices of the proposed estimators of regression parameters and the numerical study indicated that the proposed methodology works well for practical situations.
There exist several directions related to the methods proposed above for future research.One is that in the preceding sections, it has been assumed that the failure times of interest marginally follow the proportional hazards model and it is apparent that sometimes it may not fit data well or one may prefer some other models such as the additive hazards model or the semi-parametric transformation model.Thus, it would be useful to generalize the proposed method to alternative regression models.
Another assumption used above is that the joint survival function can be described by a known one-parameter Archimedean copula model and it is easy to see that this may be questionable in practice.In other words, one may want to develop some model checking procedures or develop similar estimation procedures under other copula models such as two-parameter Archimedean copula models.A third possible direction for future research is that the focus in the above has been on the situation where the observation processes are independent of the failure times of interest or the censoring is non informative.It is clear that this could be violated sometimes and the proposed method would not be applicable anymore.In other words, it would be useful to develop new estimation procedures that allow for dependent or informative censoring (Sun 2006).

A. Technical details
Let O ψ = L j , R j , ψZ, ψ; j = 1, 2 denote a single observation, where ψ indicates whether Z is observed.For the proofs, let P n denote the empirical measure based on n independent observations and P denote the true probability measure.The letter C represents a constant, and it does not necessarily represent the same value each time here and in the proofs.We define the covering number of the class L n = { (θ , O ψ ) : θ ∈ n }, where (θ , O ψ ) is the weighted log-likelihood function based on a single observation O ψ .In particular, (θ , O ψ ) is given by where w is the weight, O = {L j , R j , Z; j = 1, 2} is the complete observation.To establish the asymptotic properties, we need the following regularity conditions.Condition 1. (i) (β 0 , α 0 ) is an interior point of B ⊂ R p × R. (ii) The distribution of Z has bounded support and is not concentrated on any proper subspace of R p .

Condition 2. (i)
There exists τ j > 0 such that P(R j − L j ≥ τ j ) = 1, j = 1, 2; (ii) The union of the supports for the distributions of the L ij 's and the finite values of R ij 's is an interval [c, u] with 0 ≤ c j < c < u < u j < ∞.
Condition 3. The function j0 is continuously differentiable up to order d in [c, u]  for each of these p + 3 values, then v 1 = 0 p×1 and v 2 = v 3 = v 4 = 0.
Note that the conditions above are mild in that they are usually satisfied in practical situations.In particular, Conditions 1, 2, 3, and 5 are commonly used in the studies of interval-censored data

Figure 1 .
Figure 1.Graph illustration of case-cohort design (left) and generalized case-cohort design (right).
10 , 20 ) denote the true value of θ.The following theorems give the asymptotic consistency and normality of the proposed estimators when n → ∞.Theorem 4.1.(Strong consistency).Assume that Conditions 1 -4 given in Appendix hold.Then we have that d( θn , θ 0 ) converges almost surely towards 0. Theorem 4.2.(Rate of convergence).Assume that Conditions 1 -5 given in Appendix hold.Then we have that d

Table 1 .
Simulation results on estimation of β and τ .

Table 2 .
Simulation results on estimation of β and τ with two covariates.

Table 3 .
Comparison between the proposed approach and marginal approach with n = 2000.

Table 4 .
Simulation results on estimation of β under misspecified copula models.

Table 5 .
Analysis results for the DMO data.