# New Inference Procedures for Semiparametric Varying-Coefficient Partially Linear Cox Models

Authors: Yunbei Ma; Xuan Luo
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

In biomedical research, one major objective is to identify risk factors and study their risk impacts, as this identification can help clinicians to both properly make a decision and increase efficiency of treatments and resource allocation. A two-step penalized-based procedure is proposed to select linear regression coefficients for linear components and to identify significant nonparametric varying-coefficient functions for semiparametric varying-coefficient partially linear Cox models. It is shown that the penalized-based resulting estimators of the linear regression coefficients are asymptotically normal and have oracle properties, and the resulting estimators of the varying-coefficient functions have optimal convergence rates. A simulation study and an empirical example are presented for illustration.


---

## Body

## 1. Introduction

To balance model flexibility and specificity, a variety of semiparametric models have been proposed in survival analysis. For example, Huang [1] studied partially linear Cox models and Cai et al. [2] and Fan et al. [3] studied the Cox proportional hazard model with varying coefficients. Tian et al. [4] proposed an estimation procedure for the Cox model with time-varying coefficients. Perhaps the most appropriate model is the semiparametric varying-coefficient model; its merit includes easy interpretation, flexible structure, potential interaction between covariates, and dimension reduction of nonparametric components models. Recently, a lot of effort has been made in this direction. Cai et al. [2] considered the Cox proportional hazard model with a semiparametric varying-coefficient structure. Yin et al. [5] proposed the semiparametric additive hazard model with varying coefficients.In biomedical research and clinical trials, one major objective is to identify risk factors and study their risk impacts, as this identification can help clinicians to properly make a decision and increase efficiency of treatment and resource allocation. This is essentially a variable selection procedure in spirit. When the outcomes are censored, selection of significant risk factors becomes more complicated and challenging than with the case where all outcomes are complete. However, important progress has been made in recent literature. For instance, Fan and Li [6] extended the nonconcave penalized likelihood approach proposed by Fan and Li [7] to the Cox proportional hazards model and the Cox proportional hazards frailty model. Johnson et al. [8, 9] proposed procedures for selecting variables in semiparametric linear regression models for censored data. However, these papers did not consider variable selection with functional coefficients; while Wang et al. [10] extended the application of penalized likelihood to the varying-coefficient setting and studied the asymptotic distributions of the estimators, they only focused on the SCAD penalty for completely observed data. Du et al. [11] proposed a penalized variable selection procedure for Cox models with a semiparametric relative risk.In this paper, we study variable selection for both regression parameters and varying coefficients in the semiparametric varying-coefficient additive hazards model:(1)
λ
(
t
∣
Z
,
V
,
W
)
=
λ
0
(
t
)
+
β
T
(
W
(
t
)
)
Z
(
t
)
+
α
T
V
(
t
)
+
g
(
W
(
t
)
)
,
where V
(
t
) is a vector of covariates with a linear effect on the logarithm of the hazard function λ
(
t
), W
(
t
) is the main exposure variable of interest whose effect on the logarithm of the hazard might be nonlinear, and Z
(
·
)
=
(
Z
1
(
·
)
,
…
,
Z
p
(
·
)
)
T is a vector of covariates that may interact with the exposure covariate W
(
·
). λ
0
(
·
) is the baseline hazard ratio function and g
(
·
) is an unspecified smooth function. α is a d-vector of constant parameters. Let g
(
0
)
=
0 to assure that the model is identifiable. This model covers generally used basic semiparametric model settings. For instance, if β
(
·
)
≡
0 and g
(
·
)
≡
0, model (1) reduces to Lin and Ying’s additive hazard model [12] and Aalen’s additive hazard model when the baseline hazard function λ
0
(
t
) also equals to zero [13, 14]. Furthermore, when the exposure covariate W
(
·
) is time t, then model (1) reduces to the partly parametric additive hazard model if λ
0
(
t
)
≡
g
(
t
)
≡
0 [15, 16] and reduces to the time-varying coefficient additive hazard model [17].To select significant elements ofα and coefficient functions β
(
·
), we require an estimation procedure for the regression parameters α and β
(
·
). However, this poses a challenge because it is impossible to simultaneously get the root-n consistent penalized estimators for the nonzero components of α by penalizing β
(
·
) and α. To achieve this goal, we propose the following two-step estimation procedure for selecting the relevant variables. In Step 1, we estimate α by using the profile penalized least square function after locally approximating the nonparametric functions β
(
·
) and g
(
·
). In Step 2, given the penalized estimator of α in Step 1, we develop a penalized least square function for β
(
·
) by using the basis expansion. We will demonstrate that the proposed estimation procedures have oracle properties, and all penalized estimators of nonzero components can achieve their optimal convergence rate.This paper is organized as follows. Section2 introduces the SCAD-based variable selection procedure for the parametric components and establishes the asymptotic normality for the resulting penalized estimators. Section 3 discusses a variable selection procedure for the coefficient functions by approximating these functions by using the B-spline approach. The simulation studies and an application of the proposed methods in a real data example are included in Sections 4 and 5, respectively. The technical lemmas and proofs of the main results are given in the appendix.


## 2. Penalized Least Squares Based Variable Selection for Parametric Components

