# The Adjustment of Covariates in Cox’s Model under Case-Cohort Design

Authors: Guocai Rong; Luwei Tang; Wenting Luo; Qing Li; Lifeng Deng
Journal: Complexity
 (2020)

---

## Abstract

Case-cohort design is a biased sampling method. Due to its cost-effective and theoretical significance, this design has extensive application value in many large cohort studies. The case-cohort data includes a subcohort sampled randomly from the entire cohort and all the failed subjects outside the subcohort. In this paper, the adjustment for the distorted covariates is considered to case-cohort data in Cox’s model. According to the existing adjustable methods of distorted covariates for linear and nonlinear models, we propose estimating the distorting functions by nonparametrically regressing the distorted covariates on the distorting factors; then, the estimators for the parameters are obtained using the estimated covariates. The proof of consistency and being asymptotically normal is completed. For calculating the maximum likelihood estimates of the regression coefficients subject in Cox’s model, a minorization-maximization (MM) algorithm is developed. Simulation studies are performed to compare the estimations with the covariates undistorted, distorted, and adjusted to illustrate the proposed methods.


---

## Body

## 1. Introduction

Many survival data have the characteristics of large sample size and high censoring rate. The cost is very high when all variables of each individual in the data are measured. To reduce the cost and improve the efficiency of cohort studies, Prentice [1] firstly proposed the case-cohort design and gave a pseudolikelihood method to estimate regression parameters. Since the publication of the landmark article [1], case-cohort design has been applied more and more, especially with the development of big data in recent years. For example, some scholars have applied this design to Life Science [2–5] and natural disasters [6], and some other scholars have further improved this design [7, 8], respectively.In practice, some subjects may be interfered by other factors due to their characteristics, so their corresponding real covariates cannot be observed; only the distorted covariates and the distorting factors can be observed. For example, for the Modification of Diet in Renal Disease (MDRD) study [9, 10], glomerular filtration rate (GFR) and serum creatinine (SCr) data are distorted by the body surface area (BSA). To show the relationship of the two covariates, we need to adjust them by correcting the distorting effect of BSA. The covariate-adjusted regression was introduced for situations where both predictors and response in a regression model are not directly observable, but are contaminated with a multiplicative factor that is determined by the value of an observable factor [11]. For the nonlinear regression model, how to estimate the distorting functions was proposed by nonparametric regression the predictors and response on the distorting covariates [12]. Then, the covariate-adjusted method was extended to other models [13–15]. The data studied in these documents is all complete data. In this paper, the studied data is survival data with right censored data, in which only distorted covariates and distorting factors are included, and the nonparametric method [12] to estimate the distorting functions to obtain the adjusted value of the distorted covariates on the distorting factors, where the nonparametric method is a kernel smoothing method.In this article, inspired by the idea of Ding et al. [16] and Deng et al. [17], we construct a surrogate function with a diagonal Hessian matrix by using the convexity of the exponential function and the negative logarithm function and then maximize this surrogate function with a diagonal Hessian matrix. This algorithm is minorization-maximization (MM) algorithm. That is, the first M is to construct a mi-surrogate function and the second M is to maximize the function.The rest of the article is organized as follows. In Section2, we firstly fit data from case-cohort design to Cox’s model with distorting factors and secondly adjust the distorted covariates by using the kernel function. In Section 3, the convergence of the adjusted covariates and the asymptotic properties for the proposed estimators is completed. In Section 4, we propose a MM algorithm for implementation of the estimation and present a cross-validation (CV) to obtain the optimal bandwidth and a nonparametric bootstrap approach to get the standard error estimation. Several simulation studies are conducted to compare the estimations without distorting factors, with distorting factors, and adjusted with distorting factors in Section 5. Real data analysis of heart failure data is in Section 6. Discussion is stated in Section 7. All proofs are given in Appendix.


## 2. Design and Estimation

### 2.1. Model and Design

Suppose that there aren independent subjects in the studies cohort, where T˜i is the failure time, Ci is the censoring time, Ti=minT˜i,Ci is the observe time for the ith subject, Δi=IT˜i≤Ci is the right-censoring indicator for failure, Yit=ITi≥t is the at risk process, and Nit=ΔiITi≤t is the counting process, where I⋅ is an indicator function, τ is the end of study time, and Zi is a pth covariate for subject i; we focus on the time-independent covariate.Letβ=β1,β2,…,βp be the unknown p-dimensional vector of regression coefficients. T˜i arises from Cox’s model as the following form:(1)λt|Z=λ0texpβ⊤Zit,where λ0t is an unspecified baseline hazard function. The corresponding partial likelihood function is widely used for the inference of β [18] as the following form:(2)LFβ=∏i=1nexpZiTiβ⊤∑l=1nYlTiexpZiTiβ⊤Δi.Under the case-cohort design, letA˜ denote the subcohort, which is selected from the full cohort by simple random sampling, and ξi be an indicator, equaling 1 if the ith subject is selected into the subcohort and 0 otherwise. Let A denote the case-cohort sample which contained the subjects from the subcohort A˜ and the case outside the subcohort. Therefore, the observed data structure can be summarized as follows: Ti,Δi,ξi+Δi1−ξiZit,0≤Ti.The pseudolikelihood function [1] was proposed, and the corresponding pseudolikelihood function takes the following form:(3)Lpβ=∑i=1nωiΔiβ⊤ZiTi−log∑i=1nYlTiωlexpβ⊤ZlTi,where ωi=Δi+1−Δiξi.


### 2.2. The Adjustment of Covariates

