\documentclass[coll,linenumbers]{imsart}


\RequirePackage[OT1]{fontenc}
%\usepackage{natbib}
\usepackage{amsthm,amsmath,natbib}
\RequirePackage[colorlinks,citecolor=blue,urlcolor=blue]{hyperref}
\RequirePackage{hypernat}
\usepackage{epsfig, amsmath, amssymb, amscd, graphicx}
\usepackage{dcolumn}
\usepackage{fancybox}
\usepackage{rotating}

\usepackage{verbatim}

\pubyear{2009}
\volume{0}
%\volumetitle{IMS Collections}
\firstpage{46}
\lastpage{65}
\doi{0018/0000.0000}
\arxiv{math.PR/0000018}



\newtheorem{thm}{Theorem}
\newtheorem{lem}{Lemma}

\newtheorem{cor}{Corollary}
\newcommand{\eps}{\varepsilon}
\newcommand{\To}{\longrightarrow}
\newcommand{\h}{\mathcal{H}}
\newcommand{\s}{\mathcal{S}}
\newcommand{\A}{\mathcal{A}}
\newcommand{\J}{\mathcal{J}}
\newcommand{\M}{\mathcal{M}}
\newcommand{\W}{\mathcal{W}}
\newcommand{\X}{\mathcal{X}}
\newcommand{\R}{\mathcal{R}}
\newcommand{\I}{\mathcal{I}}
\newcommand{\B}{\bfseries}
\newcommand{\MB}{\mathbf}
\newcommand{\MI}{\mathit}
\newcommand{\littlel}{\mathcal{l}}
\newcommand{\BOP}{\mathbf{B}}
\newcommand{\BH}{\mathbf{B}(\mathcal{H})}
\newcommand{\KH}{\mathcal{K}(\mathcal{H})}
\newcommand{\eT}{\eta_0}
\newcommand{\eh}{\hat{\eta}}
\newcommand{\et}{\tilde{\eta}}
\newcommand{\bZ}{\beta^T Z}
\newcommand{\bz}{\beta^T z}
\newcommand{\bT}{\beta_0}
\newcommand{\bTZ}{{\beta_0}^T Z}
\newcommand{\bTz}{{\beta_0}^T z}
\newcommand{\bh}{\hat{\beta}}
\newcommand{\bs}{\beta^*}
\newcommand{\bt}{\tilde{\beta}}
\newcommand{\bhZ}{\hat{\beta}^T Z}
\newcommand{\bsZ}{{\beta^*}^T Z}
\newcommand{\bhz}{\hat{\beta}^T z}
\newcommand{\btZ}{\tilde{\beta}^T Z}
\newcommand{\btz}{\tilde{\beta}^T z}
\newcommand{\gh}{\hat{\gamma}}
\newcommand{\gT}{\gamma_0}
\newcommand{\gt}{\tilde{\gamma}}
\newcommand{\mh}{\hat{m}}
\newcommand{\psih}{\hat{\psi}}
\newcommand{\psit}{\tilde{\psi}}
\newcommand{\lh}{\hat{\MI l_n}}
\newcommand{\lhs}{{\MI {\hat l}_n}^*}
\newcommand{\cp}{\rightarrow_p}
\newcommand{\cd}{\rightarrow_{D}}
\newcommand{\Real}{\mathbf{R}}
\newcommand{\Complex}{\mathbf{C}}
\newcommand{\Field}{\mathbf{F}}
\newcommand{\RPlus}{\Real^{+}}
\newcommand{\Polar}{\mathcal{P}_{\s}}
\newcommand{\Poly}{\mathcal{P}(E)}
\newcommand{\EssD}{\mathcal{D}}
\newcommand{\Lom}{\mathcal{L}}
\newcommand{\States}{\mathcal{T}}
\newcommand{\abs}[1]{\left\vert#1\right\vert}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\seq}[1]{\left<#1\right>}
\newcommand{\norm}[1]{\left\Vert#1\right\Vert}
\newcommand{\essnorm}[1]{\norm{#1}_{\ess}}

%\setlength{\headsep}{8mm} \setlength{\topmargin}{-0.75in}
%\setlength{\textheight}{9.7in} \setlength{\textwidth}{6.7in}
%\setlength{\evensidemargin}{-0.3in}
%\setlength{\oddsidemargin}{-0.3in}
%\renewcommand{\baselinestretch}{1.64}
\newtheorem{proofof}{Proof of Theorem}[section]
\newcommand{\findemo}{\hfill \, \rule{1ex}{1ex}\,}
\endlocaldefs

\begin{document}
\begin{frontmatter}
\title{Proportional Hazards Regression with Unknown Link
Function}
%\protect\thanksref{T1}}
\runtitle{Proportional Hazards Regression}

\begin{aug}
\author{\fnms{Wei}
\snm{Wang}\thanksref{a1}\ead[label=e1]{wwang@partners.org}},
\author{\fnms{Jane-Ling}
\snm{Wang}\thanksref{a2}\ead[label=e2]{wang@wald.ucdavis.edu}}
\and
\author{\fnms{Qihua}
\snm{Wang}\thanksref{a3}\ead[label=e3]{qhwang@amss.ac.cn}}
\address[a1]{Harvard Medical School and Brigham and Women's Hospital, Boston, USA, email: 
	\url{wwang4@partners.org}}
\address[a2]{University of California, Davis, USA, Apartado Postal 402, email: 
	\url{wang@wald.ucdavis.edu}}
\address[a3]{Chinese Academy of Science, Beijing, China, 
email: \url{qhwang@amss.ac.cn} and University of Hongkong, Hongkong, 
email: \url{qhwang@hku.hk}}

\runauthor{W. Wang, J. Wang and Q. Wang}
\contributor{Wang, W.}{Harvard Medical School and Women's Hospital}
\contributor{Wang, J.}{University of California, Davis}
\contributor{Wang, Q.}{Chinese Academy of Science, Beijing, China and University of Hongkong, Hongkong}

\affiliation{Harvard Medical School and Brigham and Women's Hospital; 
University of California, Davis; Chinese Academy of Science and 
The University of Hong Kong}


\end{aug}

\begin{abstract}
Proportional hazards regression model assumes that the covariates
affect the hazard function through a link function and an index
which is a linear function of the covariates. Traditional
approaches, such as the Cox proportional hazards model, focus on
estimating the unknown index by assuming a known link function
between the log-hazard function and covariates. A linear link
function is often employed for convenience without any validation.
This paper provides an approach to estimate the link function, which
can then be used to guide the choice of a proper parametric link
function. This is accomplished through a two-step algorithm to
estimate the link function and the effects of the covariates
iteratively without involving the baseline hazard estimate. The link
function is estimated by a smoothing method based on a local version
of partial likelihood, and the index function is then estimated
using a full version of partial likelihood. Asymptotic properties of
the non-parametric link function estimate are derived, which
facilitates model checking of the adequacy of the Cox Proportional
hazards model. The approach is illustrated through a survival data
and simulations.
\end{abstract}

\begin{keyword}[class=AMS]
\kwd[Primary ]{62G05}
%\kwd{60K35}
\kwd[; secondary ]{62N02}
\end{keyword}

\begin{keyword}
\kwd{Partial likelihood}
\kwd{local partial likelihood}
\kwd{Nonparametric smoothing}
\kwd{Dimension reduction}
\end{keyword}
\tableofcontents
\end{frontmatter}


\section{Introduction}\label{sec1}

Proportional hazards regression model has played a pivotal role in
survival analysis since Cox proposed it in 1972. Let $T$ represent
survival time and $Z$ its associate covariate vector. Under the
proportional hazards model, the hazard function for $T$, given a
particular value $z$ for the covariate $Z$, is defined as
\begin{equation}\label{PHM}
\lambda\{t\mid z\}=\lambda_0(t) \exp\bigl\{\psi(\bTz)\bigr\},
\end{equation}
where $\lambda_0(t)$ is an unknown baseline hazard function
corresponding to $z=(0,\cdots,0)$, and $\psi(\cdot)$ is called the
link function with $\psi(0)=0$. With fully specified link function
$\psi$, the partial likelihood method was introduced in
\cite{Cox:1972, Cox:1975} to estimate the regression parameters,
$\beta_0$, with the option to  accommodate censored data. The most
common choice for $\psi$ is the identity function, which corresponds
to the time-honored Cox model. In reality the link function is
unknown and needs to be estimated. This is especially useful to
validate a preferred choice,  as an erroneous link function could
dramatically distort risk assessment or interpretation of odds
ratios. When the link function is known, such as in the Cox model,
model (\ref{PHM}) is a special case of the transformation model
first proposed in \cite{Doksum:1987} and subsequently studied in
\cite{DD:1988}, \cite{CWY:1995} etc. Our goal in this paper is to
consider model (\ref{PHM}) with an unknown link function. This
problem was first studied in an unpublished Ph.D. thesis
\cite{Wang:2001}. However, the procedure there was less efficient
and we propose an improved estimate, studying  its asymptotic
properties.

\iffalse%%%%%%%%%%%%%%%%%

Several alternative approaches have been explored to
reflect the nature of covariate effects. Hastie and Tibshirani
\cite{HT:1990} proposed the "Generalized Additive Model" (GAM) which
features an additive term of some unspecified smooth functions of
the covariates in place of the usual log-linear term of the
covariates. They estimated the unspecified smooth functions using
scatterplot smoothers.  \cite{Gray:1992} applied spline smoothing to
an extended GAM with applications to breast cancer prognosis. These
approaches assume different models from (\ref{PHM}), which is a
projection pursuit type of model with a single index. Model
(\ref{PHM}) has the advantage that it avoids the need to estimate
$q$ different link functions as in GAM.

is a projection pursuit type of model with a single index.  Other
single-index projection pursuit types of models include the
transformation models by \cite{Doksum:1987}, \cite{DD:1988}, and
\cite{CWY:1995}.  \cite{Doksum:1987} considered the Cox proportional
hazards model with time independent covariates as a special case of
transformation models and used a resampling scheme to compute the
maximum partial likelihood estimates of the parameters in the single
index. Dabrowska and Doksum \cite{DD:1988} also studied the
properties of partial likelihood methods in transformation models
with censored-data. The transformation model in \cite{CWY:1995}
assumed a known error distribution, and the parameters of the single
index were estimated through a class of estimating functions related
to ranks.

\fi%%%%%%%%%%%%%%%%%%%%%


Previous work focuses on the special case when the covariate is
one-dimensional, or equivalently when $\beta$ is known in
(\ref{PHM}).  Under this special one-dimensional case, a local
partial likelihood technique in \cite{TH:1987} and a variation of
the local scoring algorithm of \cite{HT:1986} can be used to
estimate the unknown link function in (\ref{PHM}). Gentleman and
Crowley \cite{GC:1991} proposed a local version of the full
likelihood instead of partial likelihood by alternating between
estimating the baseline hazard function and estimating the covariate
effects. The local likelihood methods in these papers were based on
data whose covariate values fall in a neighborhood of the targeted
location. Fan, Gijbels and King \cite{FGK:1997} used instead a local
polynomial method to approximate the local partial likelihood, and
derived rigorous asymptotic results for their approach. Simulation
studies there showed that the local partial likelihood method is
comparable to the local full likelihood method in \cite{GC:1991}.
Two spline approaches have also been considered, with smoothing
splines resulting  from a penalized partial likelihood in O'Sullivan
\cite{OS:1988} and regression splines from
 \cite{SH:1990}.

While the aforementioned approaches can be easily extended to
$q$-dimensional covariates by estimating a multivariate unknown link
function $\psi(z_1,\cdots,z_q)$, such  nonparametric approaches
 are subject to the curse of
dimensionality and may not be suitable for $q\ge 3$. Moreover, the
resulting model would be different from model (\ref{PHM}), which has
the attractive dimension reduction feature that the covariate
information is succinctly summarized in a single index and is a
nonparametric extension of the Cox proportional hazards model. Model
(\ref{PHM}) could also be used as an exploratory tool to guide the
choice of a suitable parametric link function.

A two-step iterative algorithm to estimate the link function and the
covariate effects is proposed in Section 2. In the first step, an
initial estimate of the regression parameter $\beta$ is plugged  in
model (\ref{PHM}) so that the link function can be estimated by a
smoothing method based on a local version of partial likelihood
(\cite{Cox:1972, Cox:1975}). The second step involves updating the
regression parameters using the full partial likelihood with the
estimated link function in step 1 inserted. These two steps will be
iterated until the algorithm converges. Asymptotic results for the
link estimators are stated in Section 2. In particular,  Theorem 2
provides the building blocks to check the link function and
inference for the individual risk, $\psi(\bz)$. It also reveals that
 the nonparametric estimate of the link function is as efficient as the one for model (1)
but with a known regression parameter $\beta$. Thus, there is no
efficiency loss to estimate the link function even if $\beta$ is
unknown in our setting. This is also reflected in the simulation
studies in Section 3. The approach in Section 2 is further
illustrated in Section 4 through a data set from the Worcester Heart
Attach Study. All the proofs of the main results are relegated to an
appendix.

\iffalse%%%%%%%%%%%%%%%%%%%%%%

While an estimated link function is of interest to correctly
reflect the risk associated with a covariate, a parametric link
function is often preferable to a nonparametric one to lend a
parsimonious model with more interpretable results. Thus, a main
incentive to estimate the link function is for exploratory model
checking so that a proper parametric link function will emerge as
the choice of link function to be employed in the proportional
hazards model (1).  In the German Breast Cancer Study Data
example in Section 4, we illustrate this by checking the
appropriateness of the conventional identity link function. This
data has been analyzed in the book by
\cite{HL:2008}. Our goal
is to check whether the identity link function employed there
fits the data, as this has not been previously explored. Initial
finding suggests that an identity link function seems appropriate
for the bulk of the data, but a convex link function may result in
some improvement.  Possible alternative parametric link functions
are suggested.  Section 5 contains a short discussion. All the
proofs of the main results are relegated in the appendix.

\fi%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

We remark here that \cite{HL:2006} also studied model (1) with a
different approach. They assumed that the link function is in a
finite dimensional subspace spanned by polynomial spline bases
functions and the dimension of this subspace is known. This leads to
a flexible parametric model where the spline coefficients
corresponding to the link functions and $\beta$ can then be
estimated directly through traditional partial likelihood
approaches. While this has the benefit of simplicity as everything
is in the parametric framework, it tends to underestimate the
standard errors of the estimates. Two sources of bias arise, one
derives from the fact that in reality the number of spline bases
depends on the data and is a function of the sample size, so the
standard errors are underestimated by the simple parametric
inference. In addition, the link estimation might be biased, as in
theory an infinite number of spine bases might be required to span
the unknown link function.  These biases could significantly affect
the asymptotic results.  In contrast, our approach provides correct
asymptotic theory and efficient estimation of the link function.

\section{Estimation Procedure and Main Results}

Since there are three different unknown parameters,
$\lambda_0(\cdot)$, $\psi(\cdot)$ and $\beta$ in model (\ref{PHM}),
we need to impose some conditions to ensure  identifiability.  To
identify $\lambda_0$, it suffices to set $\psi(v)=0$ at some point
$v$, a common choice is  $v=0$. Since only the direction of $\beta$
is identifiable if $\psi$ is unknown, we assume that $\|\beta\|$=1
(here $\|\cdot\|$ represents the Euclidean norm) and that the sign
of the first component of $\beta$ is positive. As for the sampling
plan, we assume an independent censoring scheme, in which the
survival time T and censoring time C are conditionally independent,
given the covariate vector Z. Let $X=min(T,C)$ be the observed
event-time and $\Delta=I\{T\le C\}$ be the censoring indicator. The
data $\{X_i, Z_i, \delta_i\}$ is an i.i.d. sample of $\{X, Z,
\Delta\}$. We use the notation  $t_i<\cdots<t_N$ to denote the $N$
distinctive ordered failure times, and (j) to denote the label of
the item failing at time $t_j$. The risk set at time $t_j$ is
denoted by  $\R_j=\{i: X_i\ge t_j\}$.

For a fixed parametric value $\beta$, one can estimate the link
function $\psi(\cdot)$ by any smoothing method, such as those cited
in Section 1 when $\beta$ is assumed known.  We adopt the local
partial likelihood approach in  \cite{FGK:1997} and assume, for a
given point $v$,  that the $p$-th order derivative of $\psi(v)$ at
point $v$ exists. A  Taylor expansion for $\beta^T Z$ in a
neighborhood of $v$ then yields,
%\begin{displaymath}
\begin{equation}\label{Taylor}
\begin{aligned}
\psi(\bZ) \approx& \
\psi(v)+\psi^{'}(v)(\bZ-v)+\cdots+{\frac{\psi^{(p)}(v)}{p!}}(\bZ-v)^P\\
=& \ \psi(v)+(\beta^T \MB Z)^T\gamma(v),
\end{aligned}
%\end{displaymath}
\end{equation}
where $\gamma(v)=\{\psi'(v),\cdots,\psi^{(p)}(v)/p!\}^T$ is   the
$p-$dimensional vector associated with the derivatives of $\psi$ and
${\beta^T \MB Z}=\{\beta^T Z-v,\cdots,(\beta^T Z-v)^p\}^T$.

Let $K$ be a kernel function, $h$ be a bandwidth, and define
$K_h(u)=h^{-1}K(u/h).$ Applying  kernel weights to the logarithm of
the global partial likelihood
\begin{displaymath}\label{globalLogL}
\sum_{j=1}^N \psi(\bZ_{(j)}) - \log\Biggl\{\sum_{i\in \R_j} \exp
\Bigl\{\psi(\bZ_{(j)})\Bigr\}\Biggr\}
\end{displaymath}
and replacing  $\psi(\bz)$ by the local approximation in
(\ref{Taylor}), we arrive at (similarly to \cite{FGK:1997})   the
 local version
of the log partial likelihood:
\begin{equation}\label{LocalLogPL}
\begin{aligned}[b]
\sum_{j=1}^NK_h\bigl(\bZ_{(j)}-v\bigr) &\Biggl[ \bigl( \beta^T {\MB
Z}_{(j)}\bigr)^T\gamma(v)\\
&- \log \biggl\{\sum_{i\in
\R_j}\exp\Bigl\{\bigl(\beta^T {\MB
Z}_i\bigr)^T\gamma(v)\Bigr\}K_h\{ \bZ_i-v\}\biggr\}\Biggr],
\end{aligned}
\end{equation}
where $\beta^T \MB Z_i$ and $\beta^T \MB Z_{(j)}$ are defined as
$\beta^T \MB Z$  with $\MB Z$ replaced by $\MB Z_i$ and $\MB
Z_{(j)}$ respectively.
%This version (\ref{LocalLogPL}) is similar to the local partial likelihood in equation (3.5 ) of \cite{FKG:1997}.
It can be shown that the local log partial likelihood in
(\ref{LocalLogPL}) is strictly concave with respect to
$\gamma(\cdot)$, so for a fixed $\beta$, it has a unique maximizer
with respect to $\gamma$. Let $\hat{\gamma}(v)$ be the local partial
likelihood estimate with $\gh_k(v)$ denoting its $k$-th component,
then $\psi^{(k)}(v)$ can be estimated by $\hat{\psi}^{(k)}(v)=k!$
$\gh_k(v)$, for $k=1,\cdots,p$. In principle, one could maximize
(\ref{LocalLogPL}) with respect to both $\beta$ and $\gamma$, and
this corresponds to maximizing the real local log likelihood. But we
choose to maximize
 (\ref{LocalLogPL}) only with respect to $\gamma$  for a fixed
estimated value of $\beta$, and this corresponds to maximizing a
pseudo local log likelihood as the true $\beta$ in
(\ref{LocalLogPL}) is replaced by an estimate. There are two reasons
for our choice. First (\ref{LocalLogPL}) is concave in $\gamma$, but
not necessarily in $\beta$. Second, maximizing with respect to both
parameters is probably not worth the additional computational cost,
as the local likelihood procedures mainly serve as a smoother and
the choice of the smoother is usually not crucial.

To estimate the link function, we use
\begin{displaymath}
\hat{\psi}(v)=\int_0^v \hat{\psi}'(w)dw,
\end{displaymath}
where $\psih'(v)$ is the first component of $\gh(v)$ at the last
iteration step. There are several ways to approximate this integral,
such as the trapezoidal rule or Gaussian quadrature. For
computational simplicity, we apply the trapezoidal rule in the
simulation studies, as suggested in \cite{TH:1987}, and this appears
to be satisfactory.

\subsection{Algorithm and Computational Issues}

The procedure described in the previous subsection requires a
certain choice of $\beta$ in equation (2). This can be done either
independently or iteratively as once an estimate of $\psi$ is
obtained, one can then estimate $\beta$ through the global partial
likelihood.  An iterative algorithm, as shown below, can be
established by alternatingly updating the estimates for $\beta$ and
$\psi$. Such an iteration procedure may improve the link estimate as
a better estimate of $\beta$ will lead to a better estimate of
$\psi$.

{\B\noindent Step 1}. (a) Assign a nonzero initial value to $\beta$, and call it $\hat{\beta}$.

(b) For a given $v$, plug $\hat{\beta}$ into the pseudo log local
partial likelihood and maximize
\begin{displaymath}\label{LLL}
\begin{aligned}[b]
\sum_{j=1}^N K_h \bigl\{\bh^T Z_{(j)}-v\bigr\}&\cdot \biggl[[\bh^T
\mathbf {Z}_{(j)}]^T\gamma(v) \\
&- \log \biggl\{\sum_{i\in
\R_j}\exp\{[\bh^T \mathbf {Z}_i]^T\gamma(v) \}K_h \bigl\{\bh^T Z_i- v
\bigr\}\biggr\}\biggr]
\end{aligned}
\end{displaymath}
with respect to $\gamma(v)$ to get the estimate $\hat{\gamma}(v)$.

(c) Obtain the values of $\hat{\gamma}(v)$, for $v=\bhZ_i, i=1,
\cdots, n$.

(d) Apply the trapezoidal rule to obtain $\{\psih(\bhZ_i):i=1,
\cdots, n\}$.


{\B\noindent Step 2}. Plug $\psih(\cdot)$ into the log (global)
partial likelihood
\begin{displaymath}\label{GLL}
l_G(\beta,\psih)=\sum_{j=1}^N \biggl[\psih(\bZ_{(j)})- \log
\biggl\{\sum_{i\in \R_j}\exp
\bigl\{\psih(\bZ_i)\bigr\}\biggr\}\biggr],
\end{displaymath}
and maximize it with respect to $\beta$ to update the estimate
$\bh$. We use the angle between two estimated $\bh$ at two
consecutive iterations as the convergence criterion.\\


\iffalse%%%%%%%%%%%%%%%%%%%%


If the initial estimate in 1(a) is already $\sqrt{n} $ consistent,
Theorems 1 and 2 imply that no iteration is required for the link
estimate in Step 1(d), which will converge efficiently at the usual
nonparametric rate.  Namely, the link function can be estimated with
the same efficiency as when $\beta $ is known. One possible choice
of initial estimate in Step 1(a) is the estimator in \cite{CLW:1999}
that extends the sliced inverse regression (SIR) approach to
censored data. Specifically, this approach requires a certain design
condition as listed in (2.3) there but has the advantage that it
leads to $\sqrt{n}$-consistent estimator for $\beta$ without the
need to estimate the link function. Another initial estimate which
does not rely on the design conditions (2.3) in \cite{CLW:1999} is
provided in a Ph.D. thesis \cite{WANG:2001} to be elaborated later
in Step 2*.

If a $\sqrt{n}$ consistent estimate of $\beta$ is not readily
available in Step 1(a),  one could add a step 2 below to update it
iteratively until some convergence criterion is met.

{\B\noindent Step 2}. Plug $\psih(\cdot)$ into the log (global)
partial likelihood
\begin{displaymath}\label{GLL}
l_G(\beta,\psih)=\sum_{j=1}^N \biggl[\psih(\bZ_{(j)})- \log
\biggl\{\sum_{i\in \R_j}\exp
\bigl\{\psih(\bZ_i)\bigr\}\biggr\}\biggr],
\end{displaymath}
and maximize it with respect to $\beta$ to update the estimate
$\bh$.

\noindent We use the angle between two estimated $\bh$ at two
consecutive iterations as the convergence criterion in Step 2.



\fi%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 \textbf{Remark 1}.  The Newton-Raphson method is used to find the
estimators in Step 1 and 2. The initial value of $\beta$ can be set
in different ways but cannot be zero as a nonzero value is needed in
step 1 to estimate the link function. However, this restriction does
not exclude the final $\beta$-estimate to be zero or close to zero.
A simple choice is to fit the usual Cox model and use this estimator
in the first step. To accelerate the computation, one can also use
alternative estimates as described below.

\textbf{Remark 2}. It is possible to accelerate the computation by
using a $\sqrt{n}$-consistent initial estimator, as Theorems 1 and 2
below imply that no iteration is required for the link estimate and
that it will converge efficiently at the usual nonparametric rate.
Namely, the link function can be estimated with the same efficiency
as when $\beta $ is known. In practice, we find that one iteration
helps to improve the numerical performance but further iteration is
usually not necessary. There are two choices for a
$\sqrt{n}$-consistent initial estimator, one is the estimator in
\cite{CLW:1999} that extends the sliced inverse regression (SIR)
approach to censored data. Specifically, this approach requires a
certain design condition as listed in (2.3) there but has the
advantage that it leads to a $\sqrt{n}$-consistent estimator for
$\beta$ without the need to estimate the link function. Another
initial estimate which does not rely on the design conditions (2.3)
in \cite{CLW:1999} is provided in a Ph. D. thesis \cite{Wang:2001}.
Specifically, this involves replacing the $\psi$ function in step 2
above by its local version (\ref{Taylor}), which leads to the
cancelation of the term $\psi$ and  results in a version of log
(global) partial likelihood that involves only the derivative
$\gamma(v)$ of $\psi$ but not $\psi$ itself. Thus, Step 2 above is
replaced by

{\B Step 2*}. Maximize the following approximate log (global)
partial likelihood  with respect to $\beta$:
$$\sum_{j=1}^N \biggl[[ \beta^T {Z}_{(j)}]^T{\gh\{\bhZ_{(j)}\}}-
\log \biggl\{\sum_{i\in \R_j}\exp \Bigl([\beta^T {Z}_i]^T
{\gh\{\bhZ_i\}} \Bigr)\biggr\}\biggr].$$

This approximation may result in some efficiency loss, but has
computational advantages over the estimate in Step 2, since we do
not need to estimate $\psi(v)$ and thus can skip Step 1(d). The
resulting estimate for $\beta$ was shown in  \cite{Wang:2001} to be
$\sqrt{n}$-consistent, consequently an ideal choice as the initial
estimate for $\beta$.


\textbf{Remark 3}. In step 1, the local log partial likelihood in
(\ref{LocalLogPL}) is replaced by a pseudo log partial likelihood
with $\beta$ replaced by $\bh$. As this $\bh$ approaches $\beta$,
the link estimate resulting from maximizing the pseudo log partial
likelihood can be expected to approach the true link function at the
usual nonparametric rate. This is because the parametric estimate
$\bh$ converge to its target at the $\sqrt{n}$--rate, which is faster
than the nonparametric link estimate. A rigorous proof is provided
in Theorem 1 and Theorem 2.


\textbf{Remark 4}. For large sample sizes, it is unnecessary to
estimate the link function at each data point. An alternative way is
to estimate the link function at equal-distance grid points, then
use interpolation or smoothing methods to obtain the estimated
value at each data point. Our simulation results show that this
short-cut is computationally economical while retaining similar
accuracy.

\iffalse%%%%%%%%%%%%%%%%%%


The resulting estimate is the same as the one above with Step 1(c)
replaced by:

{\B Step 1}(c$^*$) Repeat Step 1(b) over a grid of m points,
$\{v_1, v_2, \cdots, v_m\}$ and interpolate with piecewise cubic
splines to get the estimate of $\{\gh(\bhZ_i): i=1, \cdots, n\}$.


 5. Another computational concern arises when the dimension of
the covariates is high, as the estimation of the regression
parameters in Step 2 involves solving nonlinear equations without
close forms. One way to cut down this computing is to use in Step
1(a) a good initial estimate, $\bh$, which is already close to the
true $\beta$. This will reduce the computational effort in Step 2
and the number of iterations needed between Step 1 and Step 2. One
good initial estimate is the aforementioned SIR estimator in
\cite{CLW:1999}. When the design conditions in \cite{CLW:1999} is in
doubt, one could use instead an estimate following Steps 1 (or 1*)
above but replacing Step 2 by:

{\B Step 2*}. Maximize $$\sum_{j=1}^N \biggl[[ \beta^T
{Z}_{(j)}]^T{\gh\{\bhZ_{(j)}\}}- \log \biggl\{\sum_{i\in \R_j}\exp
\Bigl([\beta^T {Z}_i]^T {\gh\{\bhZ_i\}} \Bigr)\biggr\}\biggr]$$
with respect to $\beta$.

One can iterate between Step1(a)-1(c) or 1(c)* and Step 2* until convergence is reached.
This iterative algorithm, first proposed in \cite{WANG:2001}, involves only the derivatives of $\psi(v)$ through $\gamma(v)$ in Step 2*, not $\psi(v)$
itself. Thus Step 1(d) is not needed. The log (global) partial likelihood in Step 2* is an approximated value through a Taylor expansion of
the $\psi(\beta^T Z)$. This approximation may result in some
efficiency loss, but has computational advantage over the
pseudo-profile partial likelihood estimates in Step 2 since we don't need to estimate $\psi(v)$. Thus, this
would be an ideal choice as the initial estimate for $\beta$.





\textbf{Remark 5}. Newton-Raphson method is used to find the
estimators in Step 1 and 2.  The estimated covariance matrix of
$\hat{\beta}$ is given by
\begin{displaymath}
\hat{V}(\hat{\beta})=-\Bigl[\frac{\partial^2
l_G(\beta,\psih)}{\partial
\beta\partial\beta^T}\Big|_{\beta=\hat{\beta}}\Bigr]^{-1}.
\end{displaymath}
\medskip


\fi%%%%%%%%%%%%%%%%%%%%%%%

\iffalse %%%%%%%%%%%%%%%%

 For example, in the process of estimating $\beta$, finding
$\beta$ such that $l_G(\beta,\psih)$ is maximized with respect to
$\beta$ is equivalent to finding the solution $\bh$ to the
likelihood equations
%\begin{displaymath}
$\frac{\partial l_G(\beta,\psih)}{\partial \beta}=0.$
%\end{displaymath}

At the current step with $\bh^{[j]}$, update it to $\bh^{[j+1]}$
through
\begin{displaymath}
\hat{\beta}^{[j+1]}=\hat{\beta}^{[j]}-\Bigl[\frac{\partial^2
l_G(\beta,\psih)}{\partial
\beta\partial\beta^T}\Bigr]^{-1}\Big|_{\beta=\hat{\beta}^{[j]}}\frac{\partial
l_G(\beta,\psih)}{\partial \beta}\Big|_{\beta=\hat{\beta}^{[j]}}.
\end{displaymath}

 If the likelihood function evaluated at $\hat{\beta}^{[j+1]}$ is
less than that evaluated at $\hat{\beta}^{[j]}$, then
$\hat{\beta}^{[j+1]}$ is recomputed using half the step size, which
is $\frac{1}{2} \Bigl[\frac{\partial^2 l_G(\beta,\psih)}{\partial
\beta\partial\beta^T}\Bigr]^{-1}\Big|_{\beta=\hat{\beta}^{[j]}}\frac{\partial
l_G(\beta,\psih)}{\partial \beta}\Big|_{\beta=\hat{\beta}^{[j]}}.$
The iterative scheme continues until convergence is obtained.

\fi%%%%%%%%%%%%%%%%%%%%%


%\subsection{Main results}
\subsection{Main Results}

Let $f(\cdot)$ be the probability density of $\bZ$, for a given
$v$, let $P(t\mid v)=P(X\ge t\mid \bZ=v)$, $Y(t)=I\{X\ge t\}$,
$H=diag\{h,\cdots,h^p\}^T$ and ${\MB u}=\{u,\cdots,u^p\}^T$.

\begin{thm} \label{Theorem 1} Under conditions (\it C1)-(\it C5)
in the Appendix, for any $\sqrt n$ consistent estimator $\bh$ of the
true parameter $\bT$, let $\gh(\cdot)$ be the corresponding
estimator for the derivatives $\gT(\cdot)$ of the true link $\psi$
and $\psih(v)=\int_0^v \psih'(w)dw$, where $\psih'(\cdot)$ is the
first component of $\gh(\cdot)$. If
 $h\rightarrow0$, $nh/\log
n\rightarrow \infty$, $nh^4\rightarrow \infty$   then
\begin{displaymath}
\sup_v|\gh(v)-\gT(v)|\cp0,
\end{displaymath}
and
\begin{displaymath}
\sup_z|\psih(\bhz)-\psi(\bTz)|\rightarrow_p 0
\end{displaymath}
%for any $B>0$.
\end{thm}

\begin{thm} \label{Theorem 2} Under the conditions in Theorem 1 and
for bounded  $nh^{2p+3}$,
\begin{displaymath}
\begin{aligned}[b]
(a)  \sqrt{nh} \biggl\{ &H(\gh(v)-\gT(v))-
{{\psi^{(p+1)}(v)}\over{(p+1)!}}A^{-1}bh^{p+1}\biggr\}\\
&\rightarrow_D
N \biggl\{0,{{\sigma^2(v)}\over{f(v)}} A^{-1}DA^{-1}\biggr\}.
\end{aligned}
\end{displaymath}

Furthermore, we have

\begin{displaymath}
\begin{aligned}
(b) \sqrt{nh} \biggl\{ &H(\gh(\bhz)-\gT(\bTz))-
{{\psi^{(p+1)}(\bTz)}\over{(p+1)!}}A^{-1}bh^{p+1}\biggr\}\\
&\rightarrow_D
N \biggl\{0,{{\sigma^2(\bTz)}\over{f(\bTz)}}
A^{-1}DA^{-1}\biggr\},
\end{aligned}
\end{displaymath}
where $A=\int \MB u{\MB u}^T K(u)du-\nu_1\nu_1^T$, $b=\int
u^{p+1}(\MB u-\nu_1)K(u)du$, $D=\int K^2(u)(\MB
u-\nu_1)^{\otimes2}du$, $\nu_1=\int \MB u K(u)du$, and $
\sigma^2(v)=E\{\delta| \bZ=v\}^{-1}$.
\end{thm}

Theorem 1 establishes the uniform consistency of the local partial
likelihood estimator of $\gamma_0$ and Theorem 2 provides the joint
asymptotic normality of the derivative estimators.  The limiting
distribution of $\gh$ is identical to the one in \cite{FGK:1997},
where $\beta$ is assumed to be known. Thus, there is no efficiency
loss as long as $\beta$ can be estimated at the usual
$\sqrt{n}$-rate.

\subsection{Model Checking and Selection}

While an estimated link function is of interest to correctly reflect
the risk associated with a covariate, a parametric link function is
often preferable to a nonparametric one to lend a parsimonious model
with more interpretable results. Thus, a main incentive to estimate
the link function could be  for exploratory model selection to
facilitate  the choice of a proper parametric link function in the
proportional hazards model (1). If so, the $\beta$ estimate in Step
2 only aids in the link estimation and need not be the end product.
Once a suitable link function has been selected, Theorem 2 can be
used for model checking. For instance, to check the identity link
function under the Cox model, one can test $H_0$: $\psi'(v)=1$.
Since the first component of $\gamma(v)$ is $\psi'(v)$, a local
polynomial of order $p=2$ is usually employed to estimate such a
derivative, and the resulting asymptotic distribution of the
corresponding estimate is given below.
\begin{cor}\label{Corollary 1} Under the condition of Theorem 2, and with
 $p=2$ there, we have
%asyp normal dist as in FGK Corollary 1 on page 1669.
\begin{displaymath}
\begin{aligned}[b]
\sqrt{nh^3}\biggl(\psih'(v)-\psi'(v)-\frac{1}{6} \psi^{(3)}(v)h^2
\frac{\int u^4K(u)du}{\int u^2K(u)du}\biggr)\cd N
\biggl(0,\frac{\sigma^2(v)\int u^2K^2(u)du}{f(v)\int u^2
K(u)du}\biggr).
\end{aligned}
\end{displaymath}
\end{cor}

Corollary 1 facilitates the construction of  testing procedures and
asymptotic simultaneous confidence bands for the link function, but
rigorous asymptotic theory requires much further work and is not
available yet. In principle, one could check the appropriateness of
the link function at all data points $v$ that falls in the range of
$\beta^T Z$. Since the true value of $\beta$ is unknown, it is
natural to replace it with an estimate. However, one must bear in
mind the precision of this estimate as well as the low precision of
$\psi'(v)$ for $v$ in the boundary region of $\bh^TZ$. Here boundary
region is defined as within one bandwidth of the data range, where a
smoothing procedure is employed. Since the bandwidth $h$ is usually
of a higher order than $n^{-\frac{1}{2}}$, the anticipated rate of
convergence for $\hat{\beta}$, we recommend to restrict inference on
$\psi'(v)$ for $v$ that is in the interior and at least one
bandwidth $h$ away from either boundary of the range of $\bh$.

Short of  such a rigorous inference procedure for model checking,
pointwise confidence intervals have often been used as a substitute
for exploratory purposes. In the example in Section 4, we illustrate
how to check the appropriateness of the Cox model, i.e. identity
link function, using pointwise confidence intervals developed from
Corollary 1. Readers should bear in mind that this is only an
exploratory data analysis tool rather than a formal inference
procedure.

\iffalse%%%%%%%%%%%%%%%%%%%%%%%%%%%

The selection of a satisfactory parametric link function could be
implemented through a parametric class of link function, such as one
based on the Box-Cox transfromation. Once a specific link function
has been selected, the estimation of the associated $\beta$ and its
asymptotic results will be be based on such a parametric link
function.



Alternatively, if the sample size is large, say a few hundreds
observations, one could first set aside part of the data to estimate
the link function and then use the remaining data to estimate
$\beta$ applying the estimated link function. Because the link was
estimated independently from $\beta$, the usual inference for
partial likelihood could be applied conditioning on the estimated
link function. Simulation results in the next section suggests that
a sample size of 200 or even 50 yields reliable link estimates.

%\iffalse%%%%%%%%%%%%%%%%%%%%%%%%%%%


 As mentioned earlier, the main purpose to estimate the
unknown link function could be for model checking. The $\beta$
estimate in Step 2 only facilitates the link estimation and need not
be the end product. For a more interpretable result, one may prefer
to use a parametric link function, so the estimation of the
associated $\beta$ and its asymptotic results will be delayed until
such a parametric link function has been satisfactorily located.
Alternatively, if the sample size is large, say a few hundreds
observations, one could first set aside part of the data to estimate
the link function and then use the remaining data to estimate
$\beta$ applying the estimated link function. Because the link was
estimated independently from $\beta$, the usual inference for
partial likelihood could be applied conditioning on the estimated
link function. Simulation results in the next section suggests that
a sample size of 200 or even 50 yields reliable link estimates.



\fi%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Simulation Studies}

To see how the algorithm in Section 2 works for the proposed model,
we conducted simulation studies where a quadratic  link function
$\psi(\bZ)=(\bZ)^2$ with $\beta=(1,3)^T$ and a constant baseline,
$\lambda_0=0.005$, were employed. The design for the two-dimensional
covariate, $Z=(Z_1, Z_2)^T$ is: $Z_1 \sim U(-1,1)$ and $Z_2$ is a
truncated $N(0,1)$ with values in $[-1,1]$. Parameters of $\beta$
were chosen in such a way that the simulation generates a reasonable
signal to noise ratio (cf. Fig 1). If we take $\varepsilon$ to
have the standard exponential distribution, $\exp(1)$, the resulting
hazard function will be $\lambda_0\exp\{\psi(\beta^T Z)\}$, and
survival times from this model can be generated as $T
=\exp\{\psi(-\beta^T Z)\} \varepsilon /\lambda_0$.
Different uniform distributions
were utilized to generate three independent censoring times so that
the censoring rates were 0\%, and roughly 25\% and 50\%. The
Epanechnikov kernel was adopted in the link estimation. Two sample
sizes, 200 and 50, were selected to see whether the methods are
flexible for moderate to small samples. For $n=200$ we used 25
equal-distance grid points to estimate the link function to save
computational time as elucidated  in Remark 4 of Section 2.
Piece-wise spline interpolation  was then used to get the link
estimate during each iteration of the algorithm.

Due to the complication from the identifiability problem, the link
function can only be identified up to a constant. Thus, $\psih(v)$
and $\psih(v)+c$, for any constant c, are considered to be
equivalent procedures, and any measures of performance would declare
these two procedures identical. This points to selecting a measure
which measures the variation instead of  the real difference. We
adopt a measure proposed in \cite{GC:1991} which is the standard
deviation
 of the differences between the fitted values $\psih(\bhZ)$ and the
true values $\psi(\bZ)$ at all data points. More specifically, this
measure,  denoted by  $d$, is   the standard deviation of the
difference $\{\psih(\bhZ)-\psi(\bZ): i=1,\cdots,n\}$. We report in
Table 1 the average values for this measure and its standard
deviation based on 100 simulation runs.

Since at each estimating step, $\bh$ was updated and the range of
$\bhZ$ might be different, we used a bandwidth $h^*$, which took a
certain portion of the range of $\bhZ$. For instance,  an $h^*=0.3$
means that the actual bandwidth is 0.3 times the range of the values
of $\bhZ$.
\iffalse %%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Simulation 1}


The first simulation aims at the performance of the estimation
procedure when the true link is an identity link, that is, when
the true model is the Cox model. The link function was set to
$\psi(\bZ)=\bZ$ with $\beta=(6,5)^T$. Different bandwidths $h^*$
ranging from 0.1 to 1 with increments of 0.1 (Note again $h^*=0.1$
means that the actual bandwidth employed was 0.1 times the data
range of $ \bhZ$) were tried for each sample. Here we reported
only $h^*=0.6$ to 1 times of the data range, since the true link
is an identity link and smaller bandwidth did not perform as well
for such a link function. We compared our procedures with two
other procedures in Tables 1 and 2. Method 1 assumes that the
identity link function is known and the regression coefficient
estimate $\bh$ is therefore the partial likelihood estimate in the
Cox model . Method 2 assumes that $\beta$ is known and estimates
the unknown link function only. Method 3 is the new procedure
where both the link function and $\beta$ are unknown and need to
be estimated from the data.

\small

\begin{table}
\caption{\label{T3.1.1} Comparison of different methods for
estimating $\psi(\beta^T Z)=\beta^T Z$, where $\beta=(6,5)^T$ and
$Z=(Z_1, Z_2)^T$ with $Z_1\sim U(-1,1)$ and $Z_2\sim N(0,1)$
(truncated at [-1,1]), n=200. The numbers before and after ``/"
are the means and standard deviations of $d$ based on 100
simulations.} \centering
%\fbox{
\begin{tabular}{lccccccc}
\tiny
\hline\hline
&&\multicolumn{5}{c}{$h^*$}&Optimal\\
Censoring&$^\dag$& 0.6 &0.7 &0.8 &0.9 &1.0 &MSE\\
\hline No
&1&0.232/0.143&0.232/0.143&0.232/0.143&0.232/0.143&0.232/0.143&0.075\\
&2&0.237/0.147&0.229/0.148&0.224/0.148&0.221/0.149&\fbox{0.220/0.149}&0.070\\
&3&0.243/0.151&0.236/0.150&0.232/0.150&0.230/0.150&\fbox{0.228/0.150}&0.075\\

26\%
&1&0.277/0.187&0.277/0.187&0.277/0.187&0.277/0.187&0.277/0.187&0.111\\
&2&0.263/0.177&0.253/0.178&0.248/0.180&0.245/0.181&\fbox{0.244/0.181}&0.093\\
&3&0.290/0.168&0.286/0.166&0.278/0.165&0.275/0.165&\fbox{0.274/0.164}&0.102\\

53\%
&1&0.348/0.230&0.348/0.230&0.348/0.230&0.348/0.230&0.348/0.230&0.174\\
&2&0.336/0.228&0.322/0.230&0.317/0.232&0.314/0.233&\fbox{0.314/0.233}&0.153\\
&3&0.336/0.206&0.321/0.202&0.313/0.203&0.311/0.204&\fbox{0.310/0.204}&0.138\\

\hline \multicolumn{4}{l}{$^\dag$ \scriptsize  Method 1 is under
identity link and unknown $\beta$.}&
\multicolumn{4}{l}{\scriptsize Method 2 is under unknown link and true $\beta$.}\\
\multicolumn{8}{l}{\hspace{2.5mm}\scriptsize Method 3 is under unknown link and unknown $\beta$.}\\
\hline\hline
\end{tabular}%}
\end{table}

\normalsize

Table \ref{T3.1.1} summarizes the results on the link function
estimates. It lists the average (and standard deviation) of the
aforementioned standard deviation, $d$, of the differences between
the fitted values $\psih(\bhZ)$ and the true values $\psi(\bZ)$ at
each data point. The results of the best procedures together with
the corresponding optimal bandwidths are also highlighted with
boxes. The new approach, which requires the estimation of an
additional parameter $\beta$, compares relatively well to the
procedure which assumes a known $\beta$ value. This is consistent
with the theoretical findings in Theorem 2. Intriguingly, the new
link estimator is comparable to, and often better than, the link
function estimate, $\bhz$, generated directly from the Cox model.
This occurs for both $n=200$ and the small sample case, $n=50$,
reported in Table \ref{T3.1.2}. There seems to have no advantage
in knowing the linear link function. This is probably due to the
linear structure of the link function and the restriction of
$\|\beta\|$=1. The latter means that there is one less parameter to
estimate under our model than under the Cox model, giving some
advantage to our procedure. When the link function is nonlinear as
in Simulation study 2 below, knowledge of the truth link function
definitely has an advantage.

Figure \ref{F3.1} shows the results of four simulated samples for
three estimated curves from the Cox model with known link, the
procedure with known $\beta$,
and the proposed procedure with both the link function and $\beta$ unknown. %Note that $\psi-\log(\varepsilon)=\log\Lambda_0(T)$, where
%$\Lambda_0(t)=\int_0^t \lambda_0(u)du$ is the cumulative baseline hazard function. We refer $\psi-\log(\varepsilon)$ as the simulated data which was plotted at the $Y$ axis. $\bZ$ was used at the $X$ axis.
It's clear that all three curves are close to each other, and the
estimators with known $\beta$ and the new estimators are almost
identical, which again confirms that the efficiency loss is small
when the coefficients for covariates are unknown and estimated via
the profile method.




\begin{table}
\caption{\label{T3.1.2}Comparison of different methods for
estimating $\psi(\beta^T Z)=\beta^T Z$, where $\beta=(6,5)^T$ and
$Z=(Z_1, Z_2)^T$ with $Z_1\sim U(-1,1)$ and $Z_2\sim N(0,1)$
(truncated at [-1,1]), n=50. The reported quantities are the
smallest MSEs under the optimal bandwidths.}\centering
%\fbox{
\begin{tabular}{lcccccc}
\hline\hline
&\multicolumn{2}{c}{No Censoring}&\multicolumn{2}{c}{24\% Censoring}&\multicolumn{2}{c}{52\% Censoring}\\
Method& d&Optimal $h^*$ & d&Optimal $h^*$ & d&Optimal $h^*$ \\
\hline
Identity link and unknown $\beta$&0.638&-&0.675&-&1.615&- \\
Unknown link and true $\beta$&0.551&1.0&0.621&1.0&1.323&1.0\\
Unknown link and unknown $\beta$&0.553&1.0&0.558&1.0&1.392&1.0\\
\hline\hline
\end{tabular}
\end{table}

We note here that we have also conducted simulations for the
regression parameters $\beta$, and our regression estimator is
also comparable to the Cox estimator.  These simulations are not
included here, but the next simulation study contains a comparison
of the various regression estimates.

\subsection{Simulation 2}

Next we performed a simulation similar to the previous one, but
with a quadratic link function, $\psi(\bZ)=(\bZ)^2$, where
$\beta=(1,3)^T$.


\fi%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Various bandwidths were explored, but we report only the results for
bandwidth $h^*$ varying from 0.1 to 0.4 (with 0.1 increment) times
the data range of $\bZ$ at each iteration stage. Results for other
bandwidths were inferior and are not reported here.

\begin{table}
\caption{\label{T3.2.1} Comparison of different methods for
estimating $\psi(\beta^T Z)=(\beta^T Z)^2$, where $\beta=(1,3)^T$
and $Z=(Z_1, Z_2)^T$ with $Z_1\sim U(-1,1)$ and $Z_2\sim N(0,1)$
(truncated at [-1,1]), n=200. The numbers before and after ``/"
are the means and standard deviations of $d$ based on 100
simulations. }\centering
\scriptsize
\begin{tabular}{lp{.1cm}ccccp{.7cm}}
\hline\hline
&&\multicolumn{4}{c}{$h^*$}&Optimal\\
Censoring&$^\dag$& 0.1 &0.2 &0.3&0.4&MSE\\
\hline No
&1&3.152/0.181&3.152/0.181&3.152/0.181&3.152/0.181&9.970\\
&4&0.191/0.114&0.191/0.114&0.191/0.114&0.191/0.114&0.049\\
&2&1.234/2.554&\fbox{0.439/0.142}&0.843/0.150&1.417/0.151&\fbox{0.213}\\
&3&1.210/3.522&\fbox{0.490/0.144}&0.907/0.168&1.476/0.168&\fbox{0.261}\\
25\%
&1&3.151/0.181&3.151/0.181&3.151/0.181&3.151/0.181&9.961\\
&4&0.201/0.112&0.201/0.112&0.201/0.112&0.201/0.112&0.053\\
&2&1.256/2.548&\fbox{0.425/0.157}&0.684/0.172&1.197/0.178&\fbox{0.206}\\
&3&0.981/1.991&\fbox{0.468/0.156}&0.746/0.185&1.256/0.194&\fbox{0.244}\\
50\%
&1&3.149/0.181&3.149/0.181&3.149/0.181&3.149/0.181&9.947\\
&4&0.210/0.117&0.210/0.117&0.210/0.117&0.210/0.117&0.056\\
&2&1.361/2.529&\fbox{0.525/0.215}&0.535/0.161&0.729/0.210&\fbox{0.322}\\
&3&1.236/2.138&\fbox{0.536/0.201}&0.571/0.169&0.808/0.226&\fbox{0.327}\\
%\vspace{2mm}
\hline \multicolumn{4}{l}{$^\dag$ \scriptsize Method 1 is under
identity link}&
\multicolumn{3}{l}{\scriptsize Method 2 is under unknown link}\\
\multicolumn{3}{l}{\scriptsize \hspace{5pc} and unknown $\beta$.} &
\multicolumn{3}{l}{\scriptsize \hspace{11pc} and true $\beta$.}\\
\multicolumn{4}{l}{\hspace{1.3mm} \scriptsize Method 3 is under
unknown link}&
\multicolumn{3}{l}{\scriptsize Method 4 is under quadratic link}\\
\multicolumn{3}{l}{\scriptsize \hspace{5pc} and unknown $\beta$.} &
\multicolumn{3}{l}{\scriptsize \hspace{11pc} and unknown $\beta$.}\\
\hline\hline
\end{tabular}
\end{table}
\normalsize


Four procedures were compared and the results for $n=200$ are shown
in Table \ref{T3.2.1}.   Method 1 assumes that the link is identity
(which is incorrect here) and the regression coefficient estimate
$\bh$ is therefore the Cox estimate based on the partial likelihood
estimate. The aim is to see the effect of erroneously assuming the
conventional Cox proportional hazards model. Method 2 assumes that
$\beta$ is known and estimates the unknown link function as in
\cite{FGK:1997}. Method 3 is the new procedure where both the link
function and regression coefficient $\beta$ are estimated.  Method 4
assumes that the true quadratic link function is known and the
regression coefficient estimate, $\bh$, is the partial likelihood
estimate. The comparisons for the distance $d$ are reported in Table
\ref{T3.2.1}. The results of the best procedures together with the
corresponding optimal bandwidths are highlighted with boxes. It is
not surprising that the best results came from the procedures with
true quadratic link function and unknown $\beta$. Our estimators are
 close to those from method 2 (\cite{FGK:1997}) with known $\beta$,
while the estimators based on the identity link model have much
larger $d$.

To demonstrate the effect of an erroneous link function on
regression estimates, we report in Table \ref{T3.2.2} the results of
the various estimates for $\beta$. Since there is no regression
parameter estimate for method 2, only three procedures are compared
in Table \ref{T3.2.2}. There are several ways to compare the
regression estimates, one way is to set the first component of
$\beta$ to the true value, then compare the difference between
$\bh$ and $\beta$ based on the second component. Table \ref{T3.2.2}
shows the results of the difference between the true $\beta_2$ and
the estimate $\hat{\beta}_2$ for various procedures. The best
procedures for the profile estimator in method 3 are shown in boxes
under the optimal bandwidths.  Another way to compare the different
estimators is to calculate the angles between these estimators and
the true parameter.  To save space only the MSEs based on the
optimal bandwidths are listed in the last column of Table
\ref{T3.2.2} for the angle measure (degrees). Based on both optimal
MSE measures reported in the last two columns of Table \ref{T3.2.2},
the differences between the new profile estimators and the true
parameters are way smaller than those from the identity link model,
and reasonably close to those under the true link.

\begin{table}
\caption{\label{T3.2.2}Differences between the estimated $\bh_2$ and
the true $\beta_2$ where $\bh_1$ is set equal to the true $\beta_1$,
$\psi(\beta^T Z)=(\beta^T Z)^2$, where $\beta=(1,3)^T$ and $Z=(Z_1,
Z_2)^T$ with $Z_1\sim U(-1,1)$ and $Z_2\sim N(0,1)$ (truncated at
[-1,1]), n=200. The numbers before and after ``/" are the biases and
standard deviations of the $\bh_2$ based on 100 simulations.}
\tiny
\centering
\begin{tabular}{lccccccc} \hline\hline
&&\multicolumn{4}{c}{h$^*$}&Optimal&Optimal\\
Censoring&$^\dag$& 0.1 &0.2 &0.3 &0.4 &MSE$^\ddag$ &MSE$^{\S}$\\
\hline No
&1&5.344/50.847&5.344/50.847&5.344/50.847&5.344/50.847&2613.949&11742.108\\
&4&-0.015/0.142&-0.015/0.142&-0.015/0.142&-0.015/0.142&0.021&0.703\\
&3&0.023/0.276&\fbox{-0.051/0.154}&-0.096/0.158&-0.215/0.190&0.026&\fbox{0.968}\\
25\%
&1&-1.229/24.742&-1.229/24.742&-1.229/24.742&-1.229/24.742&613.684&11320.063\\
&4&-0.013/0.154&-0.013/0.154&-0.013/0.154&-0.013/0.154&0.024&0.801\\
&3&0.038/0.307&\fbox{-0.055/0.170}&-0.089/0.171&-0.162/0.195&0.032&\fbox{1.144}\\
50\%
&1&-1.472/6.486&-1.472/6.486&-1.472/6.486&-1.472/6.486&44.237&11266.339\\
&4&-0.006/0.170&-0.006/0.170&-0.006/0.170&-0.006/0.170&0.029&0.968\\
&3&0.055/0.371&\fbox{-0.050/0.188}&-0.075/0.181&-0.110/0.190&0.038&\fbox{1.362}\\
\hline 
\multicolumn{4}{l}{$^\dag$ \tiny Method 1 is under
identity link and unknown $\beta$.}&
\multicolumn{3}{l}{\tiny Method 4 is under true link and unknown $\beta$.}\\
\multicolumn{8}{l}{\hspace{1.3mm} \tiny Method 3 is under unknown link and unknown $\beta$.}\\
\multicolumn{4}{l}{\tiny $^\ddag$ is the optimal MSE for
$\bh_2$ when $\bh_1$ is set to be 1.}&
\multicolumn{4}{l}{\tiny $^{\S}$ is the optimal MSE for
the angles between $\bh$ and the true $\beta$.}\\
\hline\hline
\end{tabular}
\normalsize
\end{table}

\underline{} We can see that using the wrong link will lead to huge
bias and MSE under all censoring patterns. The average angles
between the $\beta$ estimates assuming identity link and the true
parameters are around 90$^\circ$, which suggests that the $\beta$
estimates with identity link are perpendicular to the true parameter
space, indicating a total inability to estimate the regression
parameter. This is in addition to the problem that the link function
itself has been misspecified. Both underscore the importance to
check the shape of the link function at the beginning of data
analysis.

Four typical simulated curves are shown in Fig \ref{F3.2}. The
procedure (method 4) with known quadratic link function and unknown
$\beta$ performed the best. The procedure (method 2) with known
$\beta$ and our procedure captured the shape of the true curve well,
but the procedure (method 1)  based on the Cox model failed to
capture the shape of the link function. Results for $n=50$
summarized in Table \ref{T3.2.3} are consistent with the findings
for $n=200$ in Table \ref{T3.2.1}.

%[bb=0 0 700 700, scale = 1.6, width=2.5in]

\vspace{4ex}
\begin{figure}
\centering
\begin{tabular}{cc}
\begin{minipage}{52mm}
%\includegraphics[scale=2.8, bb= 0 0 700 700, width=2.6in, height=2.6in, angle=0]
\includegraphics[width=2.25in, height=2.3in, angle=0]
	{wang_simu2_n200_g25_1.pdf}  \vspace{-5pc}
%\includegraphics[scale=2.8, bb= 0 0 700 700, width=2.6in, height=2.6in, angle=0]
\includegraphics[width=2.25in, height=2.3in, angle=0]
	{wang_simu2_n200_g25_3.pdf}\vspace{7pc}
\end{minipage}\hspace{3pc}
\begin{minipage}{52mm}
\vspace{-2pc}\includegraphics[width=2.25in, height=2.3in, angle=0]
	{wang_simu2_n200_g25_2.pdf}\vspace{ 0pc} 
\includegraphics[width=2.25in, height=2.3in, angle=0]
	{wang_simu2_n200_g25_4.pdf}
\end{minipage}
%\multicolumn{2}{c}
\end{tabular}
\begin{tabular}{c}
\hspace{-19pc}\vspace{-26pc}$\bZ$
\end{tabular}
\caption{\label{F3.2} Four typical simulated sets of data and
estimated curves from the model $\psi(\beta^T Z)=(\beta^T Z)^2$ with
$\beta=(1,3)^T, Z=(Z_1,Z_2)^T$ where $Z_1\sim U(-1,1)$ and $Z_2\sim
N(0,1)$ (truncated at [-1,1]), $25\%$ censoring and $h^*=0.2$. The
simulated data is indicated as black dots. The true curve
$\psi(\cdot)$ is indicated as a red solid line, and the estimated
curves $\hat{\psi}(\cdot)$ as a magenta dash line with identity link
and unknown $\beta$, as a black dash-dotted line with quadratic link
and unknown $\beta$, as a green dotted line with unknown link and
true $\beta$, and as a blue dash-dotted line with both unknown link
and unknown $\beta$.}
\end{figure}
\normalsize


\begin{table}
\caption{\label{T3.2.3} Comparison of different methods for
estimating $\psi(\beta^T Z)=(\beta^T Z)^2$, where $\beta=(1,3)^T$
and $Z=(Z_1, Z_2)^T$ with $Z_1\sim U(-1,1)$ and $Z_2\sim N(0,1)$
(truncated at [-1,1]), n=50. The reported quantities are the
smallest MSEs under optimal bandwidths.} \centering
%\fbox{
\scriptsize
\begin{tabular}{lcccccc}
\hline\hline
&\multicolumn{2}{c}{No Censoring}&\multicolumn{2}{c}{25\% Censoring}&\multicolumn{2}{c}{50\% Censoring}\\
Method& d&Optimal $h^*$ & d&Optimal $h^*$ & d&Optimal $h^*$ \\
\hline
Identity link and unknown $\beta$&9.199&-  &9.168&-&9.102&-    \\
Quadratic link and unknown $\beta$&0.184&-  &0.203&-&0.295&-    \\
Unknown link and true $\beta$ &0.808&0.3&0.833&0.3&1.359&0.4\\
Unknown link and unknown $\beta$&0.945&0.4&1.286&0.4&1.423&0.4\\
\hline\hline
\end{tabular}
\end{table}
\normalsize

\section{Data Analysis}

In this section, we illustrate the proposed model and estimation
algorithm in Section 2 through the Worcester Heart Attack Study
(WHAS) data. One of the goals of this study is to identify factors
associated with the survival rates following hospital admission for
acute myocardial infarction. The main data set has more than 11000
admissions, but we used only a random sample of 500 patients as
listed in \cite{HL:2008}. This data set is chosen because the
proportionality assumption has been carefully examined and is
reasonably satisfied. Our goal here is to check the adequacy of the
identity link function in the Cox proportional hazards regression
model.

There were more than 10 covariates in the data set. After detailed
model selection procedure, Hosmer et al. \cite{HL:2008} included 6
variables age (AGE), initial heart rate (HR), initial diastolic
blood pressure (DIASBP), body mass index (BMI), gender (GENDER),
congestive heart complications (CHF), and the interaction between
age and gender (AGEGENDER) in their model. After examining the
linearity assumption using fractional polynomials, they decided to
apply a two-term fractional polynomial model to the variable BMI. We
thus begin with the univariate covariate BMI.


We tried different bandwidths and found similar patterns of the
estimated link functions. In Fig \ref{F4.1} we report two of the
results, which exhibit reasonable level of smoothness. The estimated
link function in Fig \ref{F4.1} suggests clear nonlinearity. We
then constructed a 95\% point-wise approximated confidence interval
of $\gamma(v)$ (Fig \ref{F4.2}) to see whether it would cover the
constant function 1. The results suggest that the estimated link
functions have some curvature and further investigation is needed.

\begin{figure}[h!]
\centering \hspace{1pc}WHAS data\\
\makebox{\includegraphics[height=3.0in,width=3.5in]{wang_whas_BMI_psi_3in1_v2.pdf}}
\caption{\label{F4.1} The estimated link function, under two
bandwidths, for the WHAS data with BMI as a covariate.}
\end{figure}


\begin{figure}[h!]
\centering \hspace{1pc}WHAS data\\
\makebox{\includegraphics[height=3.0in,width=3.5in]{wang_whas_BMI_gamma_h6r.pdf}}
\caption{\label{F4.2} The estimated confidence interval of $\gamma$
at bandwidth $h^*=1/6$ for the WHAS data based on BMI.}
\end{figure}


Next we applied the proposed procedure to the multivariate model
with all 7 covariates.  We tried different bandwidths  ranging from
1/10, 1/8, 1/7, 1/6, 1/5, 1/4 ,1/3, to 1/2 of the single index
range, and plotted the results of three bandwidths in Fig
\ref{F4.3}. The estimated link function for $h^*=1/4$ appears
oversmoothed but all three estimates exhibit two bumps. The 95\%
point-wise approximated confidence intervals for $\gamma(v)$ as
shown in Fig \ref{F4.4} also reveal curvature away from the
constant ($=1$) horizontal line. Although it is arguable that the
Cox model could be rejected at a low level, the suitability of an
identity link function seems questionable.


\iffalse%%%%%%%%%%%%%%%%%%%%%%%%

The directions of the estimate $\bh$ based on the new procedure and
of the estimate $\hat{\beta}_{Cox}$ based on the Cox model with an
identity link were then compared. The angles between $\bh$ and
$\beta_{cox}$ range from 0.3$^{\circ}$ to $5.2^{\circ}$ as the
bandwidth decreases from the single index range to 1/10 of the
single index range.  We further adjusted the coefficient of the
age-estimate to be the same for both estimates, so that we could
compare the actual difference between the new estimate $\bh$ and the
Cox estimate $\hat{\beta}_{Cox}$ in size of the remaining four
coefficients. Since $\bh_{Cox}$ is $(0.0623, 2.2492, 0.0112,
-0.0108, -0.0439, 0.7757, -0.0320)^T$, the scale-adjusted $\bh$ is
$(0.0623, 2.2884, 0.0115,-0.0141,-0.0383, 0.8239, -0.0332)^T$ when
$h^*=1/6$. All seven estimates lie in the 95\% confidence interval
of the Cox estimates after the scale change.


\fi%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{figure}[h!]
\centering \hspace{1pc}WHAS data\\
\makebox{\includegraphics[height=3.0in,width=3.5in]{wang_whas_all_psi_4in1_v4.pdf}}
\caption{\label{F4.3} The estimated link function for the WHAS data
under three different bandwidths.}
\end{figure}

\begin{figure}[h!]
\centering \hspace{1pc}WHAS data\\
\makebox{\includegraphics[height=3.0in,width=3.5in]{wang_whas_all_gamma_h6r.pdf}}
\caption{\label{F4.4}The estimated confidence interval of $\gamma$
for the WHAS data at  bandwidth $h^*=1/6$.}
\end{figure}


\section{Conclusion and Future Research}

The proposed estimating procedures for the extended proportional
hazards regression model with unknown link function and
multi-dimensional covariates seem to be reliable for moderate to
large sample sizes.  Once the link function and the parameters of
the index have been established, one can proceed to estimate the
unknown baseline hazard function in model (\ref{PHM}) using a Breslow-type
estimate (\cite{Breslow:1974}).

The cost of a misspecified link function has been demonstrated
through the simulation studies in Section 3. As a consequence, the
risk of an individual may be misinterpreted. It is thus important to
at least estimate the link function in the initial model fitting
stage as a model checking tool or guidance to a suitable class of
parametric link functions. A rigorous  test of parametric link
function will be a worthwhile  future project, as is the asymptotic
theory for simultaneous inference of the link function and
regression parameters.

The choice of automatic smoothing parameters, the bandwidth $h$ in
this case, is a challenging problem for proportional hazards model
when a likelihood based smoother, such as the local partial
likelihood estimate, is employed in the link estimate. This is
because the components of the partial likelihood are dependent,
hence the usual automatic bandwidth selection methods for linear
models are not applicable here. The usual least square
cross-validation procedure in nonparametric regressions also cannot
be easily adapted to hazard based models such as the proportional
hazards model.  An alternative criterion, less computational
intensive than cross-validation methods, was proposed in
\cite{TH:1987}, based on a variation of the Akaike's information
criterion for span selection, when the nearest neighborhood method
was used for smoothing. However, the interpretation of AIC is not
clear here since partial likelihood involves dependent components.
The authors also acknowledged that the asymptotic correctness of the
AIC criterion has not been established.  Thus, automatic bandwidth
choice remains an open question when the link function is being
estimated. Meanwhile, we recommend to try several bandwidths and
choose one that yields a moderately smooth link function as we did
for the WHAS data. This subjective choice based on the visual degree
of smoothness is commonly adopted as an ad hoc tool.

While this paper deals with time-independent covariates, it would be
desirable to extend model (1) to time-dependent covariates as well.
One complication is that the entire history of the covariate process
would be required or some kind of imputation needs to be performed
to get even an initial estimate of $\beta$. Preliminary results were
reported in \cite{Wang:2004} by imputing the covariate process
through a functional principal components approach, and then
proceeding with the estimation of  the survival components at the
second stage. Such a two-stage procedure is prone to bias as is well
known in the joint modelling literature. Further investigations to
correct the bias would be desirable, and joint modelling the
longitudinal and survival process offers some hope if one can
resolve the additional complication of an unknown link function.
This is yet another worthwhile project to pursue in the future.

\section*{Acknowledgement}

The work of Jane-Ling Wang and Wei Wang was completed with the
partial support of National Science Foundation grants DMS04-06430,
National Institutes of Health grant R01-DK-45939, and National
Institute of Allergy and Infectious Diseases grant AI24643.
 Qihua Wang's research was supported by the  National
Science Fund for Distinguished Young Scholars in China (10725106),
the National Natural Science Foundation of China (10671198), the
National Science Fund for Creative Research Groups in China and the
Research Grants Council of the Hong Kong (HKU 7050/06P). The authors
are grateful for the invaluable suggestions of two reviewers and the
editor.

\clearpage
\section*{Appendix}

Let $f(\cdot)$ be the probability density of $\bZ$, for a given
$v$, let $P(t\mid v)=P(X\ge t\mid \bZ=v)$, $\Lambda(t,v)=\int_0^t
P(u\mid v)\lambda_0(u)du$, $Y(t)=I\{X\ge t\}$, $Y_i(t)=I\{X_i\ge
t\}$, $H=diag\{h,\cdots,h^p\}^T$ and ${\MB u}=\{u,\cdots,u^p\}$.

We begin with some regularity conditions needed for the results.

$(\it C1)\;K\ge0$ is a bounded density with compact support, and it
has bounded first and second derivative.
\par
$(\it C2)\;\psi(\cdot)$ has a continuous $(p+1)$th derivative
around $v$.
\par
$(\it C3)\;$The density $f(\cdot)$ of $\bZ$ is continuous at point
$v$ and $\inf_{w} f(w)>0$.
\par
$(\it C4)\;$ The conditional probability $P(t\mid \cdot)$ is
equicontinuous at $v$.
\par
$(\it C5)\;$ $\int_0^\tau\lambda_0(u)du<\infty$.

Denote
\par$s_0(\beta,\psi,u)=E\bigl[Y(u)\exp\{\psi(\bZ)\}\bigr],$
\par$s_1(\beta,\psi,u)=E\bigl[Y(u)\exp\{\psi(\bZ)\}\psi'(\bZ)
Z\bigr],$
\par$s_2(\beta,\psi,u)=E\bigl[Y(u)\exp\{\psi(\bZ)\}\bigl\{\psi''(\bZ)+[\psi'(\bZ)]^2\bigr\}ZZ^T\bigr],$
\par$s^*_2(\beta,\psi,u)=E\bigl[Y(u)\exp\{\psi(\bZ)\}[{\psi'}(\bZ)]^2ZZ^T\bigr].$

$(\it C6)\;$ The functions $ s_r, r=0,1, $and $2,$ and $s^*_2$ are
bounded and $s_0$ is bounded away from 0 on
$\mathcal{B}\times[0,\tau]$; the family of functions
$s_r(\cdot,\psi,u)$ and $s^*_2(\cdot,\psi,u)$ is an equicontinuous
family at $\beta_0$.

The following lemma is used repeatedly in the later proofs. The
proof will be omitted and can be found in \cite{Wang:2001}.

\begin{lem} Let $c_n(\bT,t)=n^{-1} \sum_{i=1}^n Y_i(t)
g(\bTZ_i)K_h(\bTZ_i-v)$ and $c(t)=f(v)g(v)P(t\mid v)\int K(u)du$,
Under conditions ({\it C1}) and ({\it C4}), if $g(\cdot)$ is
continuous at the point $v$, then
\begin{displaymath}
\sup_{0\le t \le \tau} | c_n(\bT,t)-c(t)| \rightarrow_P 0 ,
\end{displaymath}
provided that $h \rightarrow 0, nh/\log n \rightarrow \infty, 0 <
\tau\le +\infty$.

If furthermore, $\bh$ is a $\sqrt{n}$-consistent estimate of $\bT$
and $nh^4\rightarrow\infty$, then
\begin{displaymath}
\sup_{0\le t \le \tau} | c_n(\bh,t)-c(t)| \rightarrow_P 0.
\end{displaymath}
\end{lem}

\begin{center}
Proof of Theorem 1
\end{center}
%\begin{proof}
\noindent\textit{Proof. }For notation simplicity, we will use $\gamma$ to represent
$\gamma(v)$ and $\gT$ for $\gT(v)$.

The log local partial likelihood function at point $v$ is given as
\begin{displaymath}\label{llpl}
\begin{aligned}[e]
l\{\beta,\gamma,v\}=&\frac{1}{n}\sum_{j=1}^N K_h(\beta^T
Z_{(j)}-v)\nonumber\\
&\biggl[\bigl[(\beta^T \MB {Z}_{(j)})^T\gamma\bigr]-\log
\biggl\{\sum_{i \in \R_j}\exp\{(\beta^T
\MB{Z}_i)^T{\gamma}\}K_h(\bZ_i- v) \bigr\} \biggr].
\end{aligned}
\end{displaymath}

Using counting process notation $N(t)=I\{X\le t, \delta=1\}$ and
$N_i(t)=I\{X_i\le t,\delta_i=1\}$, under the independent
censoring,
\begin{displaymath}
M_i(t)=N_i(t)-\int_0^t Y_i(u) \exp\{\psi(\bTZ_i)\}\lambda_0(u)du,
\end{displaymath}
is a martingale with respect to the filtration
$\mathcal{F}_t=\sigma\{N(u),I_{\{X\le u,\delta=0\}}: 0\le u\le
t\}$.

The empirical counterpart of $l\{\beta,\gamma,v\}$ up to time $t$ is
\begin{displaymath}
\begin{aligned}[e]
l_n(\beta,\gamma,t,v)=&\int_0^t \frac{1}{n}\sum_{i=1}^n K_h(
\bZ_i- v)\biggl[\bigl[(\beta^T
\MB{Z}_i)^T {\gamma}\bigr]\\
&-\log \bigl\{\sum_i^n Y_i(u)\exp\{(\beta^T \MB{Z}_i)^T
{\gamma}\}K_h(\bZ_i- v) \bigr\} \biggr]dN_i(u).
\end{aligned}
\end{displaymath}

Denote
%\hspace{1mm} $S^*_{h,0}(\beta,u,v)=\frac{1}{n}\sum_{i=1}^n K_h(\bZ_i-v)Y_i(u)\exp\{\psi(\bZ_i)\}$,
%\par \hspace{8mm}$S^*_{h,1}(\beta,u,v)=\frac{1}{n}\sum_{i=1}^n K_h(\bZ_i-v)Y_i(u)\exp\{\psi(\bZ_i)\}(\beta^T \MB Z_i)$,
%\par and
$S_{h,0}(\beta,\gamma,u,v)=\frac{1}{n}\sum_{i=1}^n
K_h(\bZ_i-v)Y_i(u)\exp\{(\beta^T \MB Z_i)^T\gamma\}$.

Let $\bh$ be a $\sqrt{n}-$consistent estimate of the true parameter
$\bT$,  and $\gh$ be the corresponding estimate of the true $\gT$,
we can write

%\begin{displaymath}
%\begin{aligned}[e]
\begin{align*}
&l_n(\bh,\gamma,\tau,v)-l_n(\bT,\gT,\tau,v)\\
=&\int_0^\tau \frac{1}{n}\sum_{i=1}^n
K_h(\bTZ_i-v)\nonumber\\
&\hspace{.5in}{}\times\biggl[\bigl[(\bT^T \MB{Z}_i)^T \gamma-(\bT^T
\MB{Z}_i)^T \gT\bigr] - \log\frac{ S_{h,0}(\bT,\gamma,u,v)}{
S_{h,0}(\bT,\gT,u,v)} \biggr]dM_i(u)\\\linebreak
+&\int_0^\tau \frac{1}{n}\sum_{i=1}^n
K_h(\bhZ_i-v)\biggl[\bigl[(\bh^T \MB{Z}_i)^T \gamma-(\bT^T
\MB{Z}_i)^T \gamma\bigr]-\log \frac{S_{h,0}(\bh,\gamma,u,v)}
{S_{h,0}(\bT,\gamma,u,v)} \biggr]dM_i(u)\\
+&\int_0^\tau \frac{1}{n}\sum_{i=1}^n\Bigl(
K_h(\bhZ_i-v)-K_h(\bTZ_i-v)\Bigr)\nonumber\\
&\hspace{.5in}{}\times\biggl[(\bT^T \MB{Z}_i)^T \gamma-
\log S_{h,0}(\bT,\gamma,u,v) \biggr]dM_i(u)\\
+&\int_0^\tau \frac{1}{n}\sum_{i=1}^n
K_h(\bTZ_i-v)\biggl[\bigl[(\bT^T \MB{Z}_i)^T \gamma-(\bT^T
\MB{Z}_i)^T \gT\bigr]-\log\frac{ S_{h,0}(\bT,\gamma,u,v)}{
S_{h,0}(\bT,\gT,u,v)}\biggr]\\
&\qquad \times Y_i(u)\exp\{\psi(\bTZ_i)\}\lambda_0(u)du\\
+&\int_0^\tau \frac{1}{n}\sum_{i=1}^n
K_h(\bhZ_i-v)\biggl[\bigl[(\bh^T \MB{Z}_i)^T \gamma-(\bT^T
\MB{Z}_i)^T \gamma\bigr] - \log \frac{S_{h,0}(\bh,\gamma,u,v)}{
S_{h,0}(\bT,\gamma,u,v)}\biggr]\\
&\qquad \times
Y_i(u)\exp\{\psi(\bTZ_i)\}\lambda_0(u)du\\
+&\int_0^\tau \frac{1}{n}\sum_{i=1}^n\Bigl(
K_h(\bhZ_i-v)-K_h(\bTZ_i-v)\Bigr)
 \biggl[(\bT^T \MB{Z}_i)^T \gamma-
\log S_{h,0}(\bT,\gamma,u,v) \biggr]\\
\equiv&X_n(\bT,\gamma,\tau,v)+I+II+A_n(\bT,\gamma,\tau,v)+III+IV.
\end{align*}
%\end{aligned}
%\end{displaymath}

Under the regularity conditions and from Lemma 1, it can be shown
that

(1)\hspace{0.2cm} $X_n(\bT,\gamma,\tau,v)$ is a locally square
integrable martingale with the predictable variation process
\begin{displaymath}
\begin{aligned}[b]
&\langle X_n(\bT,\gamma,\tau,v), X_n(\bT,\gamma,\tau,v)\rangle\\
=&\int_0^\tau \frac{1}{n^2}\sum_{i=1}^n K_h^2(\bTZ_i-v)
\biggl[(\bT^T \MB{Z}_i)^T (\gamma-\gT)-\log
\frac{S_{h,0}(\bT,\gamma,u,v)} {S_{h,0}(\bT,\gT,u,v)} \biggr]^2\\
&\qquad \qquad \qquad \qquad \times Y_i(u)\exp\{\psi(\bTZ_i)\}\lambda_0(u)du\\
=&O_p(\frac{1}{nh}).
\end{aligned}
\end{displaymath}

(2) \hspace{0.2cm} $A_n(\bT,\gamma,\tau,v)\cp$ \begin{displaymath}
\begin{aligned}[e]
  &f(v)\exp\{\psi(v)\}\Lambda(\tau,v)\\
&\times\biggl[(\int \MB u K(u)du) ^T H(\gamma-\gT)-
\log\{\int \exp\{\MB u^T H(\gamma-\gT)\}K(u)du\}\biggr]+o_p(1)\\
&\equiv A(\bT,\gamma,\tau,v)+o_p(1),
\end{aligned}
\end{displaymath}


(3) \hspace{0.2cm} $I=O_p(\frac{1}{nh})$,
$II=O_p(\frac{1}{nh^2})$, $III=O_p(\frac{1}{\sqrt{nh^2}})$, and $
IV=O_p(\frac{1}{\sqrt{nh^4}}).$

This means $X_n(\bT,\gamma,\tau,v)$, I, II, III and IV converge to
zero at a faster rate than $A_n(\bT,\gamma,\tau,v)$. By Lemma
8.2.1(2) in \cite{FH:1991},
$l_n(\bh,\gamma,\tau,v)-l_n(\bT,\gT,\tau,v)$ has the same limiting
distribution as $A_n(\bT,\gamma,\tau,v)$. Thus, we have
\begin{equation}\label{Thm1}
l_n(\bh,\gamma,\tau,v)-l_n(\bT,\gT,\tau,v)\cp
A(\bT,\gamma,\tau,v).
\end{equation}

It is obvious that $A(\bT,\gamma,\tau,v)$ is strictly concave,
with a maximum at $\gamma=\gT$. Hence, the right-hand side of
(\ref{Thm1}) is maximized at $\gamma=\gT$.  The left-hand side of
(\ref{Thm1}) is maximized at $\gamma=\gh$, since $\gh$ maximizes
$l_n(\bh,\gamma,\tau,v)$. Therefore, $\sup_v
\big|\gh(v)-\gT(v)\big|\rightarrow_p0.$

By Dominated Convergence Theorem, we have
\begin{equation}\label{Con:gamma}
\sup_v\big|\psih(v)-\psi(v)\big|\rightarrow_p0.
\end{equation}
This implies
\begin{displaymath}
\begin{aligned}[b]
&\sup_z \big|\psih(\bhz)-\psi(\bTz)\big| \le\sup_z \big|
\psih(\bhz)-\psi(\bhz)\big|+\sup_z
\big|\psi(\bhz)-\psi(\bTz)\big|\cp 0,
\end{aligned}
\end{displaymath}
where the second term converges to zero by continuity of $\psi$ and
$\sqrt{n}$-consistency of $\bh$. Theorem 1 is thus proved.
%\end{proof}

\medskip
\begin{center}
Proof of Theorem 2
\end{center}

%\begin{proof}
\noindent\textit{Proof. }Let $\eta=H\gamma$, we can write the log local partial likelihood function
in terms of $\beta$ and $\eta$
\begin{displaymath}
\begin{aligned}[e]
l_n({\beta,\eta},\tau,v)=&\int_0^\tau \frac{1}{n}\sum_{i=1}^n K_h(
\bZ_i- v)\biggl[(\beta^T
\MB{Z}_i)^T H^{-1}\eta \\
&-\log \bigl\{\sum_i^n Y_i(u)\exp\{(\beta^T \MB{Z}_i)^T
H^{-1}\eta\}K_h (\bZ_i- v ) \bigr\}\biggr]dN_i(u).
\end{aligned}
\end{displaymath}


Accordingly let $S_{h,0}(\beta,\eta,u,v)=\frac{1}{n}\sum_{i=1}^n
K_h(\bZ_i-v)Y_i(u)\exp\{(\beta^T \MB Z_i)^T H^{-1}\eta\},$ \\
and
$S_{h,1}(\beta,\eta,u,v)=\frac{1}{n}\sum_{i=1}^n
K_h(\bZ_i-v)Y_i(u)\exp\{(\beta^T \MB Z_i)^T H^{-1}\eta\}(\beta^T \MB
Z_i)^T H^{-1},$ \\
and for  a $\sqrt{n}$-consistent estimate $\bh$ of
$\bT$, Lemma 1 implies
\begin{displaymath}
\sup_{0\le u<\tau}
\biggl|\frac{S_{h,1}(\bh,\eta,u,v)}{S_{h,0}(\bh,\eta,u,v)}
-\nu_1\biggr|\cp0,
\end{displaymath}
  where $\nu_1=\int {\MB
u}K(u)du$.

The derivative of $l_n(\beta,\eta,\tau,v)$ with respect to $\eta$
evaluated at $\bh$ and $\eT=H\gT$ is
\begin{displaymath}
\begin{aligned}[b]
&l'_n(\bh,\eT,\tau,v)\\
=&\int_0^\tau \frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i- v)\biggl[(\bh^T
\MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}\biggr]
dN_i(u)\\
=&\int_0^\tau \frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i- v)\biggl[(\bh^T
\MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}\biggr]
dM_i(u)\\
+&\int_0^\tau \frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i- v)\biggl[(\bh^T
\MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}\biggr]\\
&\qquad \qquad \qquad \qquad \qquad \qquad \times
Y_i(u)\exp\{\psi(\bTZ_i)\}\lambda_0(u)du\\
\equiv& U_n(\bh,\eT,\tau,v)+B_n(\bh,\eT,\tau,v).
\end{aligned}
\end{displaymath}

The first term
\begin{displaymath}\label{Un}
\begin{aligned}[e]
&U_n(\bh,\eT,\tau,v)\\ =&\int_0^t \frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i-
v)\biggl[(\bh^T \MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}\biggr]
dM_i(u)\\
=&\int_0^t \frac{1}{n}\sum_{i=1}^n K_h(\bTZ_i- v)\biggl[(\bT^T
\MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bT,\eT,u,v)}{S_{h,0}(\bT,\eT,u,v)}\biggr]
dM_i(u)\\
+&\int_0^t \frac{1}{n}\sum_{i=1}^n \Big(K_h(\bhZ_i- v)-K_h(\bTZ_i- v)\Big)\\
&\qquad \qquad \qquad \qquad \times \biggl[(\bT^T \MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bT,\eT,u,v)}{S_{h,0}(\bT,\eT,u,v)}\biggr]
dM_i(u)\\
+&\int_0^t \frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i- v)\\
&\qquad \qquad \times \biggl[(\bh^T
\MB{Z}_i-\bT^T \MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}
-\frac{S_{h,1}(\bT,\eT,u,v)}{S_{h,0}(\bT,\eT,u,v)}\biggr]
dM_i(u)\\
\equiv &U_n(\bT,\eT,\tau,v)+V+VI.
\end{aligned}
\end{displaymath}

It is clear that $\sqrt{nh} U_n(\bT,\eT,t,v)$ is a martingale with
predictable variation
\begin{displaymath}\label{Un}
\begin{aligned}[b]
&\langle
\sqrt{nh}U_n(\bT,\eT,t,v),\sqrt{nh}U_n(\bT,\eT,t,v)\rangle\\
=&\frac{nh}{n^2}\sum_{i=1}^n \int_0^\tau K_h^2(\bTZ_i-
v)\biggl[(\bT^T \MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bT,\eT,u,v)}{S_{h,0}(\bT,\eT,u,v)}\biggr]^{\otimes
2}\\
&\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \times Y_i(u)\exp\{\psi(\bTZ_i)\}\lambda_0(u)du.\\
=&f(v)\exp\{\psi(v)\}\Lambda(t,v)\int K^2(u)(\MB u-
\nu_1)^{\otimes 2}
 du +o_p(1) \equiv\Sigma_U(t,v)+o_p(1),
\end{aligned}
\end{displaymath}
where  the last step follows from Lemma 1.

\begin{comment}
since as n$\rightarrow\infty$
$$I\{\sqrt{h/n}\bigl|K_h(\bTZ_i-v) Q_{n,i_l}(u)\big|>\epsilon\}
=0\quad uniformly\; in \;u. \quad$$
\end{comment}
The Lindberg conditions are satisfied (see \cite{Wang:2001} for
details), we have thus proven that
\begin{equation}\label{AN:Un}
\sqrt{nh}U_n(\bT,\eT,\tau,v)\rightarrow_D
N\bigl(0,\Sigma_U(\tau,v)\bigr).
\end{equation}

As for the term $V$ and $VI$, similarly to the proof of Theorem 1,
we have
\begin{equation}\label{V}
\begin{aligned}[e]
V=&\int_0^\tau \frac{1}{n}\sum_{i=1}^n \Big(K_h(\bhZ_i- v)-K_h(\bTZ_i-
v)\Big) \\
&\qquad \qquad \qquad \qquad \times \biggl[(\bT^T \MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bT,\eT,u,v)}{S_{h,0}(\bT,\eT,u,v)}\biggr] dM_i(u)\\
&=
O_p(\frac{1}{nh^2}),
\end{aligned}
\end{equation}
%and similarly as shown in the proof of Theorem 1,
and
\begin{equation}\label{IV}
\begin{aligned}[e]
VI=&\int_0^\tau \frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i- v)\\
&\qquad \times \biggl[
(\bh^T\MB{Z}_i-\bT^T \MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}
-\frac{S_{h,1}(\bT,\eT,u,v)}{S_{h,0}(\bT,\eT,u,v)}\biggr]
dM_i(u)\\
=& O_p(\frac{1}{nh}).
\end{aligned}
\end{equation}

\begin{comment}
Next note that
\begin{displaymath}
\begin{aligned}[b]
&\int_0^\tau \frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i-v) \biggl[(\bh^T
\MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}\biggr]Y_i(u)\exp\{\psi(v)+(\bh^T
\MB{Z}_i)^T H^{-1}\eT\}\lambda_0(u)du =0,
\end{aligned}
\end{displaymath}

so we can rewrite
\begin{displaymath}\label{Un}
\begin{aligned}[b]
B_n(\bh,\eT,\tau,v)=&\int_0^\tau \frac{1}{n}\sum_{i=1}^n
K_h(\bhZ_i-v)\biggl[(\bh^T \MB{Z}_i)^T H^{-1}
-\frac{S_{h,1}(\bh,\eT,u,v)}{S_{h,0}(\bh,\eT,u,v)}\biggr]\\
\qquad & \times
Y_i(u)\bigl[\exp\{\psi(\bTZ_i)\}-\exp\{\psi(v)+(\bh^T \MB{Z}_i)^T
H^{-1}\eT\}\bigr]\lambda_0(u)du.
\end{aligned}
\end{displaymath}

By Taylor's expansion, we get
\begin{displaymath}
\begin{aligned}[b]
&\exp\{\psi(\bTZ_i)\}-\exp\{\psi(v)+(\bh^T \MB{Z}_i)^T H^{-1}\eT\}\\
=&\exp\{\psi(\bhZ_i)\}-\exp\{\psi(v)+(\bh^T \MB{Z}_i)^T
H^{-1}\eT\}+\big[\exp\{\psi(\bTZ_i)\}-\exp\{\psi(\bhZ_i)\}\big]\\
=&\exp\{\psi(v)\}\frac{\psi^{(p+1)}(v)}{(p+1)!}(\bhZ_i-v)^{p+1}+o_p(h^{p+1})+O_p(\frac{1}{\sqrt{n}}).\\
\end{aligned}
\end{displaymath}
\end{comment}

Applying Lemma 1 again and by Taylor expansion we get
\begin{equation}\label{Bn}
\begin{aligned}
B_n(\bh,\eT,\tau,v)
=&f(v)\exp\{\psi(v)\}\frac{\psi^{(p+1)}(v)}{(p+1)!}\Lambda(\tau,v)
\int  K(u)(\MB u- \nu_1) u^{p+1} du h^{p+1}\\
&\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad +o_p(h^{p+1})+O_p(\frac{1}{\sqrt{n}})\\
 =&b(\tau,v)+o_p(h^{p+1})+O_p(\frac{1}{\sqrt{n}}).
\end{aligned}
\end{equation}

We have thus shown that, under (\ref{AN:Un}), (\ref{V}),
(\ref{IV}) and (\ref{Bn}),
\begin{equation}\label{AN:FD}
\sqrt{nh} l'_n(\bh,\eT,\tau,v)\rightarrow_D
N\Big(b(\tau,v),\Sigma_U(\tau,v)\Bigr).
\end{equation}

Next we focus on the property of the second derivative
$l''_n(\bh,\eta,t,v)$. \begin{comment}
Let
$S_{h,2}(\beta,\eta,u,v)=\frac{1}{n}\sum_{i=1}^n
K_h(\bZ_i-v)Y_i(u)\exp\{(\beta^T \MB Z_i)^T H^{-1}\eta\}
 (\beta^T \MB Z_i) (\beta^T \MB Z_i)^T,$ we can write
\begin{displaymath}\label{l2}
\begin{aligned}[b]
&l''_n(\bh,\eta,t,v)\\
=&-\int_0^t \frac{S_{h,2}(\bh,\eta,u,v)S_{h,0}(\bh,\eta,u,v)
-S_{h,1}(\bh,\eta,u,v)^{\otimes 2} }
{S_{h,0}(\bh,\eta,u,v)^2}\frac{1}{n}\sum_{i=1}^n K_h(\bhZ_i-v)d
N_i(u).
\end{aligned}
\end{displaymath}

Again by Lemma 1, if we let $A=\int \MB u \MB u^T
K(u)du-\nu_1\nu_1^T$, then $$\sup_{0\leq
u\leq\tau}\biggl|\frac{S_{h,2}(\bh,\eta,u,v)S_{h,0}(\bh,\eta,u,v)
-S_{h,1}(\bh,\eta,u,v)^{\otimes 2} } {S_{h, 0}(\bh,\eta,u,v)^2}-A
\biggr| \cp 0.$$

So,
\begin{displaymath}
\begin{aligned}[b]
&-l''_n(\bh,\eT,\tau,v)\\
=&A \int_0^\tau\frac{1}{n}\sum_{i=1}^nK_h(\bhZ_i-v)d
N_i(u)+o_p(1)\\
=&A \frac{1}{n}\sum_{i=1}^nK_h(\bTZ_i-v)dN_i(\tau)+A
\frac{1}{n}\sum_{i=1}^n\big[K_h(\bhZ_i-v)-K_h(\bTZ_i-v)\big]dN_i(\tau)+o_p(1)\\
=&A f(v)E\{N(\tau)|\bTZ=v\}+o(\frac{1}{\sqrt{nh^4}})+o_p(1)\\
=&A f(v)\exp\{\psi(v)\}\Lambda(\tau,v)+o_p(1)\\
\equiv& \Sigma_l(\tau,v)+o_p(1).
\end{aligned}
\end{displaymath}
\end{comment}
Let $\eh=H\gh$, by Taylor expansion and Lemma 1 we have
\begin{equation}\label{TE:gamma}
0=l'_n(\bh,\eh,\tau,v)=l'_n(\bh,\eT,\tau,v)+l''_n(\bh,\eta^{**},\tau,v)(\eh-\eT),
\end{equation}
where $\eta^{**}$ lies in between $\eh$ and $\eT$. Theorem 1 implies
$\eh\cp \eT$,  hence $\eta^{**}\cp \eT$.

Using condition ({\it C1}) and boundedness of $\bhZ$, we arrive at
%it is easy to show that, for $r=0,1$ and 2,
%\begin{displaymath}
%\sup_{0\le u\le \tau}
%\big|S_{h,r}(\bh,\eta^{**},u,v)-S_{h,r}(\bh,\eT,u,v)\big|\cp0.
%\end{displaymath}

\begin{equation}\label{Con:SD}
l''_n(\bh,\eta^{**},\tau,v)=l''_n(\bh,\eT,\tau,v)+o_p(1)
=\Sigma_l(\tau,v)+o_p(1).
\end{equation}

By (\ref{AN:FD}), (\ref{TE:gamma}), (\ref{Con:SD}) and Slutsky's
theorem,
\begin{displaymath}
\begin{aligned}[b]
\sqrt{nh}(\eh-\eT)
&=\sqrt{nh}\biggl[-l''_n(\bh,\eta^{**},\tau,v)^{-1} l'_n(\bh,\eT,\tau,v)\biggr]+o_p(1)\\
&\cd
N\Big(b(\tau,v),\Sigma_l(\tau,v)^{-1}\Sigma_U(\tau,v)\Sigma_l(\tau,v)^{-1}\Big).\\
\end{aligned}
\end{displaymath}

Simple calculations lead to the result in Theorem 2.
%\end{proof}



\begin{thebibliography}{0}
\bibitem{Breslow:1974}
\textsc{Breslow, N.} (1974). 
Covariance analysis of censored survival data.
\textit{Biometrics}, \textbf{30}, 89--99. 

\bibitem{CLW:1999}
\textsc{Chen, C. H. and Li, K. C.} and \textsc{Wang, J. L.} (1999). 
Dimension reduction for censored regression data.
\textit{The Annals of Statistics}, \textbf{27}, 1--23. 

\bibitem{CWY:1995}
\textsc{Cheng, S. C., Wei, L. J.} and \textsc{Ying, Z.} (1995). 
Analysis of transformation models with censored data.
\textit{Biometrika}, \textbf{82}, 835--845. 

\bibitem{Cox:1972}
\textsc{Cox, D. R.} (1972). 
Regression models and life-table
(with discussion).
\textit{Journal of the Royal Statistcal Society, Ser.}, 
\textbf{B. 4}, 187--220. 

\bibitem{Cox:1975}
\textsc{Cox, D. R.} (1975). 
Partial likelihood.
\textit{Biometrika}, \textbf{62}, 269--276. 

\bibitem{DD:1988}
\textsc{Dabrowska, D. M.} and \textsc{Doksum, K. A.} (1988). 
Partial likelihood in transformation models with censored data.
\textit{Scandinavian Journal of Statistics}, \textbf{15}, 1--23. 

\bibitem{Doksum:1987}
\textsc{Doksum, K. A.} (1987). 
An extension of partial likelihood methods for proportional hazard models to general transformation
models.
\textit{The Annals of Statistics}, \textbf{15}, 325--345. 

\bibitem{FGK:1997}
\textsc{Fan, J., Gijbels, I.} and \textsc{King, M.} (1997). 
Local likelihood and local partial likelihood in hazard regression.
\textit{The Annals of Statistics}, \textbf{25}, 1661--1690. 

\bibitem{FH:1991}
\textsc{Fleming, T. R.} and \textsc{Harrington, D. P.} (1991). 
\textit{Counting Processes and Survival Analysis}. John Wiley \& Sons, Inc., New York.

\bibitem{GC:1991}
\textsc{Gentleman, R.} and \textsc{Crowley, J.} (1991). 
Local full likelihood estimation for the proportional hazards model.
\textit{Biometrics}, \textbf{47}, 1283--1296. 

\bibitem{Gray:1992}
\textsc{Gray, R. J.} (1992). 
Flexible methods for analyzing survival data using splines, with applications
    to breast cancer prognosis.
\textit{Journal of the American Statistical Association}, \textbf{420}, 942--951. 

\bibitem{HT:1986}
\textsc{Hastie, T.} and \textsc{Tibshirani, R.} (1986). 
Generalized additive models.
\textit{Statistical Science}, \textbf{3}, 297--318. 

\bibitem{HT:1990}
\textsc{Hastie, T.} and \textsc{Tibshirani, R.} (1990). 
Exploring the nature of covariate effects in the proportional hazards model.
\textit{Biometrics}, \textbf{46}, 1005--1016. 

\bibitem{HL:2008}
\textsc{Hosmer, D. W., Lemeshow, S.} and \textsc{May, S.} (2008). 
\textit{Applied Survival Analysis: Regression Modeling of Time to Event Data:
Second Edition}. John Wiley \& Sons, Inc., New York. 

\bibitem{HL:2006}
\textsc{Huang, J.} and \textsc{Liu, L.} (2006). 
Polynomial spline estimation and inference of proportional hazards
    regression models with flexible relative risk form.
\textit{Biometrics}, \textbf{62}, 793--802. 

\bibitem{OS:1988}
\textsc{O'Sullivan, F.} (1988). 
Nonparametric estimation of relative risk using splines and cross-validation.
\textit{SIAM Journal on Scientific and Statistical Computing}, \textbf{9}, 531--542. 

\bibitem{SH:1990}
\textsc{Sleeper, L. A.} and \textsc{Harrington, D. P.} (1990). 
Regression splines in the Cox model with application to covariate effects in
liver disease.
\textit{Journal of the American Statistical Association}, \textbf{85}, 941--949. 

\bibitem{TH:1987}
\textsc{Tibshirani, R.} and \textsc{Hastie, T.} (1987). 
Local likelihood estimation.
\textit{Journal of the American Statistical Association}, \textbf{82}, 559--567. 

\bibitem{Wang:2001}
\textsc{Wang, W} (2001). 
Proportional hazards model with unknown link function and applications
             to time-to-event longitudinal data.
\textit{Ph.D. Thesis}. 

\bibitem{Wang:2004}
\textsc{Wang, W} (2004). 
Proportional hazards regression with unknown link
function and time-dependent covariates.
\textit{Statistica Sinica}, \textbf{14}, 885--905. 

\end{thebibliography}


\end{document}