LetT
i denote the potential failure time, let C
i denote the potential censoring time, and let X
i
=
min
⁡
(
T
i
,
C
i
) denote the observed time for the ith individual, i
=
1
,
…
,
n. Let Δ
i be an indicator which equals 1 if X
i is a failure time and 0 otherwise. Let F
t
,
i represent the failure, censoring, and covariate information up to time t for the ith individual. The observed data structure is {
X
i
,
Δ
i
,
Z
i
(
t
)
,
V
i
(
t
)
,
W
i
(
t
)
,
i
=
1
,
…
,
n
}. Assume that T and C are conditionally independent given covariates and that the observation period is [
0
,
τ
], where τ is a constant denoting the time for the end of the study.LetN
i
(
t
)
=
I
(
X
i
≤
t
,
Δ
i
=
1
) denote the counting process corresponding to T
i and let Y
i
(
t
)
=
I
(
X
i
≥
t
). Let the filtration {
F
t
:
t
∈
[
0
,
τ
]
} be the history up to time t; that is, F
t
=
σ
{
X
i
≤
u
,
Z
i
(
u
)
,
V
i
(
u
)
,
W
i
(
u
)
,
Y
i
(
u
)
,
Δ
i, 0
≤
u
≤
t
,
i
=
1
,
…
,
n
}. Write M
i
(
t
)
=
N
i
(
t
)
-
∫
0
T
‍
Y
i
(
u
)
λ
i
(
u
)
d
u. Then M
i
(
t
) is a martingale with respect to F
t. For ease of presentation, we drop the dependence of covariates on time. The methods and proofs in this paper are applicable to external time-dependent covariates [18].Here we use the techniques of local linear fitting [19]. Suppose that β
(
·
) and g
(
·
) are smooth enough to allow Taylor’s expansion as follows: for each w
0
∈
W which is the support of W and for w in a neighborhood of w
0,
(2)
β
(
w
)
≈
β
(
w
0
)
+
β
′
(
w
0
)
(
w
-
w
0
)
,
g
(
w
)
≈
g
(
w
0
)
+
g
′
(
w
0
)
(
w
0
-
w
)
.
Then we can approximate the hazard ratio function given in (1) by
(3)
λ
(
t
∣
Z
i
,
V
i
,
W
i
)
≈
λ
0
*
(
t
,
w
0
)
+
η
T
(
w
0
)
Z
i
*
(
w
0
)
+
α
T
V
i
,
where η
(
w
0
)
=
(
β
T
(
w
0
)
,
(
β
′
(
w
0
)
)
T
,
g
′
(
w
0
)
)
T, Z
i
*
(
w
0
)
=
(
Z
i
T
,
Z
i
T
(
W
i
-
w
0
)
,
(
W
i
-
w
0
)
)
T, and λ
0
*
(
t
,
w
0
)
=
λ
0
(
t
)
+
g
(
w
0
).LetH be a 2
p
+
1 diagonal matrix with the first p diagonal elements 1 and the rest h. If α is given, then by using the counting-process notation, similarly to [5], we can obtain an estimator of η
(
w
0
) at each w
0
∈
W by solving the estimating equations whose objective function is given as follows:
(4)
U
n
(
η
,
w
0
)
=
∑
i
=
1
n
‍
∫
0
τ
‍
K
h
(
W
i
-
w
0
)
H
-
1
{
Z
i
*
(
w
0
)
-
Z
¯
(
t
,
w
0
)
}
n
n
n
n
n
n
n
n
n
×
[
{
η
T
(
w
0
)
Z
i
*
(
w
0
)
+
α
T
V
i
}
d
u
]
d
N
i
(
t
)
-
Y
i
(
t
)
n
n
n
n
n
n
n
n
n
n
n
n
×
{
η
T
(
w
0
)
Z
i
*
(
w
0
)
+
α
T
V
i
}
d
u
]
,
where
(5)
Z
¯
(
t
,
w
0
)
=
∑
i
=
1
n
K
h
(
W
i
-
w
0
)
Y
i
(
t
)
Z
i
*
(
w
0
)
∑
i
=
1
n
K
h
(
W
i
-
w
0
)
Y
i
(
t
)
,
K
h
(
·
)
=
K
(
·
/
h
)
/
h with K
(
·
) being a kernel function and h being a bandwidth, and τ
=
∞. To avoid the technicality of tail problems, only data up to a finite time point τ are frequently used.Denote the solution ofU
n
(
η
,
w
0
)
=
0 by η
^
(
w
0
,
α
), which can be expressed as follows:
(6)
η
^
(
w
0
,
α
)
=
[
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
K
h
(
W
i
-
w
0
)
H
-
1
n
n
n
n
n
n
n
n
n
n
n
×
{
Z
i
*
(
w
0
)
-
Z
¯
(
t
,
w
0
)
}
⊗
2
Y
i
(
t
)
d
t
∑
i
=
1
n
‍
]
-
1
×
[
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
K
h
(
W
i
-
w
0
)
H
-
1
n
n
n
n
×
{
Z
i
*
(
w
0
)
-
Z
¯
(
t
,
w
0
)
}
{
d
N
i
(
t
)
-
Y
i
(
t
)
V
i
T
α
d
t
}
∑
i
=
1
n
‍
]
≝
M
n
1
-
1
(
w
0
)
{
S
n
1
(
w
0
)
-
S
n
2
(
w
0
)
α
}
.
Here and below, a
⊗
k
=
1
,
a
,
a
a
T for any vector a and for k
=
0,1
,
2.As a consequence, we can estimateβ
(
·
) and g
(
·
) presuming that α is known:
(7)
β
^
(
W
i
,
α
)
=
(
I
p
,
0
p
×
p
,
0
p
)
M
n
1
-
1
(
W
i
)
{
S
n
1
(
W
i
)
-
S
n
2
(
W
i
)
α
}
≡
T
n
1
(
W
i
)
-
T
n
2
(
W
i
)
α
,
g
^
(
W
i
,
α
)
=
∫
0
W
i
‍
(
0
p
×
p
,
0
p
×
p
,
1
p
)
M
n
1
-
1
(
w
0
)
nnnnnn
×
{
S
n
1
(
w
0
)
-
S
n
2
(
w
0
)
α
}
d
w
0
≡
T
n
3
(
W
i
)
-
T
n
4
(
W
i
)
α
.


### 2.1. Variable Selection