We now study that the covariates are distorted by some distorting factors. The distorted observed data structure is(4)Ti,Δi,ξi+Δi1−ξiZ˜it,Ui,0≤t≤Ti,where Ui is the ith distorting factor and Z˜it is the ith distorted covariate. Here, Zit is unobservable, and ϕU is denoted as the unknown distorting functions of observable factor U.Firstly, we give some basic assumptions [12] as follows:(i)
Z˜r=ϕrUZr(ii)
EϕrU=1,r=1,2,…,p(iii)
Zr,U are independent of each other, r=1,2,…,pUnder these assumptions, the working-likelihood function (3) can be rewritten as the following form:(5)l˜pβ=∑i=1nωiΔiβ⊤Z˜iTi−log∑l=1nYlTiωiexpβ⊤Z˜lTi,Z˜ri=ϕrUiZri.Now, our object is to estimate the unknown coefficientsβ based on observation (4) and function (5). From condition (A3), we have(6)ϕrU=EZ˜r|UEZr,1≤r≤p.Many scholars have exposed some methods to estimateϕrU [19]. So, we construct the function ϕrU by using the classic kernel method, for 1≤r≤p:(7)ϕ^rU=1/nh∑i=1nKu−Ui/hZ˜ir1/nh∑i=1nKu−ωi/h×1Z˜¯r,where Z˜¯r=1/n∑i=1nZ˜ri, K⋅ is a kernel function, and h is bandwidth.LetZ^ri=Z˜ri/ϕ^rUi denote the adjusted covariate. That is, we obtain the adjusted covariates by removing the distorting factors from the distorted covariates. Then, function (4) can be abbreviated as(8)l^pβ=∑i=1nωiΔiZ^iTiβ⊤−log∑l=1nYlTiωlexpZ^lTiβ⊤.Under the case-cohort design, the estimateβ^ about the adjustment of covariates is defined as follows:(9)β^=argmaxβl^pβ.


## 2.1. Model and Design

Suppose that there aren independent subjects in the studies cohort, where T˜i is the failure time, Ci is the censoring time, Ti=minT˜i,Ci is the observe time for the ith subject, Δi=IT˜i≤Ci is the right-censoring indicator for failure, Yit=ITi≥t is the at risk process, and Nit=ΔiITi≤t is the counting process, where I⋅ is an indicator function, τ is the end of study time, and Zi is a pth covariate for subject i; we focus on the time-independent covariate.Letβ=β1,β2,…,βp be the unknown p-dimensional vector of regression coefficients. T˜i arises from Cox’s model as the following form:(1)λt|Z=λ0texpβ⊤Zit,where λ0t is an unspecified baseline hazard function. The corresponding partial likelihood function is widely used for the inference of β [18] as the following form:(2)LFβ=∏i=1nexpZiTiβ⊤∑l=1nYlTiexpZiTiβ⊤Δi.Under the case-cohort design, letA˜ denote the subcohort, which is selected from the full cohort by simple random sampling, and ξi be an indicator, equaling 1 if the ith subject is selected into the subcohort and 0 otherwise. Let A denote the case-cohort sample which contained the subjects from the subcohort A˜ and the case outside the subcohort. Therefore, the observed data structure can be summarized as follows: Ti,Δi,ξi+Δi1−ξiZit,0≤Ti.The pseudolikelihood function [1] was proposed, and the corresponding pseudolikelihood function takes the following form:(3)Lpβ=∑i=1nωiΔiβ⊤ZiTi−log∑i=1nYlTiωlexpβ⊤ZlTi,where ωi=Δi+1−Δiξi.


## 2.2. The Adjustment of Covariates

We now study that the covariates are distorted by some distorting factors. The distorted observed data structure is(4)Ti,Δi,ξi+Δi1−ξiZ˜it,Ui,0≤t≤Ti,where Ui is the ith distorting factor and Z˜it is the ith distorted covariate. Here, Zit is unobservable, and ϕU is denoted as the unknown distorting functions of observable factor U.Firstly, we give some basic assumptions [12] as follows:(i)
Z˜r=ϕrUZr(ii)
EϕrU=1,r=1,2,…,p(iii)
Zr,U are independent of each other, r=1,2,…,pUnder these assumptions, the working-likelihood function (3) can be rewritten as the following form:(5)l˜pβ=∑i=1nωiΔiβ⊤Z˜iTi−log∑l=1nYlTiωiexpβ⊤Z˜lTi,Z˜ri=ϕrUiZri.Now, our object is to estimate the unknown coefficientsβ based on observation (4) and function (5). From condition (A3), we have(6)ϕrU=EZ˜r|UEZr,1≤r≤p.Many scholars have exposed some methods to estimateϕrU [19]. So, we construct the function ϕrU by using the classic kernel method, for 1≤r≤p:(7)ϕ^rU=1/nh∑i=1nKu−Ui/hZ˜ir1/nh∑i=1nKu−ωi/h×1Z˜¯r,where Z˜¯r=1/n∑i=1nZ˜ri, K⋅ is a kernel function, and h is bandwidth.LetZ^ri=Z˜ri/ϕ^rUi denote the adjusted covariate. That is, we obtain the adjusted covariates by removing the distorting factors from the distorted covariates. Then, function (4) can be abbreviated as(8)l^pβ=∑i=1nωiΔiZ^iTiβ⊤−log∑l=1nYlTiωlexpZ^lTiβ⊤.Under the case-cohort design, the estimateβ^ about the adjustment of covariates is defined as follows:(9)β^=argmaxβl^pβ.


## 3. Consistency and Asymptotic Normality