Notice that ifα is given, then from (7) we can estimate Λ
0
(
t
,
W
i
)
=
∫
0
T
‍
λ
0
(
u
,
W
i
)
d
u by
(8)
Λ
^
0
(
t
,
W
i
,
α
)
=
∑
i
=
1
n
‍
∫
0
t
‍
[
d
N
i
(
u
)
-
Y
i
(
u
)
{
Z
i
T
T
n
1
(
W
i
)
+
T
n
3
(
W
i
)
}
d
u
]
∑
i
=
1
n
Y
i
(
u
)
-
∑
i
=
1
n
‍
∫
0
t
‍
Y
i
(
u
)
[
V
i
T
-
Z
i
T
T
n
2
(
W
i
)
-
T
n
4
(
W
i
)
]
α
d
u
∑
i
=
1
n
Y
i
(
u
)
≝
∑
i
=
1
n
‍
∫
0
t
‍
∑
i
=
1
n
[
d
N
i
(
u
)
-
Y
i
(
u
)
{
G
n
i
+
V
i
*
(
W
i
)
α
}
d
u
]
∑
i
=
1
n
Y
i
(
u
)
.
Given the observed data and based on (7) and (8), β
(
·
), g
(
·
), and the cumulative function Λ
0
(
·
) can be approximately expressed as functions of α. These expressions inspire us to estimate α by minimizing the following objective function subject to α:
(9)
L
n
(
α
)
=
1
2
\
d
i
s
p
l
a
y
s
t
y
l
e
∑
i
=
1
n
‍
(
∫
0
τ
‍
[
d
N
i
(
t
)
-
Y
i
(
t
)
d
Λ
^
0
(
t
,
W
i
,
α
)
n
n
n
n
n
n
n
-
Y
i
(
t
)
{
Z
i
T
β
^
(
W
i
,
α
)
+
g
^
(
W
i
,
α
)
n
n
n
n
n
n
n
n
n
n
n
n
n
+
V
i
T
α
}
d
t
]
∫
0
τ
‍
{
d
N
i
(
t
)
-
Y
i
(
t
)
d
Λ
^
0
(
t
,
W
i
,
α
)
)
2
.
Accordingly, our penalized objective function is defined as follows:
(10)
L
n
*
(
α
)
=
L
n
(
α
)
+
n
∑
j
=
1
d
‍
p
λ
(
|
α
j
|
)
.
Thus, minimizing L
n
*
(
α
) with respect to α results in a penalized least squares estimator of α.Leta
n
=
max
⁡
{
p
λ
n
′
(
|
α
0
j
|
)
:
α
0
j
≠
0
} and b
n
=
max
⁡
{
p
λ
n
′′
(
|
α
0
j
|
)
:
α
0
j
≠
0
}. Now we establish the oracle properties for the resulting estimators. Let α
0
=
(
α
10
T
,
α
20
T
)
T be the true value of α, the first subvector α
10
T of length s, containing all nonzero elements of α and α
20
T
=
0
d
-
s.We introduce an additional notation as follows. LetZ
~
=
(
Z
T
,
0
p
T
,
0
)
T. For k
=
0,1, let μ
k
=
∫
‍
x
k
K
(
x
)
d
x,   ν
k
=
∫
‍
x
k
K
2
(
x
)
d
x,  ρ
k
(
u
,
z
,
v
,
w
0
)
=
P
(
X
≥
u
∣
Z
=
z
,
V
=
v
,
W
=
w
0
)
λ
k
(
t
∣
z
,
v
,
w
0
). Denote
(11)
φ
k
,
l
(
t
,
w
0
)
=
f
(
w
0
)
E
{
ρ
0
(
t
,
Z
,
v
,
w
0
)
Z
~
(
V
⊗
l
)
T
⊗
k
∣
W
=
w
0
}
,
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
nn
n
n
n
k
=
0,1
,
l
=
0,1
,
φ
2,0
(
t
,
w
0
)
=
f
(
w
0
)
E
{
(
Z
⊗
2
0
p
×
p
0
p
0
p
×
p
μ
2
Z
⊗
2
Z
μ
2
0
p
T
Z
T
μ
2
μ
2
)
ρ
0
(
t
,
Z
,
v
,
w
0
)
n
n
n
n
n
n
n
n
n
n
n
n
n
×
(
Z
⊗
2
0
p
×
p
0
p
0
p
×
p
μ
2
Z
⊗
2
Z
μ
2
0
p
T
Z
T
μ
2
μ
2
)
∣
W
=
w
0
}
,
where f
(
·
) is the density of W. Denote t
1
(
w
)
=
(
I
p
,
0
p
×
p
,
0
p
)
m
1
-
1
(
w
)
m
2
(
w
)  and t
2
(
w
)
= (
0
p
×
p
,
0
p
×
p
,
1
p
) ∫
0
w
‍
m
1
-
1
(
w
0
)
m
2
(
w
0
)
d
w
0, where
(12)
m
1
(
w
0
)
=
∫
0
τ
‍
{
φ
2,0
(
t
,
w
0
)
-
φ
1,0
⊗
2
(
t
,
w
0
)
φ
0,0
(
t
,
w
0
)
}
d
t
,
m
2
(
w
0
)
=
∫
0
τ
‍
{
φ
1,1
(
t
,
w
0
)
-
φ
1,0
(
t
,
w
0
)
φ
0,0
(
t
,
w
0
)
φ
0,1
(
t
,
w
0
)
}
d
t
.
For k
=
0,1
,
2, l
=
0,1, let κ
k
,
l
(
t
)
=
E
[
ρ
l
(
t
,
Z
,
V
,
W
)
{
V
-
Z
t
1
T
(
W
)
-
t
2
T
(
W
)
}
⊗
k
], and V
~
i
(
W
i
)
=
V
i
T
-
Z
i
T
t
1
(
W
i
)
-
t
2
(
W
i
),  i
=
1
,
…
,
n
.
∥
a
∥
=
(
a
T
a
)
1
/
2 for any vector a.Theorem 1.
If Assumptions (A.i)–(A.vi) are satisfied andb
n
→
0 as n
→
∞, then there exist local minimizers α
^ of L
n
*
(
α
) such that ∥
α
^
-
α
0
∥
=
O
p
(
n
-
1
/
2
+
a
n
).Remark 2.
Theorem1 indicates that by choosing a proper λ
n, which leads to a
n
=
O
(
n
-
1
/
2
), there exists a root-n consistent penalized least squares estimator of α.Theorem 3.
Assume that(13)
∞
>
lim
⁡
inf
⁡
n
→
∞
lim
⁡
inf
⁡
θ
→
0
+
p
λ
n
′

(
θ
)
λ
n
>
0
.
Under Assumptions (A.i)–(A.vi), if λ
n
→
0 and n
λ
n
→
∞ as n
→
∞, then the local root-n consistent local minimizer α
^
=
(
α
^
1
T
,
α
^
2
T
)
T in Theorem 1 with probability tending to 1 satisfies(1)
α
^
2
=
0; and
(2)
asymptotic normality(14)
n
(
I
1
(
τ
)
+
Γ
)
×
{
α
^
1
-
α
10
-
(
I
1
(
τ
)
+
Γ
)
-
1
b
}
⟶
D
N
(
0
,
Σ
1
(
τ
)
)
,where I
1
(
τ
) and Σ
1
(
τ
) are the first s
×
s submatrix of I
(
τ
)n and Σ
(
τ
), respectively, and
(15)
b
=
(
p
λ
n
′
(
|
α
10
|
)
sgn
⁡
(
α
10
)
,
…
,
p
λ
n
′
(
|
α
s
0
|
)
sgn
⁡
(
α
s
0
)
)
T
,
Γ
=
diag
⁡
(
p
λ
n
′′
(
|
α
10
|
)
,
…
,
p
λ
n
′′
(
|
α
s
0
|
)
)
.


## 2.1. Variable Selection

Notice that ifα is given, then from (7) we can estimate Λ
0
(
t
,
W
i
)
=
∫
0
T
‍
λ
0
(
u
,
W
i
)
d
u by
(8)
Λ
^
0
(
t
,
W
i
,
α
)
=
∑
i
=
1
n
‍
∫
0
t
‍
[
d
N
i
(
u
)
-
Y
i
(
u
)
{
Z
i
T
T
n
1
(
W
i
)
+
T
n
3
(
W
i
)
}
d
u
]
∑
i
=
1
n
Y
i
(
u
)
-
∑
i
=
1
n
‍
∫
0
t
‍
Y
i
(
u
)
[
V
i
T
-
Z
i
T
T
n
2
(
W
i
)
-
T
n
4
(
W
i
)
]
α
d
u
∑
i
=
1
n
Y
i
(
u
)
≝
∑
i
=
1
n
‍
∫
0
t
‍
∑
i
=
1
n
[
d
N
i
(
u
)
-
Y
i
(
u
)
{
G
n
i
+
V
i
*
(
W
i
)
α
}
d
u
]
∑
i
=
1
n
Y
i
(
u
)
.
Given the observed data and based on (7) and (8), β
(
·
), g
(
·
), and the cumulative function Λ
0
(
·
) can be approximately expressed as functions of α. These expressions inspire us to estimate α by minimizing the following objective function subject to α:
(9)
L
n
(
α
)
=
1
2
\
d
i
s
p
l
a
y
s
t
y
l
e
∑
i
=
1
n
‍
(
∫
0
τ
‍
[
d
N
i
(
t
)
-
Y
i
(
t
)
d
Λ
^
0
(
t
,
W
i
,
α
)
n
n
n
n
n
n
n
-
Y
i
(
t
)
{
Z
i
T
β
^
(
W
i
,
α
)
+
g
^
(
W
i
,
α
)
n
n
n
n
n
n
n
n
n
n
n
n
n
+
V
i
T
α
}
d
t
]
∫
0
τ
‍
{
d
N
i
(
t
)
-
Y
i
(
t
)
d
Λ
^
0
(
t
,
W
i
,
α
)
)
2
.
Accordingly, our penalized objective function is defined as follows:
(10)
L
n
*
(
α
)
=
L
n
(
α
)
+
n
∑
j
=
1
d
‍
p
λ
(
|
α
j
|
)
.
Thus, minimizing L
n
*
(
α
) with respect to α results in a penalized least squares estimator of α.Leta
n
=
max
⁡
{
p
λ
n
′
(
|
α
0
j
|
)
:
α
0
j
≠
0
} and b
n
=
max
⁡
{
p
λ
n
′′
(
|
α
0
j
|
)
:
α
0
j
≠
0
}. Now we establish the oracle properties for the resulting estimators. Let α
0
=
(
α
10
T
,
α
20
T
)
T be the true value of α, the first subvector α
10
T of length s, containing all nonzero elements of α and α
20
T
=
0
d
-
s.We introduce an additional notation as follows. LetZ
~
=
(
Z
T
,
0
p
T
,
0
)
T. For k
=
0,1, let μ
k
=
∫
‍
x
k
K
(
x
)
d
x,   ν
k
=
∫
‍
x
k
K
2
(
x
)
d
x,  ρ
k
(
u
,
z
,
v
,
w
0
)
=
P
(
X
≥
u
∣
Z
=
z
,
V
=
v
,
W
=
w
0
)
λ
k
(
t
∣
z
,
v
,
w
0
). Denote
(11)
φ
k
,
l
(
t
,
w
0
)
=
f
(
w
0
)
E
{
ρ
0
(
t
,
Z
,
v
,
w
0
)
Z
~
(
V
⊗
l
)
T
⊗
k
∣
W
=
w
0
}
,
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
nn
n
n
n
k
=
0,1
,
l
=
0,1
,
φ
2,0
(
t
,
w
0
)
=
f
(
w
0
)
E
{
(
Z
⊗
2
0
p
×
p
0
p
0
p
×
p
μ
2
Z
⊗
2
Z
μ
2
0
p
T
Z
T
μ
2
μ
2
)
ρ
0
(
t
,
Z
,
v
,
w
0
)
n
n
n
n
n
n
n
n
n
n
n
n
n
×
(
Z
⊗
2
0
p
×
p
0
p
0
p
×
p
μ
2
Z
⊗
2
Z
μ
2
0
p
T
Z
T
μ
2
μ
2
)
∣
W
=
w
0
}
,
where f
(
·
) is the density of W. Denote t
1
(
w
)
=
(
I
p
,
0
p
×
p
,
0
p
)
m
1
-
1
(
w
)
m
2
(
w
)  and t
2
(
w
)
= (
0
p
×
p
,
0
p
×
p
,
1
p
) ∫
0
w
‍
m
1
-
1
(
w
0
)
m
2
(
w
0
)
d
w
0, where
(12)
m
1
(
w
0
)
=
∫
0
τ
‍
{
φ
2,0
(
t
,
w
0
)
-
φ
1,0
⊗
2
(
t
,
w
0
)
φ
0,0
(
t
,
w
0
)
}
d
t
,
m
2
(
w
0
)
=
∫
0
τ
‍
{
φ
1,1
(
t
,
w
0
)
-
φ
1,0
(
t
,
w
0
)
φ
0,0
(
t
,
w
0
)
φ
0,1
(
t
,
w
0
)
}
d
t
.
For k
=
0,1
,
2, l
=
0,1, let κ
k
,
l
(
t
)
=
E
[
ρ
l
(
t
,
Z
,
V
,
W
)
{
V
-
Z
t
1
T
(
W
)
-
t
2
T
(
W
)
}
⊗
k
], and V
~
i
(
W
i
)
=
V
i
T
-
Z
i
T
t
1
(
W
i
)
-
t
2
(
W
i
),  i
=
1
,
…
,
n
.
∥
a
∥
=
(
a
T
a
)
1
/
2 for any vector a.Theorem 1.
If Assumptions (A.i)–(A.vi) are satisfied andb
n
→
0 as n
→
∞, then there exist local minimizers α
^ of L
n
*
(
α
) such that ∥
α
^
-
α
0
∥
=
O
p
(
n
-
1
/
2
+
a
n
).Remark 2.
Theorem1 indicates that by choosing a proper λ
n, which leads to a
n
=
O
(
n
-
1
/
2
), there exists a root-n consistent penalized least squares estimator of α.Theorem 3.
Assume that(13)
∞
>
lim
⁡
inf
⁡
n
→
∞
lim
⁡
inf
⁡
θ
→
0
+
p
λ
n
′

(
θ
)
λ
n
>
0
.
Under Assumptions (A.i)–(A.vi), if λ
n
→
0 and n
λ
n
→
∞ as n
→
∞, then the local root-n consistent local minimizer α
^
=
(
α
^
1
T
,
α
^
2
T
)
T in Theorem 1 with probability tending to 1 satisfies(1)
α
^
2
=
0; and
(2)
asymptotic normality(14)
n
(
I
1
(
τ
)
+
Γ
)
×
{
α
^
1
-
α
10
-
(
I
1
(
τ
)
+
Γ
)
-
1
b
}
⟶
D
N
(
0
,
Σ
1
(
τ
)
)
,where I
1
(
τ
) and Σ
1
(
τ
) are the first s
×
s submatrix of I
(
τ
)n and Σ
(
τ
), respectively, and
(15)
b
=
(
p
λ
n
′
(
|
α
10
|
)
sgn
⁡
(
α
10
)
,
…
,
p
λ
n
′
(
|
α
s
0
|
)
sgn
⁡
(
α
s
0
)
)
T
,
Γ
=
diag
⁡
(
p
λ
n
′′
(
|
α
10
|
)
,
…
,
p
λ
n
′′
(
|
α
s
0
|
)
)
.


## 3. Variable Selection for the Varying Coefficients

When some variables ofZ are not relevant in the additive hazard regression model, the corresponding functional coefficients are zero. In this section, by using basis function expansion techniques, we try to estimate those zero functional coefficients as identically zero via nonconcave penalty functions.For a simplified expression, we denoteg
(
·
) as β
0
(
·
) and Z
0
≡
1. Hence, we can rewrite the additive hazard regression model (1) as
(16)
λ
(
t
∣
Z
i
*
,
V
i
,
W
i
)
=
λ
0
(
t
)
+
β
*
T
(
W
i
)
Z
i
*
(
t
)
+
α
T
V
i
(
t
)
,
where Z
i
*
=
(
1
,
Z
i
T
)
T and β
*
(
W
i
)
=
(
β
0
(
W
i
)
,
β
T
(
W
i
)
)
T. Assume α
* is a given root-n consistent estimator of α. From the arguments in Section 2, such an α
* is available.


### 3.1. Penalized Function

We approximate each element ofβ
*
(
w
) by the basis expansion; that is, β
k
(
w
)
≈
∑
l
=
1
L
k
‍
γ
k
l
B
k
l
(
w
) for k
=
0,1
,
…
,
p. This approximation indicates that the hazard function given in (16) can be approximated as follows:
(17)
λ
(
t
∣
Z
i
*
,
V
i
,
W
i
)
≈
λ
0
(
t
)
λ
+
∑
k
=
0
p
‍
Z
i
k
{
∑
l
=
1
L
k
‍
γ
k
l
B
k
l
(
w
)
}
+
α
T
V
i
.
Let γ
k
=
(
γ
k
1
,
…
,
γ
k
l
k
)
T and γ
=
(
γ
0
T
,
γ
1
T
,
…
,
γ
p
T
)
T. Write (18)
B
(
W
i
)
=
(
B
01
(
W
i
)
⋯
B
0
L
0
(
W
i
)
0
⋯
0
⋯
0
⋯
0
0
⋯
0
B
11
(
W
i
)
⋯
B
1
L
1
(
W
i
)
⋮
0
⋯
0
⋮
⋮
⋱
⋮
0
⋯
0
0
⋯
0
⋯
B
p
1
(
W
i
)
⋯
B
p
L
p
(
W
i
)
)
,and U
i
(
W
i
)
=
Z
i
*
T
B
(
W
i
). Then, by using similar arguments as in Section 2, we estimate the basic hazard function Λ
0
(
t
) by
(19)
Λ
~
0
(
t
)
=
∑
i
=
1
n
‍
∫
0
T
‍
d
N
i
(
u
)
-
Y
i
(
u
)
{
U
i
(
W
i
)
γ
+
α
*
T
V
i
}
d
u
∑
i
=
1
n
Y
i
(
u
)
and estimate γ by minimizing the following least square function with respect to γ:
(20)
L
~
(
γ
)
n
=
1
2
∑
i
=
1
n
‍
∫
0
τ
‍
[
{
λ
~
0
(
t
)
+
U
i
(
W
i
)
γ
+
α
*
T
V
i
}
d
t
d
N
i
(
t
)
-
Y
i
(
t
)
n
n
n
n
n
n
n
n
n
n
n
n
×
{
λ
~
0
(
t
)
+
U
i
(
W
i
)
γ
+
α
*
T
V
i
}
d
t
]
2
.
Suppose that β
k
=
0,  k
=
q
+
1
,
…
,
p, as mentioned before. We are trying to correctly identify these zero functional coefficients via nonconcave penalty functions, that is, via minimizing the following penalized least square function with respect to γ:
(21)
P
L
n
(
γ
)
=
L
~
(
γ
)
n
+
n
∑
k
=
0
p
‍
p
λ
(
∥
γ
k
∥
k
)
,
where ∥
γ
k
∥
k
2
=
γ
k
T
R
k
γ
k and R
k
=
(
r
k
i
j
)
L
k
×
L
k with r
k
i
j
=
∫
‍
B
k
i
(
w
)
B
k
j
(
w
)
d
w.Letγ
^ be the minimizer of penalized least square function (21) and define β
^
*
(
w
)
=
B
(
w
)
γ
^.Remark 4.
Various basis systems, including Fourier basis, polynomial basis, and B-spline basis, can be used in the basis expansion. We focus on the B-spline basis and examine asymptotic properties ofβ
^
*
(
w
).


### 3.2. Asymptotic Properties

DefineU
¯
(
t
)
=
∑
i
=
1
n
‍
Y
i
(
t
)
U
i
(
W
i
)
/
∑
i
=
1
n
‍
Y
i
(
t
) and let γ
~ be the minimizer of L
~
(
γ
)
n in (20). Hence, we can obtain γ
~ by resolving
(22)
U
¯
n
(
γ
)
=
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
W
i
)
-
U
¯
(
t
)
}
nnnnnnnnn
×
[
d
N
i
(
t
)
-
Y
i
(
t
)
{
U
i
T
(
W
i
)
γ
+
α
*
T
V
i
}
]
d
t
=
0
.
U
¯
n
(
γ
) is the derivative of L
~
(
γ
)
n with respect to γ, the usual score function of γ without the penalty.For any square integrable functionsb
(
w
) and c
(
w
) on W, let ∥
b
(
w
)
∥
L
2 denote the L
2 norm of b
(
w
) on W, and define ∥
b
(
w
)
-
c
(
w
)
∥
L
∞
=
sup
⁡
w
∈
W
|
b
(
w
)
-
c
(
w
)
| as the L
∞ distance between b
(
w
) and c
(
w
). Let G
k denote all functions that have the form ∑
l
=
1
L
k
‍
γ
k
l
B
k
l
(
w
),  k
=
0,1
,
…
,
p for B-spline base. L
max
⁡
=
max
⁡
0
≤
k
≤
p
L
k and ρ
n
=
max
⁡
0
≤
k
≤
p
inf
⁡
c
∈
G
k
∥
β
k
-
c
∥
L
∞. Denote β
~
*
(
w
)
=
B
(
w
)
γ
~.Leta
n
*
=
max
⁡
{
p
λ
n
′
(
∥
β
k
∥
L
2
)
:
∥
β
k
∥
L
2
≠
0
} and b
n
*
=
max
⁡
{
p
λ
n
′′
(
∥
β
k
∥
L
2
)
:
∥
β
k
∥
L
2
≠
0
}.Theorem 5.
Suppose that (13) and Assumptions (A.i) and (A.vii)–(A.iv) hold. If λ
n
→
0 and λ
n
/
max
⁡
{
ρ
n
,
a
n
*
,
(
L
max
⁡
/
n
)
1
/
2
}
→
∞ as n
→
∞, then one has(a)
β
^
k
=
0,  k
=
q
+
1
,
…
,
p, with probability approaching 1; and
(b)
∥
β
^
k
-
β
k
∥
L
2
=
O
p
(
max
⁡
{
ρ
n
,
a
n
*
,
(
L
max
⁡
/
n
)
1
/
2
}
),  k
=
0
,
…
,
q.Remark 6.
LetG
k be a space of splines with a degree no less than 1 and with L
k equally spaced interior knots, where L
k
≍
n
1
/
5,  k
=
0,1
,
…
,
p. Note that ρ
n
=
O
(
L
max
⁡
-
2
) [20, Theorem 6.21]. Hence, if a
n
*
=
o
p
(
ρ
n
), Theorem 5 implies that ∥
β
^
k
-
β
k
∥
L
2
=
O
p
(
n
-
2
/
5
),  k
=
0
,
…
,
q.Remark 7.
Suppose thatL
k
≍
n
1
/
5,  k
=
0,1
,
…
,
p. Because a
n
*
=
max
⁡
{
p
λ
n
′
(
∥
β
k
∥
)
:
∥
β
k
∥
≠
0
}, which implies that for the hard thresholding and SCAD penalty function, the convergent rate of penalized least square estimator is n
-
2
/
5. This is the optimal rate for nonparametric regression [21] if λ
n
→
0. However, for the L
1 penalty, a
n
*
=
λ
n; hence λ
n
=
o
p
(
ρ
n
) is required to lead to the optimal rate. On the other hand, the oracle property in part (b) of Theorem 5 requires λ
n
/
ρ
n
→
∞, which contradicts λ
n
=
o
p
(
ρ
n
). As a consequence, for the L
1 penalty, the oracle property does not hold.Next, we demonstrate the asymptotic normality ofβ
^
*. First, we analyze β
^
(
1
)
=
(
β
^
0
,
…
,
β
^
q
)
T, that is, the penalized estimator of β
(
1
)
=
(
β
0
,
…
,
β
q
)
T. Let Z
i
(
1
), U
i
(
1
), and U
¯
(
1
) denote the selected columns of Z
i
*, U
i, and U
¯, respectively, corresponding to β
(
1
). Similarly, let B
(
1
)
(
W
i
) denote selected diagonal blocks of B
(
W
i
),  i
=
1
,
…
,
n. Define γ
(
1
)
=
(
γ
0
T
,
γ
1
T
,
…
,
γ
q
T
)
T and γ
^
(
1
)
=
(
γ
^
0
T
,
γ
^
1
T
,
…
,
γ
^
q
T
)
T. Then by part (a) of Theorem 5 and (21), γ
^
(
1
) is the local solution of
(23)
P
U
n
*
(
γ
(
1
)
)
=
-
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
(
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
)
T
n
n
n
n
n
n
n
n
n
n
n
n
×
[
{
U
i
(
1
)
(
W
i
)
γ
(
1
)
+
α
*
T
V
i
}
d
N
i
(
t
)
-
Y
i
(
t
)
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
×
{
U
i
(
1
)
(
W
i
)
γ
(
1
)
+
α
*
T
V
i
}
]
d
t
+
∑
k
=
0
q
‍
∂
p
λ
(
∥
γ
k
∥
k
)
∂
γ
(
1
)
.
Recall that α
* is root-n consistent. It follows from (23) that
(24)
0
=
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
d
M
i
(
t
)
nnnnnnnn
×
(
1
+
O
p
(
n
-
1
/
2
)
)
-
∑
k
=
0
q
‍
∂
p
λ
(
∥
γ
^
k
∥
k
)
∂
γ
(
1
)
-
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
Y
i
(
t
)
nnnnnnnn
×
{
U
i
(
1
)
(
W
i
)
γ
^
(
1
)
-
Z
i
(
1
)
T
β
(
1
)
(
W
i
)
}
d
t
=
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
d
M
i
(
t
)
nnnnnnn
×
(
1
+
O
p
(
n
-
1
/
2
)
)
-
∑
k
=
0
q
‍
∂
p
λ
(
∥
γ
^
k
∥
k
)
∂
γ
(
1
)
-
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
Z
i
(
1
)
T
Y
i
(
t
)
nnnnnnnnn
×
{
β
^
(
1
)
(
W
i
)
-
β
(
1
)
(
W
i
)
}
d
t
.Thus we can obtain the following theorem.Theorem 8.
Suppose Assumptions (A.i) and (A.vii)–(A.iv) hold,lim
⁡
n
→
∞
⁡
ρ
n
=
0, if λ
n
→
0 and λ
n
/
max
⁡
{
ρ
n
,
a
n
*
,
(
L
max
⁡
/
n
)
1
/
2
}
→
∞ as n
→
∞. For any w
∈
W, let γ
¯
(
1
) and β
¯
(
1
)
(
w
) be the conditional means of γ
^
(
1
) and β
^
(
1
)
(
w
) given {
Z
i
,
V
i
,
W
i
,
i
=
1,2
,
…
,
n
}. We then have
(25)
n
L
max
⁡
-
1
{
B
(
1
)
(
w
)
(
Ξ
1
+
Γ
*
)
-
1
Ω
*
(
Ξ
1
+
Γ
*
)
-
1
{
B
(
1
)
(
w
)
}
T
}
-
1
/
2
×
{
β
^
(
1
)
(
w
)
-
β
¯
(
1
)
(
w
)
}
⟶
L
N
(
0
,
I
q
+
1
)
,
where Γ
*
=
diag
⁡
(
(
L
max
⁡
/
2
)
(
p
λ
′
(
∥
β
k
∥
L
2
)
/
∥
β
k
∥
L
2
)
R
k
)
0
≤
k
≤
q, Ξ
1
=
lim
⁡
n
→
∞
⁡
(
L
max
⁡
/
n
)
∑
i
=
1
n
∫
0
τ
‍
{
U
i
(
1
)
(
w
)
-
U
¯
(
1
)
(
t
)
}
⊗
2 Y
i
(
t
)
d
t, and Ω
*
=
lim
⁡
n
→
∞
⁡
(
L
max
⁡
/
n
) ∑
i
=
1
n
∫
0
τ
‍
U
i
(
1
)
(
W
i
) -
U
¯
(
1
)
(
t
)
}
⊗
2
Y
i
(
t
)
λ
i
(
t
)
d
t.


## 3.1. Penalized Function

We approximate each element ofβ
*
(
w
) by the basis expansion; that is, β
k
(
w
)
≈
∑
l
=
1
L
k
‍
γ
k
l
B
k
l
(
w
) for k
=
0,1
,
…
,
p. This approximation indicates that the hazard function given in (16) can be approximated as follows:
(17)
λ
(
t
∣
Z
i
*
,
V
i
,
W
i
)
≈
λ
0
(
t
)
λ
+
∑
k
=
0
p
‍
Z
i
k
{
∑
l
=
1
L
k
‍
γ
k
l
B
k
l
(
w
)
}
+
α
T
V
i
.
Let γ
k
=
(
γ
k
1
,
…
,
γ
k
l
k
)
T and γ
=
(
γ
0
T
,
γ
1
T
,
…
,
γ
p
T
)
T. Write (18)
B
(
W
i
)
=
(
B
01
(
W
i
)
⋯
B
0
L
0
(
W
i
)
0
⋯
0
⋯
0
⋯
0
0
⋯
0
B
11
(
W
i
)
⋯
B
1
L
1
(
W
i
)
⋮
0
⋯
0
⋮
⋮
⋱
⋮
0
⋯
0
0
⋯
0
⋯
B
p
1
(
W
i
)
⋯
B
p
L
p
(
W
i
)
)
,and U
i
(
W
i
)
=
Z
i
*
T
B
(
W
i
). Then, by using similar arguments as in Section 2, we estimate the basic hazard function Λ
0
(
t
) by
(19)
Λ
~
0
(
t
)
=
∑
i
=
1
n
‍
∫
0
T
‍
d
N
i
(
u
)
-
Y
i
(
u
)
{
U
i
(
W
i
)
γ
+
α
*
T
V
i
}
d
u
∑
i
=
1
n
Y
i
(
u
)
and estimate γ by minimizing the following least square function with respect to γ:
(20)
L
~
(
γ
)
n
=
1
2
∑
i
=
1
n
‍
∫
0
τ
‍
[
{
λ
~
0
(
t
)
+
U
i
(
W
i
)
γ
+
α
*
T
V
i
}
d
t
d
N
i
(
t
)
-
Y
i
(
t
)
n
n
n
n
n
n
n
n
n
n
n
n
×
{
λ
~
0
(
t
)
+
U
i
(
W
i
)
γ
+
α
*
T
V
i
}
d
t
]
2
.
Suppose that β
k
=
0,  k
=
q
+
1
,
…
,
p, as mentioned before. We are trying to correctly identify these zero functional coefficients via nonconcave penalty functions, that is, via minimizing the following penalized least square function with respect to γ:
(21)
P
L
n
(
γ
)
=
L
~
(
γ
)
n
+
n
∑
k
=
0
p
‍
p
λ
(
∥
γ
k
∥
k
)
,
where ∥
γ
k
∥
k
2
=
γ
k
T
R
k
γ
k and R
k
=
(
r
k
i
j
)
L
k
×
L
k with r
k
i
j
=
∫
‍
B
k
i
(
w
)
B
k
j
(
w
)
d
w.Letγ
^ be the minimizer of penalized least square function (21) and define β
^
*
(
w
)
=
B
(
w
)
γ
^.Remark 4.
Various basis systems, including Fourier basis, polynomial basis, and B-spline basis, can be used in the basis expansion. We focus on the B-spline basis and examine asymptotic properties ofβ
^
*
(
w
).