To present the asymptotic results, we will introduce some notations. Letβ0 denote the true value of β. For d=0,1,2, define(10)Sdβ,t=1n∑i∈AYitexpZiTβZi⊗d,S˜dβ,t=1n˜∑i∈A˜YitexpZi⊤βZi⊗d,Qdβ,t,ω=1n∑i∈AYitYiωexp2Zi⊤βZi⊗d,Q˜dβ,t,ω=1n˜∑i∈A˜YitYiωexp2Zi⊤βZi⊗d,where a⊗0=1,a⊗1=a,anda⊗2=aa⊤ for a vector a.Define(11)Eβ,t=S1β,tS0β,t,E˜β,t=S˜1β,tS˜0β,t.We impose the following conditions throughout the paper. Note that convergence properties involven⟶∞ and ⋅ refers to the Euclidean norm in these conditions:(A1)
∫0τλ0tdt<∞.(A2)
There exists someδ>0 such that n−1/2supi,t∈0,τYitZiIZiTβ0>−δZi⟶P0.(A3)
There exists a neighborhoodβ of β0, and functions sdβ,t, d=0,1,2, defined on β×0,τ, which satisfy(i)
supβ∈ℬ,t∈0,τSdβ,t−sdβ,t⟶P0,d=0,1,2.(ii)
sdβ,t, d=0,1,2, are continuous of β on ℬ uniformly in t∈0,τ; sdβ,t, d=0,1,2, are bounded on ℬ×0,τ, and s0β,t is bounded away from zero; and s1β,t=∇βs0β,t and s2β,t=∇β2s0β,t for all β∈ℬ, t∈0,τ.(iii)
The matrix(12)Σβ0=∫0τνβ0,ts0β,tλ0tdtis positive definite, whereνβ0,t=s2β0,t/s0β0,t−s1β0,t/s0β0,t⊗2.(B1)
n˜/n⟶α for some α∈0,1.(B2)
The sequence of distributions ofn1/2E˜β0,t−Eβ0,t is tight on ℬ×0,τ of left-continuous functions with right-hand limits equipped with the product Skorohod topology.(B3)
The existing functionsqdβ,t,ω, d=0,1,2, defined on ℬ×0,τ2 which satisfy(i)
supβ∈ℬ,t,w∈0,τ2Qdβ,t,ω−qdβ,t,ω⟶P0,d=0,1,2.(ii)
qdβ,t,ω, d=0,1,2, are continuous of β on ℬ uniformly in t,w∈0,τ2, and qdβ,t,ω, d=0,1,2, are bounded on ℬ×0,τ2.(iii)
supn≥1EQdβ,t,ω, d=0,1,2, are bounded sequences.(B4)
Ford=0,1,2,(13)supβ∈ℬ,t∈0,τS˜dβ,t−sdβ,t⟶P0,supβ∈ℬ,t,w∈0,τ2Q˜dβ,t,ω−qdβ,t,ω⟶P0.(C)
AllgrU=ϕrUpU,1≤r≤p, and ϕrU and pU are greater than a positive constant. They are differential and their derivatives satisfy the condition that there exists a neighborhood of the origin. For example, Δ and a constant c>0 such that, for any ξ∈Δgr3U+ξ−gr3U≤cξ,1≤r≤p, p3U+ξ−p3U≤cξ.(D1)
The continuous kernel functionK⋅ satisfies the following properties:(i)
The support ofK⋅ is the interval [−1, 1].(ii)
K⋅ is symmetric about zero.(iii)
∫−11KUdU=1 and ∫−11UdkUdU=0,d=1,2,3.(D2)
Asn⟶∞, the bandwidth h is in the range from On−1/4logn to On−1/8.(E)
EZr is bounded away from 0, EZr2<∞. These conditions are mild and can be satisfied in most circumstances.Conditions (A1)–(A3) are the regularity conditions for the asymptotic results of Cox’s model [20]. Conditions (B1)–(B4) are the regular conditions under the case-cohort design [21]. Conditions (C) are related to smoothness of the function gr⋅ and the density function p⋅ of U. Conditions (D1)-(D2) are commonly assumed for the root n consistency of kernel-based estimation [22, 23]. Condition (E) is special for this problem [12, 22].Under these conditions, we have the following result with detailed proofs given in Appendix.Theorem 1.
(asymptotic properties ofβ^). Under conditions (A1)–(A3), (B1)–(B4), (C), (D1)-(D2), and (E), β^ converges in probability to β0, and(14)nβ^−β0⟶dN0,Ωβ0,where the asymptotic variance matrix Ωβ0=Σ0β0+Σ1β0, where(15)Σ1β0=∫0τ∫0τΨβ0,t,ωs0β0,ts0β0,ωλ0tλ0ωdtdω,Ψβ0,t,w=1−αα−1s0β0,ts0β0,w−1⋅q2β0,t,w−q1β0,w,teβ0,w′−eβ0,tq1β0,t,w′+q0β0,t,weβ0,teβ0,w′,eβ0,t=s1β0,ts0β0,t,eβ0,ω=s1β0,ωs0β0,ω.


## 4. The Implementation of Algorithm

InSection 3, we have finished the asymptotic properties for the adjusted estimation β^. To obtain β^ by solution of function (9), we can use the Newton–Raphson algorithm. However, in the process of calculation, we find that its Hessian matrix is complicated and easily irreversible. So, we proposed the MM algorithm [16] to get β^ for the adjustment of distorted covariates under case-cohort design.


### 4.1. Construct Surrogate

The key of the MM algorithm to construct surrogate function: combine the ideas of [16]; we expose the surrogate function Q^β|βm for l^pβ by using the convexity of the exponential function ex and the negative logarithm function −logx:(16)Q^β|βm=c^0+∑i=1nωiΔiβ⊤Z^iTi−∑l=1n∑r=1pYlTiωlλ^lrg^lrβr|βm∑l=1nYlTiωlexpZ^l⊤Tiβm,where βm is the mth approximate of the adjusted maximum likelihood estimation β^ defined in (9), and(17)c^0=∑i=1nωiΔi1−log∑l=1nYlTiωlexpZ^lTi⊤βm,g^lrβr|βm=expλ^lr−1Z^lrTiβr−βrm+Z^l⊤Tiβm,λ^lr=Z^lrTi∑r=1pZ^lrTi,with λ^lr−1=0 if Z^l⊤Ti=0, where βr, βrm, and Z^lr are the rth components of β, βm, and Z^l, respectively.With the surrogate function (16), we can transfer the optimization problem (9) into the problem as follows:(18)βm+1=argmaxβQ^β|βm.Following the construction process of the surrogate function (16), we can immediately conclude that(19)l^pβ≥Q^β|βm,where l^pβ=Q^β|βm if and only if β=βm.Then, withQ^βm+1|βm≥Q^βm|βm, we have(20)l^pβm+1=l^pβm+1−Q^βm+1|βm+Q^βm+1|βm≥l^pβm−Q^βm|βm+Q^βm|βm=l^pβm,where l^pβ is strictly ascending when βm+1≠βm.With the surrogate function (16), we obtain the derivatives with respect to β if β=βm as follows.The score vector:(21)∂Q^βm|βm∂β=∑i=1nΔiZ^Ti−∑l=1nYlTiωlZ^lTiexpZ^l⊤Tiβm∑l=1nYlTIωlexpZ^l⊤Tiβm.The negative Hessian matrix:(22)−∂2Q^βm|βm∂β∂β⊤=diag−∂2Q^βm|βm∂2β1,…,−∂2Q^βm|βm∂2βp,where(23)−∂2Q^βm|βm∂βr2=∑i=1nΔiωi∑l=1nYlTiωlZ^lr2/λ^lrexpZ^lTiβm∑l=1nYlTiωlexpZ^lTiβm,r=1,2,…,p.Here, the Hessian matrix is diagonal, the process of calculation changes simply and is reversible by using Newton–Raphson algorithm.


### 4.2. Cross-Validation and Bootstrap

#### 4.2.1. Cross-Validation

In the process of the adjustment for the distorted covariates, we use the cross-validation (CV) method to obtain the optimal bandwidthh to construct the kernel smoothing methods. The idea of the CV method is to calculate the coverage of error between some naive responses and the estimated responses obtained from the model fitted by the other predictors. The specific expression of the function in this article is given as follows:(24)CVh=∑j=1n1nh∑i=1nKuj−uih2Δu−2n∑j=1n2nh∑j≠iKuj−u−ih,where u−i is the set of removed ith element, Δu=uj−uj−1.Then, the optimal bandwidth is defined as follows:(25)h^=argminhCVh.


#### 4.2.2. Bootstrap

Bootstrap has been widely used since it was first proposed [24–26]. We adopt this method to calculate the standard error of β^; the basic idea of the nonparametric approach is to construct an empirical distribution function by repeatedly sampling from the observed data. The steps are summarized as follows:Step 1:X^1,X^2,…,X^n are the adjusted observations under the case-cohort design, where X^i=Ti,Δi,Ui,ξi+1−ξiΔiZ^i,1≤i≤n. We use the samples with Δi=1 as the failure cases and the samples with Δi=0 as the censoring objects to constitute the bootstrap case-cohort sample X^1∗,X^2∗,…,X^n∗.Step 2: a bootstrap replicationβ^∗ can be obtained based on X^1∗,X^2∗,…,X^n∗.Step 3: repeating step 1 and step 2B times, we can obtain B bootstrap estimate β^1∗,…,β^B∗. Therefore, the standard error of the rth component of β^ can be estimated by the following form:(26)Se^β^r=1B−1∑b=1Bβ^r∗b−1B∑b=1Bβ^r∗b2,r=1,2,…,p.


## 4.1. Construct Surrogate

The key of the MM algorithm to construct surrogate function: combine the ideas of [16]; we expose the surrogate function Q^β|βm for l^pβ by using the convexity of the exponential function ex and the negative logarithm function −logx:(16)Q^β|βm=c^0+∑i=1nωiΔiβ⊤Z^iTi−∑l=1n∑r=1pYlTiωlλ^lrg^lrβr|βm∑l=1nYlTiωlexpZ^l⊤Tiβm,where βm is the mth approximate of the adjusted maximum likelihood estimation β^ defined in (9), and(17)c^0=∑i=1nωiΔi1−log∑l=1nYlTiωlexpZ^lTi⊤βm,g^lrβr|βm=expλ^lr−1Z^lrTiβr−βrm+Z^l⊤Tiβm,λ^lr=Z^lrTi∑r=1pZ^lrTi,with λ^lr−1=0 if Z^l⊤Ti=0, where βr, βrm, and Z^lr are the rth components of β, βm, and Z^l, respectively.With the surrogate function (16), we can transfer the optimization problem (9) into the problem as follows:(18)βm+1=argmaxβQ^β|βm.Following the construction process of the surrogate function (16), we can immediately conclude that(19)l^pβ≥Q^β|βm,where l^pβ=Q^β|βm if and only if β=βm.Then, withQ^βm+1|βm≥Q^βm|βm, we have(20)l^pβm+1=l^pβm+1−Q^βm+1|βm+Q^βm+1|βm≥l^pβm−Q^βm|βm+Q^βm|βm=l^pβm,where l^pβ is strictly ascending when βm+1≠βm.With the surrogate function (16), we obtain the derivatives with respect to β if β=βm as follows.The score vector:(21)∂Q^βm|βm∂β=∑i=1nΔiZ^Ti−∑l=1nYlTiωlZ^lTiexpZ^l⊤Tiβm∑l=1nYlTIωlexpZ^l⊤Tiβm.The negative Hessian matrix:(22)−∂2Q^βm|βm∂β∂β⊤=diag−∂2Q^βm|βm∂2β1,…,−∂2Q^βm|βm∂2βp,where(23)−∂2Q^βm|βm∂βr2=∑i=1nΔiωi∑l=1nYlTiωlZ^lr2/λ^lrexpZ^lTiβm∑l=1nYlTiωlexpZ^lTiβm,r=1,2,…,p.Here, the Hessian matrix is diagonal, the process of calculation changes simply and is reversible by using Newton–Raphson algorithm.


## 4.2. Cross-Validation and Bootstrap

### 4.2.1. Cross-Validation

In the process of the adjustment for the distorted covariates, we use the cross-validation (CV) method to obtain the optimal bandwidthh to construct the kernel smoothing methods. The idea of the CV method is to calculate the coverage of error between some naive responses and the estimated responses obtained from the model fitted by the other predictors. The specific expression of the function in this article is given as follows:(24)CVh=∑j=1n1nh∑i=1nKuj−uih2Δu−2n∑j=1n2nh∑j≠iKuj−u−ih,where u−i is the set of removed ith element, Δu=uj−uj−1.Then, the optimal bandwidth is defined as follows:(25)h^=argminhCVh.


### 4.2.2. Bootstrap

Bootstrap has been widely used since it was first proposed [24–26]. We adopt this method to calculate the standard error of β^; the basic idea of the nonparametric approach is to construct an empirical distribution function by repeatedly sampling from the observed data. The steps are summarized as follows:Step 1:X^1,X^2,…,X^n are the adjusted observations under the case-cohort design, where X^i=Ti,Δi,Ui,ξi+1−ξiΔiZ^i,1≤i≤n. We use the samples with Δi=1 as the failure cases and the samples with Δi=0 as the censoring objects to constitute the bootstrap case-cohort sample X^1∗,X^2∗,…,X^n∗.Step 2: a bootstrap replicationβ^∗ can be obtained based on X^1∗,X^2∗,…,X^n∗.Step 3: repeating step 1 and step 2B times, we can obtain B bootstrap estimate β^1∗,…,β^B∗. Therefore, the standard error of the rth component of β^ can be estimated by the following form:(26)Se^β^r=1B−1∑b=1Bβ^r∗b−1B∑b=1Bβ^r∗b2,r=1,2,…,p.


## 4.2.1. Cross-Validation

In the process of the adjustment for the distorted covariates, we use the cross-validation (CV) method to obtain the optimal bandwidthh to construct the kernel smoothing methods. The idea of the CV method is to calculate the coverage of error between some naive responses and the estimated responses obtained from the model fitted by the other predictors. The specific expression of the function in this article is given as follows:(24)CVh=∑j=1n1nh∑i=1nKuj−uih2Δu−2n∑j=1n2nh∑j≠iKuj−u−ih,where u−i is the set of removed ith element, Δu=uj−uj−1.Then, the optimal bandwidth is defined as follows:(25)h^=argminhCVh.