## 3.2. Asymptotic Properties

DefineU
¯
(
t
)
=
∑
i
=
1
n
‍
Y
i
(
t
)
U
i
(
W
i
)
/
∑
i
=
1
n
‍
Y
i
(
t
) and let γ
~ be the minimizer of L
~
(
γ
)
n in (20). Hence, we can obtain γ
~ by resolving
(22)
U
¯
n
(
γ
)
=
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
W
i
)
-
U
¯
(
t
)
}
nnnnnnnnn
×
[
d
N
i
(
t
)
-
Y
i
(
t
)
{
U
i
T
(
W
i
)
γ
+
α
*
T
V
i
}
]
d
t
=
0
.
U
¯
n
(
γ
) is the derivative of L
~
(
γ
)
n with respect to γ, the usual score function of γ without the penalty.For any square integrable functionsb
(
w
) and c
(
w
) on W, let ∥
b
(
w
)
∥
L
2 denote the L
2 norm of b
(
w
) on W, and define ∥
b
(
w
)
-
c
(
w
)
∥
L
∞
=
sup
⁡
w
∈
W
|
b
(
w
)
-
c
(
w
)
| as the L
∞ distance between b
(
w
) and c
(
w
). Let G
k denote all functions that have the form ∑
l
=
1
L
k
‍
γ
k
l
B
k
l
(
w
),  k
=
0,1
,
…
,
p for B-spline base. L
max
⁡
=
max
⁡
0
≤
k
≤
p
L
k and ρ
n
=
max
⁡
0
≤
k
≤
p
inf
⁡
c
∈
G
k
∥
β
k
-
c
∥
L
∞. Denote β
~
*
(
w
)
=
B
(
w
)
γ
~.Leta
n
*
=
max
⁡
{
p
λ
n
′
(
∥
β
k
∥
L
2
)
:
∥
β
k
∥
L
2
≠
0
} and b
n
*
=
max
⁡
{
p
λ
n
′′
(
∥
β
k
∥
L
2
)
:
∥
β
k
∥
L
2
≠
0
}.Theorem 5.
Suppose that (13) and Assumptions (A.i) and (A.vii)–(A.iv) hold. If λ
n
→
0 and λ
n
/
max
⁡
{
ρ
n
,
a
n
*
,
(
L
max
⁡
/
n
)
1
/
2
}
→
∞ as n
→
∞, then one has(a)
β
^
k
=
0,  k
=
q
+
1
,
…
,
p, with probability approaching 1; and
(b)
∥
β
^
k
-
β
k
∥
L
2
=
O
p
(
max
⁡
{
ρ
n
,
a
n
*
,
(
L
max
⁡
/
n
)
1
/
2
}
),  k
=
0
,
…
,
q.Remark 6.
LetG
k be a space of splines with a degree no less than 1 and with L
k equally spaced interior knots, where L
k
≍
n
1
/
5,  k
=
0,1
,
…
,
p. Note that ρ
n
=
O
(
L
max
⁡
-
2
) [20, Theorem 6.21]. Hence, if a
n
*
=
o
p
(
ρ
n
), Theorem 5 implies that ∥
β
^
k
-
β
k
∥
L
2
=
O
p
(
n
-
2
/
5
),  k
=
0
,
…
,
q.Remark 7.
Suppose thatL
k
≍
n
1
/
5,  k
=
0,1
,
…
,
p. Because a
n
*
=
max
⁡
{
p
λ
n
′
(
∥
β
k
∥
)
:
∥
β
k
∥
≠
0
}, which implies that for the hard thresholding and SCAD penalty function, the convergent rate of penalized least square estimator is n
-
2
/
5. This is the optimal rate for nonparametric regression [21] if λ
n
→
0. However, for the L
1 penalty, a
n
*
=
λ
n; hence λ
n
=
o
p
(
ρ
n
) is required to lead to the optimal rate. On the other hand, the oracle property in part (b) of Theorem 5 requires λ
n
/
ρ
n
→
∞, which contradicts λ
n
=
o
p
(
ρ
n
). As a consequence, for the L
1 penalty, the oracle property does not hold.Next, we demonstrate the asymptotic normality ofβ
^
*. First, we analyze β
^
(
1
)
=
(
β
^
0
,
…
,
β
^
q
)
T, that is, the penalized estimator of β
(
1
)
=
(
β
0
,
…
,
β
q
)
T. Let Z
i
(
1
), U
i
(
1
), and U
¯
(
1
) denote the selected columns of Z
i
*, U
i, and U
¯, respectively, corresponding to β
(
1
). Similarly, let B
(
1
)
(
W
i
) denote selected diagonal blocks of B
(
W
i
),  i
=
1
,
…
,
n. Define γ
(
1
)
=
(
γ
0
T
,
γ
1
T
,
…
,
γ
q
T
)
T and γ
^
(
1
)
=
(
γ
^
0
T
,
γ
^
1
T
,
…
,
γ
^
q
T
)
T. Then by part (a) of Theorem 5 and (21), γ
^
(
1
) is the local solution of
(23)
P
U
n
*
(
γ
(
1
)
)
=
-
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
(
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
)
T
n
n
n
n
n
n
n
n
n
n
n
n
×
[
{
U
i
(
1
)
(
W
i
)
γ
(
1
)
+
α
*
T
V
i
}
d
N
i
(
t
)
-
Y
i
(
t
)
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
×
{
U
i
(
1
)
(
W
i
)
γ
(
1
)
+
α
*
T
V
i
}
]
d
t
+
∑
k
=
0
q
‍
∂
p
λ
(
∥
γ
k
∥
k
)
∂
γ
(
1
)
.
Recall that α
* is root-n consistent. It follows from (23) that
(24)
0
=
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
d
M
i
(
t
)
nnnnnnnn
×
(
1
+
O
p
(
n
-
1
/
2
)
)
-
∑
k
=
0
q
‍
∂
p
λ
(
∥
γ
^
k
∥
k
)
∂
γ
(
1
)
-
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
Y
i
(
t
)
nnnnnnnn
×
{
U
i
(
1
)
(
W
i
)
γ
^
(
1
)
-
Z
i
(
1
)
T
β
(
1
)
(
W
i
)
}
d
t
=
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
d
M
i
(
t
)
nnnnnnn
×
(
1
+
O
p
(
n
-
1
/
2
)
)
-
∑
k
=
0
q
‍
∂
p
λ
(
∥
γ
^
k
∥
k
)
∂
γ
(
1
)
-
1
n
∑
i
=
1
n
‍
∫
0
τ
‍
{
U
i
(
1
)
(
W
i
)
-
U
¯
(
1
)
(
t
)
}
T
Z
i
(
1
)
T
Y
i
(
t
)
nnnnnnnnn
×
{
β
^
(
1
)
(
W
i
)
-
β
(
1
)
(
W
i
)
}
d
t
.Thus we can obtain the following theorem.Theorem 8.
Suppose Assumptions (A.i) and (A.vii)–(A.iv) hold,lim
⁡
n
→
∞
⁡
ρ
n
=
0, if λ
n
→
0 and λ
n
/
max
⁡
{
ρ
n
,
a
n
*
,
(
L
max
⁡
/
n
)
1
/
2
}
→
∞ as n
→
∞. For any w
∈
W, let γ
¯
(
1
) and β
¯
(
1
)
(
w
) be the conditional means of γ
^
(
1
) and β
^
(
1
)
(
w
) given {
Z
i
,
V
i
,
W
i
,
i
=
1,2
,
…
,
n
}. We then have
(25)
n
L
max
⁡
-
1
{
B
(
1
)
(
w
)
(
Ξ
1
+
Γ
*
)
-
1
Ω
*
(
Ξ
1
+
Γ
*
)
-
1
{
B
(
1
)
(
w
)
}
T
}
-
1
/
2
×
{
β
^
(
1
)
(
w
)
-
β
¯
(
1
)
(
w
)
}
⟶
L
N
(
0
,
I
q
+
1
)
,
where Γ
*
=
diag
⁡
(
(
L
max
⁡
/
2
)
(
p
λ
′
(
∥
β
k
∥
L
2
)
/
∥
β
k
∥
L
2
)
R
k
)
0
≤
k
≤
q, Ξ
1
=
lim
⁡
n
→
∞
⁡
(
L
max
⁡
/
n
)
∑
i
=
1
n
∫
0
τ
‍
{
U
i
(
1
)
(
w
)
-
U
¯
(
1
)
(
t
)
}
⊗
2 Y
i
(
t
)
d
t, and Ω
*
=
lim
⁡
n
→
∞
⁡
(
L
max
⁡
/
n
) ∑
i
=
1
n
∫
0
τ
‍
U
i
(
1
)
(
W
i
) -
U
¯
(
1
)
(
t
)
}
⊗
2
Y
i
(
t
)
λ
i
(
t
)
d
t.