## 4.2.2. Bootstrap

Bootstrap has been widely used since it was first proposed [24–26]. We adopt this method to calculate the standard error of β^; the basic idea of the nonparametric approach is to construct an empirical distribution function by repeatedly sampling from the observed data. The steps are summarized as follows:Step 1:X^1,X^2,…,X^n are the adjusted observations under the case-cohort design, where X^i=Ti,Δi,Ui,ξi+1−ξiΔiZ^i,1≤i≤n. We use the samples with Δi=1 as the failure cases and the samples with Δi=0 as the censoring objects to constitute the bootstrap case-cohort sample X^1∗,X^2∗,…,X^n∗.Step 2: a bootstrap replicationβ^∗ can be obtained based on X^1∗,X^2∗,…,X^n∗.Step 3: repeating step 1 and step 2B times, we can obtain B bootstrap estimate β^1∗,…,β^B∗. Therefore, the standard error of the rth component of β^ can be estimated by the following form:(26)Se^β^r=1B−1∑b=1Bβ^r∗b−1B∑b=1Bβ^r∗b2,r=1,2,…,p.


## 5. Simulation

As we mentioned before, if the distorting factors on covariates were ignored, the result of inference may be misled [11, 12]. In this article, we will construct several simulation studies to compare the three estimators of regression coefficient β^ for the covariates without distorting factors, with distorting factors, and adjusted with distorting factors in Cox’s model. Under the case-cohort design, where β^W denotes the estimator calculated based on the covariates without distorting factors, β^C denotes the estimator calculated based on the covariates with distorting factors, and β^A denotes the estimator calculate based on the adjusted covariates with distorting factors, given (Z1,Z2), we consider the hazard function of the failure time T˜i as follows:(27)λt|Z1,Z2=λ0texpβ1Z1+β2Z2.We set the true values of the parametersβ to (–0.25, 0.5), (–0.25, 0.693), and (0, 0.5). Z1 is generated from a normal distribution with mean Z and variance 1.44. Z2 is generated from a uniform distribution U7/3−19/2,7/3+19/2. The distorting factor U is generated from a uniform distribution U∼U4−7,4+7. The baseline hazard function λ0t is set to be 1 and 2t. Thus, the failure time T˜ satisfies an exponential distribution with failure rate expβ1Z1+β2Z2 and scale parameter expβ1Z1+β2Z2−1/2, respectively. The censoring time C is generated from a uniform distribution U0,c with c being chosen to desire the censoring rate ρ=80%, 85%, and 90%.Under the case-cohort design, 1000 full cohort samples are generated from function (27), and the corresponding distorting factor is generated. The two parts are merged to constitute the observed sample. Then, we randomly selected n˜=150, 200, and 300 subcohort from the full cohort.We compare the sample bias (Bias), the sample MSE (SMSE), the sample standard deviations of estimates (SD), the means of estimated standard errors (SE), and the coverage probabilities of 95% nominal confidence intervals (CP) of the three estimators toβ^C,β^W, and β^A based on 1000 independent simulated datasets and especially apply the nonparametric bootstrap approach in Section 4.2.2 with B=500. We set only the covariate Z2 to be distorted and adjusted. The distribution of the distorting factor U is set to be Φu=u+102/194.9160. The kernel function is set to be an Epanechnikov kernel [27] Kt=15/323−7t21−t2t≤1. The optimal bandwidth h is selected by the CV method in Section 4.2.1. The criteria are stopped with ε=0.0001. The simulation results are summarized in Table 1 with β=−0.25,0.5, Table 2 with β=−0.25,0.693, and Table 3 with β=0,0.5, respectively.Table 1
β=−0.25,0.5,Z1∼N2,1.44,Z2∼U7/3−19/2,7/3+19/2,u∼U4−7,4+7.
n˜ρ (%)Typeλ0t=1β1=−0.25β2=0.5BiasSMSESDSECPBiasSMSESDSECP15080β^W−0.01020.09650.09600.09920.9650.01460.10940.10850.11180.961β^C−0.00740.09610.09580.09780.959−0.11580.14990.09510.09750.742β^A−0.01000.09600.09550.09900.9640.00500.11050.11040.11260.95585β^W−0.01450.10510.10420.10530.9500.01150.11670.11620.11890.960β^C−0.00030.09680.09680.10300.956−0.11990.15750.10220.10240.724β^A−0.00140.09710.09710.10410.964−0.00460.12010.12010.11970.95090β^W−0.00980.11140.11100.11530.9610.02370.13060.12850.13280.957β^C−0.00940.10960.10930.11410.962−0.11290.15630.10810.11160.802β^A−0.00960.11130.11090.11520.9630.01060.12920.12880.13260.95420080β^W−0.00900.08150.08100.08780.9690.00840.09910.09880.09920.943β^C−0.00510.08070.08060.08680.971−0.12270.14970.08580.08520.669β^A−0.00910.08180.08140.08770.9680.00060.10120.10130.09940.94485β^W−0.00490.08990.08980.09260.9510.00050.10070.10080.10550.963β^C−0.00270.09280.09280.09140.940−0.12740.15600.09000.08970.666β^A−0.00480.09500.09490.09230.941−0.00260.10630.10630.10610.94990β^W−0.00940.10350.10310.10360.9540.01170.11740.11690.11960.962β^C−0.00890.10380.10350.10270.949−0.12290.15630.09650.09860.720β^A−0.00850.10340.10310.10340.953−0.00050.11510.11520.11880.96030080β^W−0.00290.07320.07320.07370.9510.00600.08480.08470.08500.954β^C−0.00080.07250.07260.07300.941−0.12410.14320.07140.07120.561β^A−0.00230.07270.07270.07360.950−0.00460.08440.08430.08480.95285β^W−0.00540.08040.08030.07950.9480.00690.09020.09000.09290.959β^C−0.00250.07680.07680.07920.953−0.13250.15320.07680.07610.568β^A−0.00360.07640.07630.07980.954−0.01140.09270.09210.09210.94790β^W−0.00340.09120.09120.09100.9580.00620.11050.11030.10780.950β^C−0.00230.09020.09020.09040.961−0.13080.15750.08770.08610.632β^A−0.00370.09080.09070.09090.957−0.00690.10930.10910.10630.947λ0t=2t15080β^W−0.00670.11070.11050.10920.9480.01030.12050.12010.12300.958β^C−0.00370.11100.11100.10780.943−0.11700.15710.10490.10820.773β^A−0.00600.10980.10970.10910.9480.00010.11980.11990.12380.95885β^W−0.00400.11110.11110.11430.9660.01610.13160.13060.13070.952β^C−0.00010.11040.11050.11270.958−0.11830.16060.10870.11390.798β^A−0.00360.11090.11090.11410.9640.00690.13030.13020.13140.95190β^W−0.00750.12420.12400.12540.9540.01020.14450.14420.14330.954β^C−0.00620.12280.12270.12380.954−0.12590.17370.11970.12160.784β^A−0.00660.12340.12330.12510.952−0.00370.14340.14340.14330.94720080β^W0.00160.09110.09110.09440.9510.00760.10220.10200.10740.957β^C0.00560.09160.09150.09320.952−0.12130.15210.09170.09340.710β^A0.00210.09180.09180.09430.954−0.00180.10240.10250.10780.95985β^W−0.00620.09710.09700.10110.9650.01160.11390.11330.11600.956β^C−0.00330.09730.09730.10000.959−0.11980.15420.09710.09970.747β^A−0.00600.09700.09690.10100.9640.00010.11410.11420.11640.94690β^W−0.00630.10600.10590.11240.9560.01170.12780.12730.12970.958β^C−0.00500.10610.10600.11130.957−0.12120.16080.10560.10850.746β^A−0.00690.10640.10620.11230.959−0.00070.12530.12540.12960.95630080β^W−0.00080.07630.07640.07900.9640.00360.08800.08790.09090.963β^C0.00240.07580.07580.07830.961−0.12380.14580.07700.07740.604β^A−0.00030.07660.07670.07890.963−0.00480.08860.08860.09090.95985β^W−0.00330.08710.08710.08540.9480.01200.09910.09850.09890.951β^C0.00020.08620.08620.08470.949−0.12270.14790.08270.08250.660β^A−0.00280.08700.08700.08530.9480.00180.09790.09790.09870.94290β^W−0.01270.09350.09270.09780.9690.01170.11490.11440.11410.945β^C−0.01100.09330.09270.09700.961−0.12720.15750.09290.09300.676β^A−0.01190.09290.09220.09770.967−0.00140.11250.11260.11310.946Note:β^W denotes the estimator calculated based on the covariates without distorting factors, β^C denotes the estimator calculated based on the covariates with distorting factors, and β^A denotes the estimator calculated based on the adjusted covariates with distorting factors. Algorithm: MM algorithm.Table 2
β=−0.25,0.693,Z1∼N2,1.44,Z2∼U7/3−19/2,7/3+19/2,u∼U4−7,4+7.
n˜ρ (%)Typeλ0t=1β1=−0.25β2=0.693BiasSMSESDβ^ASECPBiasSMSESDSECP15080β^W−0.00480.09980.09980.10380.9650.01470.11900.11810.12280.968β^C0.00050.09870.09880.10160.969−0.19210.21700.10090.10420.513β^A−0.00380.10050.10050.10380.963−0.00100.12370.12370.12360.94485β^W−0.01410.10820.10740.11080.9590.01100.13200.13160.13180.955β^C−0.00910.10420.10380.10870.964−0.19800.22580.10860.10880.517β^A−0.01370.10780.10700.11060.959−0.00580.13510.13510.13220.94790β^W−0.00600.11270.11260.12070.9740.01100.14620.14590.14460.953β^C−0.00380.11000.11000.11850.975−0.20560.23630.11650.11600.537β^A−0.00570.11390.11380.12040.967−0.01590.14530.14450.14430.94620080β^W−0.01050.09410.09360.09140.9490.01050.10890.10850.10920.953β^C−0.00170.09290.09290.08970.954−0.19690.21500.08630.09060.382β^A−0.00950.09370.09330.09130.948−0.00720.11110.11090.10970.94985β^W−0.00830.09730.09700.09790.9580.01150.10880.10830.11770.964β^C−0.00350.09550.09550.09600.954−0.19900.21780.08860.09490.446β^A−0.00790.09790.09770.09770.956−0.00870.11250.11220.11720.95390β^W−0.00460.10640.10640.10840.9560.01220.12640.12590.13260.971β^C−0.00220.10560.10570.10670.953−0.20180.22780.10570.10400.461β^A−0.00360.10600.10600.10820.957−0.01400.12810.12730.13120.95930080β^W−0.00250.07760.07760.07650.9530.00880.09320.09280.09280.955β^C0.00380.07660.07660.07520.955−0.20000.21390.07580.07450.260β^A−0.00110.07730.07740.07630.958−0.00900.09570.09530.09230.94585β^W−0.00370.08330.08330.08320.9470.00840.10410.10380.10240.947β^C0.00110.08220.08220.08190.945−0.20780.22230.07910.07910.257β^A−0.00270.08380.08380.08310.951−0.01450.10510.10420.10110.93490β^W−0.00350.09590.09590.09510.9540.00260.11530.11530.11810.949β^C0.00010.09610.09610.09400.949−0.21350.23020.08610.08850.326β^A−0.00210.09560.09570.09500.953−0.01950.11860.11700.11600.936λ0t=2t15080β^W−0.00900.11230.11200.11370.9520.01310.12750.12690.13380.954β^C0.00060.11100.11110.11150.951−0.18900.22090.11440.11510.579β^A−0.00730.11110.11090.11350.953−0.00320.13020.13030.13500.95185β^W−0.01090.12070.12030.12080.9490.01770.13960.13860.14340.962β^C−0.00850.11910.11880.11800.953−0.18640.22250.12160.12050.589β^A−0.01060.11980.11940.12050.946−0.00120.13870.13880.14390.95090β^W−0.01240.12870.12810.13230.9610.01580.16310.16240.15650.949β^C−0.01380.12770.12700.12960.955−0.19390.23050.12480.12920.621β^A−0.01510.12940.12860.13170.956−0.00720.15570.15560.15720.95520080β^W−0.01150.09940.09880.09860.9510.01620.12030.11930.11840.959β^C−0.00250.09410.09410.09630.956−0.19090.21540.09990.09910.492β^A−0.01090.09880.09830.09830.953−0.00110.12140.12150.11880.93985β^W−0.00070.10090.10100.10500.9600.00910.12600.12580.12690.954β^C0.00360.10080.10080.10310.957−0.19550.22180.10480.10420.489β^A−0.00090.10120.10120.10500.955−0.00810.12800.12780.12720.94990β^W−0.00510.11050.11050.11760.9620.01570.14840.14770.14290.943β^C−0.00110.10820.10820.11550.970−0.19960.23000.11440.11420.540β^A−0.00450.11030.11030.11730.961−0.01110.14340.14300.14160.94130080β^W−0.00670.07990.07970.08180.9520.00560.09340.09330.09890.962β^C0.00080.07890.07900.08050.948−0.19690.21140.07690.08070.322β^A−0.00580.08000.07980.08170.947−0.01120.09620.09560.09880.95285β^W−0.00830.08600.08570.08940.9570.00810.10630.10610.10820.961β^C−0.00040.08540.08550.08790.965−0.20180.21940.08610.08640.352β^A−0.00720.08670.08640.08930.959−0.01150.10790.10730.10780.95090β^W−0.00820.09660.09630.10070.9610.01250.12230.12170.12590.963β^C−0.00340.09480.09480.09910.951−0.20120.22190.09350.09640.428β^A−0.00600.09640.09630.10060.955−0.01120.12060.12010.12390.954Note:β^W denotes the estimator calculated based on the covariates without distorting factors, β^C denotes the estimator calculated based on the covariates with distorting factors, and β^A denotes the estimator calculate based on the adjusted covariates with distorting factors. Algorithm: MM algorithm.Table 3
β=0,0.5,Z1∼N2,1.44,Z2∼U7/3−19/2,7/3+19/2,u∼U4−7,4+7.
n˜ρ (%)Typeλ0t=1β1=0β2=0.5BiasSMSESDSECPBiasSMSESDSECP15080β^W0.00160.08990.08990.09300.9610.01060.10320.10270.10840.968β^C0.00160.09000.09010.09200.963−0.11840.14900.09060.09540.759β^A0.00170.09030.09030.09290.9570.00140.10580.10580.10910.96285β^W−0.00130.10100.10100.09990.9590.00790.11070.11040.11510.964β^C−0.00070.09910.09910.09890.961−0.11910.15320.09640.09980.728β^A−0.00110.10120.10120.09970.958−0.00520.11140.11130.11540.95290β^W−0.00030.10780.10790.10830.9540.01230.12420.12360.12820.965β^C−0.00040.10830.10840.10740.946−0.12180.15760.10000.10760.772β^A−0.00030.10780.10790.10820.956−0.00020.12250.12250.12820.96520080β^W0.00160.08070.08070.08230.9630.00730.09730.09710.09630.951β^C0.00290.07970.07960.08130.962−0.12390.14880.08240.08270.633β^A0.00170.08030.08040.08220.966−0.00180.09840.09850.09650.94385β^W0.00070.08620.08630.08810.9620.00710.10240.10220.10360.953β^C0.00100.08450.08460.08730.960−0.12300.15150.08840.08780.672β^A0.00050.08650.08650.08810.965−0.00480.10320.10310.10370.94890β^W0.00370.09400.09400.09860.9540.00910.11480.11450.11680.959β^C0.00440.09450.09440.09770.958−0.12700.15800.09410.09600.688β^A0.00340.09450.09450.09850.956−0.00360.11580.11580.11630.95330080β^W0.00160.07060.07060.07030.9510.00120.08030.08040.08270.964β^C0.00100.07050.07060.06960.942−0.12940.14610.06790.06920.526β^A0.00190.07060.07060.07020.945−0.00810.08160.08130.08260.95585β^W−0.00030.07670.07670.07610.947−0.00060.08680.08680.09080.960β^C0.00010.07650.07650.07560.947−0.13210.15040.07190.07450.545β^A−0.00020.07670.07670.07610.946−0.01110.08840.08780.09040.95090β^W0.00070.08550.08550.08650.9600.00950.10460.10420.10540.954β^C0.00130.08450.08450.08610.969−0.12740.15210.08310.08370.638β^A0.00100.08560.08570.08650.959−0.00270.10380.10380.10400.950λ0t=2t15080β^W0.00050.10120.10130.10350.9570.00910.11440.11400.11950.953β^C0.00120.09980.09980.10240.959−0.11910.15700.10240.10570.755β^A0.00100.10070.10070.10340.9560.00200.11580.11580.12010.95485β^W−0.00200.10590.10590.10940.9650.01320.12100.12030.12720.971β^C−0.00310.10480.10480.10800.965−0.11850.15900.10620.11090.798β^A−0.00200.10570.10580.10910.9650.00160.11980.11980.12790.96690β^W−0.00380.11570.11570.11890.9620.00690.13760.13750.13890.957β^C−0.00400.11340.11340.11750.964−0.12690.16980.11280.11820.781β^A−0.00340.11600.11600.11870.959−0.00360.13870.13880.13930.95520080β^W0.00420.09060.09050.09050.9440.01000.10570.10530.10610.962β^C0.00440.09030.09020.08920.944−0.11960.15090.09210.09240.696β^A0.00470.09040.09030.09030.942−0.00020.10500.10500.10650.95685β^W−0.00180.09500.09500.09540.9620.00930.10890.10880.11240.962β^C−0.00160.09490.09500.09420.952−0.12580.15660.09330.09600.700β^A−0.00170.09440.09450.09520.957−0.00430.10920.10920.11250.95890β^W−0.00010.10070.10080.10650.9720.01440.12380.12310.12630.965β^C0.00090.09910.09910.10540.968−0.12110.15870.10260.10550.763β^A0.00000.10040.10050.10640.9740.00260.12390.12390.12640.95930080β^W−0.00100.07320.07320.07540.9590.00890.08770.08730.08900.965β^C−0.00110.07260.07260.07460.962−0.12160.14390.07690.07570.605β^A−0.00110.07330.07340.07540.960−0.00120.08790.08790.08910.96585β^W0.00200.08100.08100.08130.9540.00530.09630.09620.09690.959β^C0.00110.08000.08000.08080.952−0.12750.15080.08070.08120.628β^A0.00150.08130.08130.08120.952−0.00530.09560.09550.09670.95290β^W0.00440.09260.09260.09250.9530.01320.10920.10840.11180.961β^C0.00450.09230.09230.09190.954−0.12570.15300.08720.09060.687β^A0.00410.09220.09220.09240.9490.00130.10840.10850.11080.960Note:β^W denotes the estimator calculated based on the covariates without distorting factors, β^C denotes the estimator calculated based on the covariates with distorting factors, and β^A denotes the estimator calculate based on the adjusted covariates with distorting factors. Algorithm: MM algorithm.According to the simulation results, the estimatorsβ^W and β^A for β2 are all unbiased, the three values of SMSE, SD, and SE are all closed, and the CP values are reasonable. However, the estimators β^C for β2 are all biased, the difference between the two values of SMSE and SD is large, and the CP value is low. These facts confirm that ignoring the distorting factors may lead to the wrong results and conclusion and also explain that the proposed method about the adjustment of covariates is effective. For example, the estimator β^C and β^A for β=−0.25,0.5, ρ=80%, n˜=150, for β2, the bias is –0.1178, SMSE are 0.1521 and 0.1095, and CP are 0.735 and 0.947, respectively. Furthermore, we can conclude that the estimators for β1 are all reasonable.