## 4. Simulation Studies

We carried out two sets of simulation studies to examine the finite-sample properties of the proposed methods.w was generated from a uniform distribution over [
0,3
] and λ
n was selected via cross-validation [7]. We set a
=
2
+
3 and reported sampling properties based on 200 replications. We consider two scenarios: (I) the performance of the regression coefficient estimates. Without loss of generality, we set p
=
1. We generated failure times from the partially linear additive hazards model (1) with β
(
w
)
=
1.2
+
sin
(
2
w
), α
=
(
1.5,0
,
1,0
,
0
)
T, g
(
w
)
=
0.3
w
2, and λ
0
(
t
)
=
0.5. Covariate Z was generated from a uniform distribution over [
0,1
], V
1 was a Bernoulli random variable taking a value of 0 or 1 with a probability of 0.5, V
2
,
V
3 were generated from uniform distributions over [
0,1
] and [
0,2
], respectively, and (
V
4
,
V
5
)
T was generated from a normal distribution with mean zero, variance 0.25, and covariance 0.25. We generated the censoring time C from a uniform distribution over [
c
c
/
2,3
c
c
/
2
]. We took c
c
=
0.86 to yield an approximate censoring rate of CR
=
20
% and c
c
=
0.35 to yield an approximate censoring rate of CR
=
30
%. We used sample size n
=
200. We selected the optimal bandwidth h
^
opt by using the method of [22] and found that the value is approximately 0.203. We present the average number of zero coefficients in Table 1, in which the column labeled “correct” presents the average restricted only to the true zero coefficients, while the column labeled “incorrect” depicts the average of coefficients erroneously set to 0. We also report the standard deviations (SD), the average standard errors (SE), and the coverage probability (CP) of the 95% confidence interval for the nonzero parameters in Table 2.Table 1
Results of the simulation study for scenario (I). Here “correct” is the average true positive rate and “incorrect” is the average true negative rate.
SCAD
HARD
L
1
SCAD
HARD
L
1
n
=
200, CR
=
20%
n
=
200, CR
=
30%
Correct
2.73
2.76
2.67
Correct
2.85
2.85
2.85
Incorrect
0.12
0.11
0.21
Incorrect
0.13
0.13
0.15Table 2
Simulation results for scenario (I): bias, standard error (SE), standard deviation (SD), and coverage probability (CP).
α
1
α
3
SCAD
HARD
L
1
SCAD
HARD
L
1
n
=
200, CR
=
20%
bias
0.0738
0.0729
0.0527
0.0552
0.0475
0.0324
SE
0.0767
0.0767
0.0757
0.0567
0.0566
0.0558
SD
0.0635
0.0645
0.0631
0.0495
0.0501
0.0558
CP
0.96
0.91
0.91
0.90
0.89
0.88
n
=
200, CR
=
30%
bias
0.0726
0.0726
0.0402
0.0480
0.0466
0.0280
SE
0.0824
0.0824
0.0812
0.0612
0.0612
0.0602
SD
0.0701
0.0721
0.0726
0.0618
0.0602
0.0613
CP
0.93
0.89
0.92
0.94
0.90
0.91As shown in Tables1 and 2, SCAD, HARD, and L
1 perform well and select about the same correct number of significant variables. The coverage probability (CP) based on SCAD and HARD, however, is better than that based on L
1.In scenario (II), we focused on functional coefficientsβ
(
·
). We set d
=
1, g
(
w
)
=
0.3
w
2, and α
=
1. We generated failure times from the partially linear additive hazards model (1) with β
1
(
w
)
=
1.2
+
sin
(
2
w
), β
3
(
w
)
=
0.5
w
2, and β
2
(
w
)
=
β
4
(
w
)
=
β
5
(
w
)
≡
0. Both covariates V and Z
1 were generated from a uniform distribution over [
0,1
], respectively, (
V
2
,
V
3
)
T was generated from a normal distribution with a mean of (
0.5,0.5
)
T, variance (
0.25,0.25
)
T and covariance 0, V
4 was a Bernoulli random variable taking a value of 0 or 1 with a probability of 0.5, and V
5 was uniform over [
0,2
]. In simulation II, we took c
c
=
1.24 to yield an approximate censoring rate of CR
=
10
% and c
c
=
0.45 to yield an approximate censoring rate of CR
=
30
%. We used sample size n
=
200. The average number of zero functional coefficients based on SCAD and HARD is reported in Table 3, and fitted curves of nonzero functional coefficients are presented in Figures 1 and 2 based on SCAD and HARD, respectively.Table 3
Results from the simulation study for scenario (II). The legend is the same as in Table1.
SCAD
HARD
SCAD
HARD
n
=
200, CR
=
10%
n
=
200, CR
=
30%
Correct
2.61
2.79
Correct
2.625
2.745
Incorrect
0.165
0.195
Incorrect
0.48
0.48Results of the simulation study for scenario (I): the estimated curves of the nonzero functional coefficientsβ
1
(
w
), β
3
(
w
), and g
(
w
). The solid, dashed, and dotted curves represent the true function, estimated functions, and 95% pointwise confidence intervals, respectively. Left panel—HARD penalty; right panel—SCAD penalty.
(a)
HARD
(b)
SCAD
(c)
HARD
(d)
SCAD
(e)
HARD
(f)
SCADResults of the simulation study for scenario (II). The legend is the same as in Figure1.
(a)
HARD
(b)
SCAD
(c)
HARD
(d)
SCAD
(e)
HARD
(f)
SCADNotice that penalized estimates of functional coefficients are based on the results in Step 1. Because the oracle property does not hold for theL
1 penalty according to Remark 7, we did not take L
1 penalty into account.We can see from Table3 and Figures 1 and 2 that the penalized spline estimation of the functional coefficient that we proposed performed well. Furthermore, for the varying coefficient part, the “correct” numbers of significant variables selected based on SCAD and HARD penalties are close. In comparing the results between different censoring rates, it is also shown that our method is extremely robust against the censoring rate.


## 5. Real Data Example

In this section, we apply the proposed method to a SUPPORT (Study to Understand Prognoses Preferences Outcomes and Risks of Treatment) dataset [23]. This study was a multicenter study designed to examine outcomes and clinical decision making for seriously ill-hospitalized patients. One thousand patients suffering from one of eight different diseases, including acute respiratory failure (ARF), chronic obstructive pulmonary disease (COPD), congestive heart failure (CHF), cirrhosis, coma, colon cancer, lung cancer, and multiple organ system failure (MOSF) were followed up to 5.56 years. 672 observations were collected, including age, sex, race, coma score (scoma), number of comorbidities (num.co), charges, average of therapeutic intervention scoring system (avtiss), white blood cell count (wbc), heart rate (hrt), respiratory rate (resp), temperature (temp), serum bilirubin level (bili), creatinine (crea), serum sodium (sod), and imputed ADL calibrated to surrogate (adlsc). For more details about this study, refer to [23].We suppose that binary and multinary risk factors have constant coefficients, which yields a 12-dimensional parameter after transforming multinary variables into binary. Conversely, other risk factors are supposed to have varying coefficients, that is, a 13-dimension functional coefficient. Here the main exposure variableW corresponds to age. Censoring rate of 672 patients is 32.7%. The model we considered is
(26)
λ
(
t
∣
Z
i
,
V
i
,
ag
e
i
)
=
λ
0
(
t
)
+
∑
j
=
1
13
‍
β
j
(
ag
e
i
)
Z
i
j
+
∑
k
=
1
12
‍
α
k
V
i
k
+
g
(
ag
e
i
)
,
where V
i
j,  j
=
1,2
,
…
,
13, denote binary variables of the ith patient and Z
i
k
,
k
=
1,2
,
…
,
12, denote the others of the ith patient. Here, baseline hazard function λ
0
(
t
) is the failure risk of patients who are white females suffering from MOSF with other variables all equal to zero.By fitting model (26) with a SCAD penalty, the parameter of “other race” is identified as zero, and functional coefficients of num.co, charges, wbc, hrt, temp, crea, sod, adlsc, and g
(
·
) are identified as zero functions. Identified significant risk factors and results of their parametric or functional coefficients are shown in Figure 3 and Table 4.Table 4
Applied to SUPPORT data. Estimation results of significant binary variables using the SCAD penalty.
Risk factor
Estimator
Variance × 10−3
Sex
0.0006
0.0178
Disease group
ARF
0.0007
0.0247
CHF
0.0011
0.0309
Cirrhosis
0.0018
0.0450
Colon cancer
0.0034
0.0580
Lung cancer
0.0042
0.0585
Coma
0.0074
0.1410
COPD
0.0017
0.0344
Race
Asian
−0.0002
0.0414
Black
0.0006
0.0232
Hispanic
0.0008
0.0421Results of the real data example. The estimated curves of those nonzero functional coefficients. The solid and dashed curves represent the estimated curve and 95% pointwise confidence intervals (×
1
0
4).
(a)
scoma
(b)
avtisst
(c)
meanbp
(d)
resp
(e)
biliThe identified zero parameter of “other race” shows that, besides Asian, black, and Hispanic, “other race” has no significant different impact on risk in contrast with white people. Table4 demonstrates that Asians have the lowest failure probability, as do people who are suffering from MOSF. Conversely, coma is the most dangerous among these eight different diseases. Figure 3 shows the fitted coefficients and their 95% pointwise confidence intervals of the five risk factors which were identified as significant.



---

*Source: 360249-2014-05-25.xml*