## 6. Real Data Analysis

The dataset is about the medical records of heart failure patients, from “Machine Learning Repository” (http://archive.ics.uci.edu/ml/datasets/Heart+failure+clinical+records). The original dataset version was collected by Tanvir Ahmad et al. (Government College University, Faisalabad, Pakistan) and made available by them on FigShare under the Attribution 4.0 International copyright in July 2017. The current version of the dataset was elaborated by Davide Chicco et al. and donated to the University of California Irvine Machine Learning Repository under the same Attribution 4.0 International copyright in January 2020.This dataset contains the medical records of 299 patients who had heart failure, collected during their follow-up period, where each patient profile has 13 clinical features, including age, anaemia, creatinine phosphokinase, diabetes, ejection fraction, high blood pressure, platelets, serum creatinine, serum sodium, sex, smoking, time, and death event. The goal is to assess the association between clinical features and heart failure and its happened time. This real data is survival data which meets the requirements of the article. Death event and time are regarded as heart failure and its happened time, respectively. For the remaining 11 clinical features, through pairwise correlation analysis, we conclude that ejection fraction and serum sodium and anaemia and Creatinine phosphokinase are significantly related. Referring to relevant medical knowledge, we know that ejection fraction may distort serum sodium. And in the same way, anaemia may distort creatinine phosphokinase. The main purpose of our research is to evaluate the association between these clinical features under smoothed serum sodium and creatinine phosphokinase by ejection fraction and anaemia.There are 96 deaths in this data, and the censoring rate is nearby67.9%. Due to the different dimensions of the covariates, we firstly standardize each covariate. Then, we set the standardized age as Z1, the standardized creatinine phosphokinase as Z^2, the standardized diabetes as Z3, the standardized high blood pressure as Z4, the standardized platelets as Z5, the standardized serum creatinine as Z6, the standardized serum sodium as Z^7, the standardized sex as Z8, the standardized smoking as Z9, the standardized ejection fraction as a distorting factor U1, and the standardized anaemia as a distorting factor U2, where Z^2 and Z^7 indicate that the covariate is distorted.The case-cohort sample is composed of the subcohort randomly selected 100 from the whole data and the dead objects in the whole data expect the subcohort. We perform regression analysis on the sample using the proportional hazard function:(28)xt|Z=x0texp∑i=19βiZi.Because the size of this data is small, to improve the reliability of the sampling method, the sampling process is repeated 100 times. The average value of the 100 distorted estimates is2.7224,14.0226,0.1603,0.4295,−0.2667,3.0520,11.0269,−0.0742,0.1933⊤.Then, for the same 100 case-cohort sample, we obtain the adjusted covariatesZ^2 and Z^7 by using the proposed kernel smoothing method on Z˜2 and Z˜7 and replace Z˜2 and Z˜7 with Z^2 and Z^7, respectively. On this basis, the average of adjusted estimates is 2.5077,1.0344,0.1136,0.5474,−0.3552,2.8577,−2.4303,0.2115⊤.The kernel and bandwidth selector chosen are the same as those inSection 5. After the covariates were adjusted, the estimator of creatinine phosphokinase is adjusted from 1.0344 to 14.0226, and the estimator of serum sodium is adjusted from −2.4303 to 11.0269; the estimators of other covariates changed little after adjustment. The result is reasonable and feasible in medicine and also shows that the proposed method is effective for real data.In practice, the covariates are measured after the case-cohort sampling to reduce costs and improve efficiency. In this real data, the size of the data is small, and the each covariate has been measured. To reduce the error, we repeat 100 samplings and set the average value as the estimators. For example, based on the 10th case-cohort sample from the original data, the distorted estimator is3.5393,−0.5883,0.1934,0.3712,0.4114,6.5332,−3.9763,0.1236,0.3585⊤; then, after the adjustment of covariates in the same sample, the adjusted estimator is 4.5125,79.7599,0.1690,0.2902,1.0486,5.9035,12.7945,0.1399,0.1720⊤. Therefore, the estimators based on same case-cohort sample are quite different from the average of estimators, but the trends are the same. So, for real data with a large sample size, we only need to draw a case-cohort sample and measure their covariates.


## 7. Conclusion

We study the estimation in Cox’s model with the adjusted covariates by using the nonparametric method based on the kernel function under case-cohort design. Consistency and asymptotic normality of the proposed estimators are derived. We use MM algorithm for calculating the regression coefficients, where the surrogate function with a diagonal Hessian matrix is established to overcome the operation difficulty in the Newton–Raphson algorithm. Simulation and real data studies suggest that the case-cohort design can be used to reduce the cost for cohort studies and the development of the adjustment for distorted covariates has nice performance.Here, the adjustment for survival data is discussed in Cox’s model. The method can be extended to other models, such as an accelerated failure time model [28], an additive-multiplicative model [29], and an accelerated hazards model [30]. The cost-effective case-cohort design is mainly applied to the data with the high censoring rate. However, in the sampling process, some censoring subjects have not been measured. To improve the efficiency, our future works include developments of inference for estimation with adjusted covariates under a generalized case-cohort design [31] and an outcome-dependent sampling design.



---

*Source: 8884665-2020-12-24.xml*