---
abstract: 'We propose a novel approach to estimate the Cox model with temporal covariates. Our new approach treats the temporal covariates as arising from a longitudinal process which is modeled jointly with the event time. Different from the literature, the longitudinal process in our model is specified as a bounded variational process and determined by a family of Initial Value Problems associated with an Ordinary Differential Equation. Our specification has the advantage that only the observation of the temporal covariates at the time to event and the time to event itself are required to fit the model, while it is fine but not necessary to have more longitudinal observations. This fact makes our approach very useful for many medical outcome datasets, like the New York State’s Statewide Planning and Research Cooperative System (SPARCS) and the National Inpatient Sample (NIS), where it is important to find the hazard rate of being discharged given the accumulative cost but only the total cost at the discharge time is available due to the protection of patients’ information. Our estimation procedure is based on maximizing the full information likelihood function. The resulting estimators are shown to be consistent and asymptotically normally distributed. Variable selection techniques, like Adaptive LASSO, can be easily modified and incorporated into our estimation procedure. The oracle property is verified for the resulting estimator of the regression coefficients. Simulations and a real example illustrate the practical utility of the proposed model. Finally, a couple of potential extensions of our approach are discussed.'
author:
- 'Xiaoqi Zhang[^1]'
- 'Xiaobing Zhao[^2]'
- 'Yanqiao Zheng[^3]'
title: A novel approach to estimate the Cox model with temporal covariates and application to medical cost data
---

[***Keywords—*** Cox regression; longitudinal process; joint model; maximum full likelihood; adaptive LASSO; Gaussian process; semi-martingale]{}

Introduction
============

In the proportional hazards model ([@cox1972regression; @andersen1982cox]), the hazard function of the event time $T$ takes the form

$$\label{cox-proportional-hazard}
\lambda\left(t\mid Z\right)=\lambda_{0}\left(t\right)\exp\left(b_{0}^{T}Z\right)$$

where $\lambda\left(t\mid Z\right)$ is the conditional hazard function of $T$ given the $p\times1$ covariate vector $Z$, $\lambda_{0}\left(t\right)$ is an unspecified baseline hazard function and $b_{0}$ is a $p\times1$ vector of unknown regression coefficients. Although in the original paper, the covariate $Z$ is viewed as random vectors independent from time $t$, the model can be easily extended to the case where $Z$ is time-dependent and $Z\left(t\right)$ is assumed to be given through an unknown stochastic process.

The main-stream procedure used to estimate the model is the maximum partial likelihood (MPL) procedure, which applies well whether or not the covariates are time-dependent. However, when time-dependent covariates being involved, the consecutive observation of the covariates is required in the sense that every subject must have its covariates observed at all the failure time prior to its own failure. In the other words, let $T_i<T_j$ be the observed failure time for two different subjects $i$ and $j$, then for the dying later subject $j$, its covariate $Z_j$ must have values observed at both of $T_i$ and $T_j$. Otherwise, the MPL procedure won’t work. Although some approximation methods were proposed to relax consecutive observation requirement as discussed in [@andersen1992repeated], the MPL procedure can’t be applied effectively to the medical cost datasets, like the New York State’s Statewide Planning and Research Cooperative System (SPARCS), or the National Inpatient Sample (NIS), where it is important to find the hazard rate of being discharged given the accumulative cost while only the total cost for each inpatient observed at the discharge time is available.

In this paper, we shall propose a novel estimation procedure for the model , which can generate consistent estimates for parameters and the baseline hazard in the model even if only the observations $\left\{\left(Z_{i,T_i},T_i\right):i=1,\dots,n\right\}$ are available, where all the observations $Z_{i,T_j}$ with $T_j<T_i$ are missing. Our procedure is based on joint modelling the longitudinal process that generates the time-dependent covariates and the time to event. The topic of joint model has been widely discussed ([@henderson2000joint; @song2002semiparametric; @hsieh2006joint; @ye2008semiparametric; @rizopoulos2011dynamic; @kim2013joint; @lawrence2015joint]). A comprehensive review is also available in [@tsiatis2004joint; @sousa2011review; @ibrahim2010basic]. However, all these works require the condition that the number of observations of the longitudinal measurement before the time to event is greater than the dimension of the longitudinal process, which restricts their usefulness for the input data with only one observation of the covariate value at the event time. In this paper, we propose an alternative specification of the joint model. Formally, we only assume that the longitudinal measurements follow a bounded variational process that can be expressed as a stochastic integral as below: $$\label{potential-growth-process}
\mathcal{Z}(t):=Z_{0}+\int_{0}^{t}Y(s)\epsilon(s)ds.$$ where $Z_{0}$ represents the initial value and $Y(t):=I\left(T>t\right)$. The model consists of two components:

\(1) the longitudinal process: $$\label{longitudinal}
Z(t)=Z_0+\int_{0}^{t}\epsilon(s)ds$$ which characterizes the evolution of the longitudinal measurements, we assume the conditional expectation of the increment rate $\epsilon(t)$ has the following parametric form $$\label{conditional-expectation-of-growth-rate}
q\left(z,t\mid a_{0}\right) := E\left(\epsilon(t)\mid Z(t)=z,a_{0}\right);$$ (2) the event process $Y$ which determines the time to event and its conditional expectation has the form , i.e. $$\label{conditonal-hazard}
\lambda\left(t\mid z\right)=
\lambda_{0}\left(t\right)\exp\left(b_{0}^{T}z\right):=
E\left(-dY(t)\mid Y(t^{-})=1, Z(t^{-})=z\right).$$ It turns out by [@XQZhang2017], combining and yields a complete specification of the joint model in the sense that if two joint models share a common pair of the conditional hazard function and the conditional expectation , all the distributions in interest arising from the two joint models are identical. The equation is the key to derive an explicit expression of the joint probability density function (pdf) of $Z_T$ and $T$ which helps design our estimation procedure. To our best knowledge, there has not been any previous works attempting to specify the longitudinal process as . We hope our work could provide some hints to the future development of this field.

In model and , there are three sets of parameters, $a_{0}$, $b_{0}$ and $\lambda_{0}$. Among them, $\lambda_{0}$ has infinite dimension. Our procedure will estimate the three kinds of parameters through maximizing the full information likelihood function, where the likelihood is constructed from the joint probability density function (pdf) of $Z_T$ and $T$. By [@XQZhang2017], this joint pdf is expressed as below by using the function $q$ and $\lambda$:

$$\label{joint-pdf}
{\scalebox{0.8}{$
pdf\left(z,t\mid a_{0},b_{0},\lambda_{0}\right)=\tilde{p}\left(z,t\right)\times\exp\left(-\int_{0}^{t}\lambda_{0}\left(t-s\right)\exp\left(b_{0}^{T}g\left(z,t,s\right)\right)ds
\right)\times\lambda_{0}\left(t\right)\exp
\left(b_{0}^{T}z\right)$}}$$

where the function $\tilde{p}$ is the time-dependent pdf induced by the longitudinal process $\left\{ Z_{t}\right\} $ and by [@XQZhang2017] it can be expressed as $$\label{p-tilde}
\tilde{p}(z,t)=p\left(g\left(z,t,t\right),0\right)\cdot \mathcal{J}_{z,t},$$ The function $p(,.0)$ is the initial pdf induced by $Z_0$ and for every $t$, $\mathcal{J}_{z}(t)$ denotes the Jacobian of the function $g\left(.,t,t\right)$ evaluated at the point $z$. The function $g$ is solely determined by $q$ through solving a family of initial value problems (IVPs). Namely for every fixed $z$ and $t$, $g\left(z,t,.\right)$ is the solution to the following ordinary differential equation (ODE) for $s\in\left(0,t\right)$: $$\label{ode}
z'\left(s\right) = -q\left(z\left(s\right),t-s\mid a_{0}\right)$$ subject to the initial condition $g\left(z,t,0\right)=z$.

In addition to estimating model parameters, in practice it is also important to select the significant covariates. Variable selection approaches has been extensively studied by many authors. The least absolute shrinkage and selection operator (LASSO) was presented by [@tibshirani1996regression]. [@fan2001variable] developed the nonconcave penalized approach (SCAD) for variable selection which applies to likelihood-based estimation procedures, including the MPL procedure for the Cox model ([@fan2002variable]). [@zou2006adaptive] developed an Adaptive LASSO approach and showed its oracle property under a general set of conditions. The estimation procedure proposed in the current paper can be easily combined with those variable selection approaches. In particular, we will incorporate a modified version of the Adaptive LASSO into our procedure and verify its oracle property.

The rest of this paper is organized as follows. In Section 2, we will sketch the estimation procedures in detail. The large sample properties of resulting estimators are stated in Section 3. Simulation results and the application to real world data are presented in Section 4. Section 5 discusses some extensions of our model and concludes. All proofs are collected in Appendix.

ESTIMATION PROCEDURE
====================

The estimation procedure is based on maximizing the full information likelihood function which is formed as the product of joint pdf of the failure time $T\in[t,t+dt)$ and the observation of the longitudinal measure $Z_{t}$ at time $t$. To deal with the non-parametric $\lambda_{0}$, we adopt the method that approximates $\lambda_{0}$ through a sequence of finite dimensional step-wise functions, denoted as $\lambda_{n}^{s}$, with the number of steps given by the sample size $n$.

Throughout this section, we assume the data input for the estimation procedure only has the observation at the time to event, such as $\left\{\left(Z_{T_i},T_i\right):\,i=1,\dots,n\right\}$, while in the remark section, we will briefly discuss the adjustment of our procedure to deal with the case where more longitudinal observations are available.

Likelihood Function {#likelihood function}
-------------------

Define $A\subset\mathbb{R}^{d}$ as the domain of all possible values of the parameter $a$, $a_{0}\in A$ as the true parameter. Similarly, Define $B\subset\mathbb{R}^{p}$ as the domain of $b$, $b_{0}\in B$ as the true parameter. For every fixed $a\in A$, define $g\left(.\mid a\right)$ as the solution trajectories to the IVPs conditional on $a$. When the analytic form of $g\left(.\mid a\right)$ is not available, we can use its numerical approximation in place. There are many efficient numerical solvers to the IVPs . In this paper, we pick up the Euler’s method \cite{} for the purpose of being simple and illustrative. Similarly, write $\mathcal{J}_{z|a}(t)$ as the Jacobian of $g\left(.,t,t\mid a\right)$ and when necessary it can be replaced by its numerical version.

By [@XQZhang2017], the initial pdf $p(.,0)$ is uniquely determined by the function and the joint pdf . In particularly, given the joint pdf , there is a well defined map $a\mapsto p_{a}(.,0)$ from the parameter space $A$ to the space of the pdfs over $\mathbb{R}^{p}$. Therefore, given the input data and a fixed parameter $a$, we can estimate $p_a(.,0)$ by the Gaussian kernel density method as below: $$\label{initial-appro}
p_{a,n}(z,0):=\frac{1}{n}\sum_{i=1}^{n}G_{n^{-1/4}}\left(z-g\left(z_{t_i},t_i,t_i\mid a\right)\right)$$ where $G_h$ denote the Gaussian kernel function with kernel width $h$. In this paper, we simply select the kernel width as $\frac{1}{n^4}$ as can guarantee the function converges to $p_{a}(.,0)$ in the $L^1$ norm for all $a$ \cite{}.

For the baseline hazard $\lambda$, without loss of generality, we set $t_{0}=0$, $\theta_{0}=1$ and let $0\leq t_{1}<t_{2}<\cdots<t_{n}<\infty$ be the ordered statistics of the $n$ observed failure time. A step-wise version of the non-parametric baseline hazard is constructed as below: $$\label{step-wise-lambda}
\lambda_{n}^{s}\left(t\right):=\sum_{i=1}^{n}\theta_{i}\cdot I_{[t_{i-1},t_{i})}\left(t\right)$$ where $\theta=\left\{ \theta_{i}\geq0:i=1,\dots,n\right\} $ is a set of parameters to be estimated. For each profile of the parameters $\Omega_{n}=\left(\lambda_{n}^{s},a,b\right)$, we can define the log likelihood as below: $$\label{likelihood1}
l_n\left(a,b,\lambda_{n}^{s}\right):=\frac{1}{n}\sum_{i=1}^{n}\left(
\begin{aligned}
&\log p_{a,n}\left(g\left(z_{t_i},t_i,t_i\mid a\right),0\right)
+\log \mathcal{J}_{z_{t_i}|a}\left(t_i\right)+\log\lambda_{n}^{s}\left(t_{i}\right)\\
&+b^{T}z_{t_{i}}-\int_{0}^{t_{i}}\lambda_{n}^{s}\left(t_{i}-\tau\right)\exp\left(b^{T}g\left(z_{t_{i}},t_{i},\tau|a\right)\right)d\tau
\end{aligned}
\right).$$ The estimator resulting from maximizing is denoted as $\hat{\lambda}^{s}$ and $\hat{b}$ and $\hat{a}$.

\[remark1\] The first order condition of the optimization problem indicates the relation $$\label{constraint-for-lambda}
\theta_{i}=\frac{1}{\sum_{j=i+1}^{n}\int_{0}^{t_{i}-t_{i-1}}\exp\left(b^{T}g\left(z_{t_{j}},
t_{j},t_{j}-t_{i}+\tau|a\right)\right)d\tau}$$ at the optimal point $\hat{\lambda}_{n}^{s}$ and $\hat{b}$ and $\hat{a}$. Relation can be inserted as a set of constraints back into the optimization problem , which helps sharply reduce the dimension of the original problem.

Variable Selection  {#variable selection}
-------------------

A penalty function can be naturally incorporated into the estimation and determine the non-zero component of the coefficients $b$ automatically. There exist many types of penalty functions in the literature, such as LASSO ([@tibshirani1996regression]), adaptive LASSO ([@zou2006adaptive]), SCAD ([@fan2001variable]), MCP. For convenience of calculation, we choose the adaptive LASSO penalty function defined by:

$$\label{penality-funcion}
P_{\Lambda}\left(b_{j}\right):=\Lambda\frac{\left|b_{j}\right|}{\left|\hat{b}_{j}\right|^{r}}$$

where $b_{j}$ is the $j$-th component of the the vector $b$, $\hat{b}_{j}$ is a consistent estimator of $b_{j}$ (say the estimator generated by maximizing Eq. ), $\Lambda$ is a tuning parameter and $r>0$ (for simplicity of calculation, we let $r=2$ throughout the paper). The penalty function is added into the likelihood function and forms the following maximization problem : $$\label{likelihood4}
\max_{a,u,\lambda^{s}}\left(l_{n}\left(a,\hat{b}+u,\lambda^{s}\right)+\sum_{j=1}^{p}P_{\Lambda_{n}}
\left(\hat{b}_{j}+\frac{u}{\sqrt{n}}\right)\right)$$ where the choice of tuning parameter $\Lambda_{n}$ depends on the sample size $n$. Denote $\hat{a}_{AL}$, $\hat{\lambda}_{AL}^{s}$ and $\hat{b}_{AL}$ ($:=\hat{b}+\hat{u}_{AL}$ with $\hat{u}_{AL}$ being the solution to ) as the estimator corresponding to maximize Eq. .

In the original work ([@zou2006adaptive]), a wide range of candidate values can be selected for the tuning parameter $\Lambda_n$ as long as the asymptotic property $\Lambda_n\cdot n^{-\frac{1}{2}}\rightarrow 0$ and $\Lambda_n\cdot n^{\frac{1}{2}}\rightarrow \infty$ hold. In practice, there are multiple ways to select the value of $\Lambda_n$ such as the Bayesian criterion method \cite{} which is based on iterative calculation of the values of the Bayesian criterion and the likelihood function. Instead of using the iterative procedure, we simply set $\Lambda_n=n^{-\frac{1}{4}}$ in order to reduce the computation load.

The original adaptive LASSO method is designed for the OLS procedure ([@zou2006adaptive]), but it turns out this method applies very well to the likelihood-based estimation procedures as ours. The oracle property of the adaptive LASSO estimator, $\hat{b}_{AL}$, will be proved in Appendix \[theorem-ALASSO\].

LARGE SAMPLE PROPERTIES
=======================

The consistency and asymptotic normality of the estimators, $\hat{a}$, $\hat{b}$ and $\hat{\lambda}_{n}^{s}$, constructed in section 2.1 will be established in this section. We will also show the oracle property of the adaptive LASSO estimator $\hat{b}_{AL}$ and the consistency and asymptotic normality of $\hat{\lambda}_{AL}^{s}$ and $\hat{a}_{AL}$.

Large Sample Properties of $\hat{a}$, $\hat{b}$ and $\hat{\lambda}^{s}$ {#large sample property}
-----------------------------------------------------------------------

Let $\Omega:=\left(a,b,\lambda\right)$ be a given profile of parameter and in particular, $\Omega_{0}$ be the true parameter. Define $$\label{func V}
{\scalebox{0.8}{$
V_{\Omega}:=p_{a}\left(g\left(Z_{T},T,s\mid a\right),0\right)\mathcal{J}_{Z_T|a}(T)\exp\left(-\int_{0}^{T}\lambda\left(T-s\right)\exp\left(b^{T}g\left(Z_{T},T,s\mid a\right)\right)ds\right)\lambda\left(T\right)
\exp\left(b^{T}Z_{T}\right)
$}}$$ as a random variable with $Z_{T}$ and $T$ being the random variables following the joint pdf associated with $\Omega_{0}$, denoted as $pdf_{\Omega_0}$. The following technical conditions are needed for the consistency result:

$C1$. The domain $A$ and $B$ are compact. $B$ has open interior with $b_{0}\in B^{\circ}$. The domain of $\lambda$, denoted as $L$, has $\lambda_0\in L$ and is a set of uniformly bounded right-continuous functions satisfying that $\lambda(0)=1$ for all $\lambda\in L$ (which means $L\subset L^{\infty}\left([0,\infty)\right)$ and under the weak-$\ast$ topology, the closure of $L$ is compact).

$C2$. $E_{\Omega_{0}}\left(\left|\log\left(V_{\Omega}\right)\right|\right)$, $E_{\Omega_{0}}\left(\left|\nabla_{b_{i}}\log\left(V_{\Omega}\right)\nabla_{b_{j}}\log\left(V_{\Omega}\right)\right|\right)$, $E_{\Omega_{0}}\left(\left|\nabla_{b_{i}b_{j}}\log\left(V_{\Omega}\right)\right|\right)$ are finite for all $i,j=1,\dots,p$ and all $\Omega\in A\times B\times L$; and the matrix $I=\left\{ E_{\Omega_{0}}\left(\nabla_{b_{i}}\log\left(V_{\Omega}\right)\nabla_{b_{j}}\log\left(V_{\Omega}\right)\right)\right\} _{1\leq i,j\leq p}$ and $H=\left\{ E_{\Omega_{0}}\left(\nabla_{b_{i}b_{j}}\log\left(V_{\Omega}\right)\right)\right\} _{1\leq i,j\leq p}$ are positive definite.

$C3$. There exists a positive function $d\in L^{1}(pdf_{\Omega_0})-$ such that for all $\Omega\in A\times B\times L$, $\left|\log V_{\Omega}\right|\leq d\left(Z_T,T\right)$ almost surely with respect to the probability measure $pdf_{\Omega_0}$. $C4$. For all $a\in A$, $q\left(z,t\mid a\right)\in C^{2}\left(\mathbb{R}^{p}\times\mathbb{R}_{+}\right)$ and the map given through $q\left(.\mid.\right):A\rightarrow C^{2}\left(\mathbb{R}^{p}\times\mathbb{R}_{+}\right)$ is continuous with respect to the $C^{2}$ topology.

$C5$. For every $a\in A$, there is an $p\times p$ matrix $M_a$, such that $\frac{\partial q\left(z,0\mid a\right)}{\partial z}\rightarrow M_a$ as $\left\Vert z\right\Vert\rightarrow \infty$ ($\left\Vert.\right\Vert$ is the Euclidean Norm of a vector). And for different $a$ and $a'$, $M_a-M_{a'}$ has at least one eigenvalue with non-zero real part.

$C6$. $pdf_{\Omega_0}$ has the full support $\mathbb{R}^{p}\times \mathbb{R}_{+}$. The true initial $p(.,0)$ satisfies that $\int_{\mathbb{R}^{p}}\exp(c\cdot z)p(z,0)dz\not=1$ for every $c\in \mathbb{R}^{p}$. Condition $C1$-$C3$ are standard for the consistency and asymptotic normality of maximum likelihood estimators. $C4$ is the regularity condition that guarantees the trajectories $g\left(.\mid a\right)$ depends on $a$ smoothly. $C5$ is the key to guarantee the identification of the model and , although it turns out that $C5$ can be discarded without any impact on the consistency of $\hat{b}$ and $\hat{\lambda}$, and both of $C5$ and $C6$ can be discarded when the event time and the longitudinal process satisfy a kind of Markovian property and the extra longitudinal observations are available. We will go back to these extensions in the section \[extension\].

\[theorem-identification\] Under Condition $C5$ and $C6$, model and are identifiable. And $E_{\Omega_{0}}\left(\log\left(V_{\Omega}\right)\right)$ has the unique maximal point, $\Omega_{0}$. In addition, if $C4$ holds, $E_{\Omega_{0}}\left(\log\left(V_{\Omega}\right)\right)$ is continuous with respect to the variable $\Omega$.

\[theorem-consistency\] (1). Under $C1$-$C6$, the estimator $\hat{a}$, $\hat{b}$ are consistent and $\hat{b}-b_{0}\rightarrow_{d}N\left(0,I^{-1}\right)$;

(2). the estimator $\hat{\lambda}^{s}$ converges to $\lambda_{0}$ according to the weak-$*$ topology and $\sqrt{n}\left(\int_{0}^{t}\hat{\lambda}^{s}\left(\tau\right)-\lambda_{0}\left(\tau\right)d\tau\right)$ converges weakly to a Gaussian Process;

\[theorem-ALASSO\] Under $C7$, the estimator $\hat{a}_{AL}$ and $\hat{\lambda}_{AL}^{s}$ has the same properties as $\hat{a}$ and $\hat{\lambda}^{s}$ as stated in Theorem 2, and the estimator $\hat{b}_{AL}$ has the following oracle property:

(1). denote $\mathcal{A}\subset\left\{ 1,\dots,p\right\} $ as the set of indices with $b_{0,j}\not=0$ for $j\in\mathcal{A}$ and $\hat{\mathcal{A}}\subset\left\{ 1,\dots,p\right\} $ as the set of indices with $\hat{b}_{0,j}\not=0$ for $j\in\hat{\mathcal{A}}$, then $\hat{b}_{j}\rightarrow_{p}b_{0,j}$ for all $j\in\mathcal{A}$ and $Prob\left(\hat{\mathcal{A}}=\mathcal{A}\right)\rightarrow1$;

(2). denote $I_{\mathcal{A}}:=\left\{ E_{\Omega_{0}}\left(\nabla_{b_{i}}\log\left(V_{\Omega}\right)\nabla_{b_{j}}\log\left(V_{\Omega}\right)\right)\right\} _{i,j\in\mathcal{A}}$, $\hat{b}_{L,\mathcal{A}}=\left(\hat{b}_{L,j}\right)_{j\in\mathcal{A}}$ and $b_{0,\mathcal{A}}=\left(b_{0,j}\right)_{j\in\mathcal{A}}$, $\sqrt{n}\left(\hat{b}_{L,\mathcal{A}}-b_{0,\mathcal{A}}\right)\rightarrow_{d}N\left(0,I_{\mathcal{A}}^{-1}\right)$ .

Extension
---------

When the longitudinal observations are available at the observation time before failure occurs at $T$, i.e. the input data has the form $\left\{\left(Z_{i,j},t_{i,j}\right):\,j=1,\dots,m_i;\,i=1,\dots,n\right\}$ with $m_{i}>1$ for $i=1,\dots,n$. A two-step procedure can be applied to estimate the parameter $\Omega_{0}$, and the resulting estimator turns out to be consistent and have asymptotically normal distribution even without the assumption $C5$ and $C6$. Instead, the following Markovian-style condition are required: $$\label{markovian}
E\left(\epsilon(t-s)\mid Z(t),T\geq t\right)\equiv E\left(\epsilon(t-s)\mid Z(t)\right)$$ where $\epsilon(t)$ is the instantaneous variational rate of the longitudinal process as specified in model model . Eq. implies that the conditional mean trajectory that reaches a given realization, $Z(t)$, at the observational time $t$ won’t be affected by whether or not the event has already occurred. Formally, the two-step algorithm is stated as following:

Step 1: estimate the parameter $a$ through minimizing the empirical mean of the $L^{2}$ distance between the empirical longitudinal trajectories observed for each individual $i$ and the theoretical mean trajectories passing through the point $\left(Z_{t_{i,m_{i}}},t_{i,m_i}\right)$:

$$\label{mean-L2-distance}
\min_{a\in A}\frac{1}{n}\sum_{i=1}^{n}\frac{1}{m_i}\sum_{j=1}^{m_{i}}\left(g\left(z_{t_{i,m_{i}}},t_{i,m_{i}},t_{i,m_{i}}-t,_{i,j}\mid a\right)-z_{t_{i,j}}\right)^{2}.$$

It turns out when $n\rightarrow\infty$, the estimator $\hat{a}_{E}$ of solving the problem is consistent.

Step 2: replace $a$ by $\hat{a}_{E}$ and maximize the likelihood function or for the parameter $b$ and $\lambda$.

It turns out that the resulting estimators do have the same properties as stated in Theorem \[theorem-consistency\] or \[theorem-ALASSO\]. The two-step procedure separates the estimation of $a_{0}$ from the estimation of $b_{0}$ and $\lambda_{0}$. Thanks to this separation, in the second step, the initial pdf and the Jacobian can be completely removed from the likelihood function or because they only depends the parameter $a$ and the fixed underlying true distribution. In the other words, once if the parameter $a$ is replaced by its first-step estimator $\hat{a}_E$, the component of $\tilde{p}$ becomes constants, and can be deleted from the second-step maximization problem without any impact on the final estimators. As a consequence, there is no need to calculate the Jacobian $\mathcal{J}_{z_{t_{m_i}}|a}\left(t_{i,m_i}\right)$, which makes the two-step procedure running much faster than the original procedure because the computation of the Jacobian is the most time-consuming part.

Due to the fact that the conditional density function $$\rho(z,t):=\exp\left(-\int_{0}^{t}\lambda\left(t-s\right)
\exp\left(b^{T}g\left(z,t,s\mid a\right)\right)ds\right)\lambda\left(t\right)
\exp\left(b^{T}z\right)$$ has already had the full support $\mathbb{R}^{p}\times \mathbb{R}_{+}$, and the second-step optimization is equivalent to the optimization of a conditional log-likelihood function formed by the sum of $\log\rho\left(z_{t_{m_i}},t_{m_i}\right)$ for all $m_i$’s, the full support condition in $C6$ can be relaxed.

In addition, from the proof of the theorem \[theorem-identification\], it is clear that the main difficulty to achieve the injectivity of the map from the parameter space to the space of all joint pdfs lies in the exclusion of the possibility that different $a$ could lead to the same joint pdf, which is exactly the condition $C5$ and $C6$ designed for. In contrast, when $a_0$ is given, the identification of $b_0$ and $\lambda_0$ becomes trivial and doesn’t require any further conditions like $C5$ and $C6$. So when the estimation of $a_0$ can be separated out, $C5$ and $C6$ are redundant.

Proof of the validity of the two-step procedure is in Appendix \[proof-two-step\]. A latent assumption behind the proof is that the observational time is uninformative and the total number of observations, $m_i$, is at least $2$ for all subjects. Unlike the joint model discussed in [@tsiatis2004joint], we don’t have to assume $m_i$ greater than the dimension of the covariates. This fact makes our two-step procedure more attractive to the scenarios where there are only a few longitudinal observations but a large set of covariates.

Numerical Studies
=================

Simulation Studies {#simulation}
------------------

In this section, simulation studies are conducted to evaluate the finite-sample performance of the estimation procedures proposed in section \[likelihood function\]. Consider the following examples:

\[example1\] 50 samples, each consisting of $n=400,\,800$ subjects, are generated from simulating the version of model that has covariate dimension $p=16$, coefficients $b_0=(1,1,-1,0,\dots,0)$ with $3$ non-zero covariate effects. Given $$a=(1,0.5,-1,0.3,1,0.5,-1,0.3,1,0.5,-1,0.3,1,0.5,
-1,0.3),$$ the conditional expectation function is specified as the constant function as below: $$\label{parametrized-conditional-growth-rate}
q\left(z,t\mid a\right)=a.$$ The baseline hazard is specified through the function $$\label{lambda0}
\lambda_0(t)=\frac{\exp(10)+\exp^{-t}}{\exp(10)+1}.$$ The initial $Z_0\sim N(0,I_{16})$ with $I_{16}$ being the $16$-dimensional identity matrix.

The simulation results are presented in the terms of the following criteria:

\(1) Figure \[fig: fitting\] shows the the estimated cumulative hazard $\int_{0}^{t}\hat{\lambda}_{n}^{s}(\tau)d\tau$ versus the true cumulative hazard \[lambda0\] for Example \[example1\]. The bias and standard deviation of the estimated non-zero $b_0$ are given in Table \[table: non-zero param\]. (2) We also conduct variable selection for both of the example \[example1\] by the adaptive LASSO method , the result are summarized in Table \[table: variable selection\] for Example \[example1\]. The result are reported by:

i\. The average number of the true zero coefficients of $b_0$ that are correctly set to zero, denoted by $C(b_0)$.

ii\. The average number of the true non-zero coefficients of $b_0$ that are incorrectly set to zero, which is given by $IC(b_0)$.

iii\. The proportion of samples that excluding any non-zero coefficients, denoted by $U-$fit.

iv\. The proportion of samples selecting the exact subset models (correct-fit) and the proportion of smaples including all the variables (over-fit), labeled by $C-$fit and $O-$fit respectively.

From Table \[table: non-zero param\], it is clear that for both of the two cases $N=400$ and $N=800$, the fitting to the non-zero coefficients are very good while the fitting accuracy in the case of $N=800$ is even better. As for variable selection, Table \[table: variable selection\] shows that in most of the samples, the set of zero variables can be exactly identified by our procedure. In particular, as the sample size increases, the identification accuracy is risen up as well. Even in the rare samples where some zero variables are misclassified, the misclassification happened sparsely as $C(b_0)>12$ and that value is close to the true value, $13$.

Real Example {#real example}
------------

The New York State’s Statewide Planning and Research Cooperative System (SPARCS) 2013 is a system initially created to collect information on discharges from hospitals within New York State. SPARCS currently collects patient level detail on patient characteristics, diagnoses and treatments, services, and charges for each hospital inpatient stay and outpatient visit; and each ambulatory surgery and outpatient services visit to a hospital extension clinic and diagnostic and treatment center licensed to provide ambulatory surgery services. In 2013, the SPARCS contains nearly 2.5 million inpatient discharges from 218 facilities and 58 counties in New York State. Patient demographics in the SPARCS include age group at admission, gender, race, source of payment and zip code. Patient clinical characteristics include type of admission, diagnosis codes (MDC code, DRG code, CCS diagnosis code etc.) and treatment procedures undergone (CCS Procedure Code).

An important property of the SPARCS data is that there is not any other longitudinal observation available for time-dependent variables, like the cumulative charge, than the observation at the discharge time. Therefore, neither the traditional maximum partial likelihood method nor the estimation procedures designed for the joint models as discussed in [@kim2013joint; @zeng2007maximum] can be well applied to the SPARCS data. In contrast, the approach proposed in this paper can effectively address the data issue as it is designed for.

In this paper, we consider the discharge time $T$, with the time-dependent covariate, the logarithm of the cumulative charge $Z_1$, and the stationary covariates consisting of the categorical variables, $Z_2,\dots,Z_{25}$, associated with 25 Major Diagnosis Code (MDC) and the degree ($1\sim 4$) of the Severity of Illness, $Z_{26}$. Our analysis is conducted on a subsample of the entire SPARCS 2013 database with sample size $400$. The summary statistics of our subsample are presented in Table \[table: 3.1\].

The penalized maximum likelihood estimators $\hat{b}_{AL}$ are reported in Table \[table: real\]. The non-parametric estimator $\int_0^t\hat{\lambda}^{s}_{n}(s)ds$ for the cumulative baseline hazard are plotted in Figure \[fig: real data study\].

In Table \[table: real\], the significant negative coefficients for log-charge indicates the strong positive correlation between the total charge and los. In addition, it seems that there does not exist robust connection between the los and the severity/mortality of illness.

By Figure \[fig: real data study\], the day 5 seems to be relatively special because the variation of the slope of the cumulative hazard turns from increasing to decreasing around this time, which implies that for patients who have already stayed in hospital for 5 days, they are more probable to have a longer stay.

Remarks and Conclusion {#conclusion}
======================

In this paper, we proposed a maximum full information likelihood procedure to estimate the Cox model with temporal covariates. The most significant advantage of our procedure is that it can generate well-performed estimation without requiring the extra longitudinal observations before the time to event. There are also three potential extensions to the current work.

Censoring
---------

Although censoring is not discussed in the current framework, it can be added in the standard way such that censoring is (1) independent from the occurrence of the interested event, or (2) conditionally independent from the event given the covariates at the observational time. In both of the two cases, the consistency and asymptotic normality of the resulting estimators still hold and their proof is straightforward from the proof of Theorem \[theorem-consistency\] and \[theorem-ALASSO\].

Forecast Long Term Survival Rate
--------------------------------

In addition to the hazard function , the estimators proposed in section \[likelihood function\] indicates a consistent estimator to the long term survival rate (LTSR): $$\label{LTSR}
S\left(z,t,t'\right):=Prob\left(T\in \left[t,t+t'\right)\mid Z_t=z\right).$$ Using the notation $g^{-1}$ as in Eq. , the estimator to can be given as below: $$\label{hat-LTSR}
\hat{S}\left(z,t,t'\right):=\exp\left(
-\int_{t}^{t+t'}\exp\left(\hat{b}^{T}g^{-1}\left(
z,t,\tau\mid \hat{a}\right)\right)\hat{\lambda}_n^s
\left(\tau\right)d\tau\right),$$ where $\hat{a}$, $\hat{b}$ and $\hat{\lambda}_n^s$ are the estimators derived in section \[likelihood function\], which can be replaced by their penalized version as well. The consistency and asymptotic normality of the estimator is just a direct result of the theorem \[theorem-consistency\] and/or \[theorem-ALASSO\]. It is worthwhile to mention that is not possible to be constructed from the maximum partial likelihood estimators of the Cox model when temporal covariates are included. Because it is clear from that $\hat{S}$ relies on the information of the temporal covariates $Z$ within the forecast interval $\left[t,t+t'\right)$, which is not available from the maximum partial likelihood estimators.

Semi-martingale Longitudinal Processes
--------------------------------------

Although in the current discussion, the longitudinal process is assumed to have bounded variation and absolutely continuous with respect to the Lebesgue measure on $\mathbb{R}_+$, the same framework should be extensible to more general cases where the longitudinal process may not have bounded variation (for example, given by a semi-martingale process). In a series of related works, the authors construct an explicit expression of the joint pdf of the event time and a semi-martingale longitudinal process, which enables us to construct the full information likelihood function. But the challenges to extend the current work to the situation with the semi-martingale longitudinal measurements are the identification of the resulting model and the challenge in computation. For the identification issue, it is clear from the proof \[proof-identification\] that in the current framework, the identification relies on detailed analysis of the solution trajectories of ODE system induced by the function . In the case of semi-martingale longitudinal measurements, the ODE system will be replaced with a more complicated partial differential equation (PDE) system. Although it seems that there is no barrier to make the same trick in proof \[proof-identification\] invalid, the details to transplant the proof \[proof-identification\] to the semi-martingale case is open to future studies. In the aspect of computation, we have to apply numerical method to a PDE system in place of an ODE system, while,as known, the numerical method to solve PDE system is much more time consuming. A potential solution to the computation issue is to utilize the relation between PDE systems and the semi-martingale processes, through which simulating the underlying process could yield exactly the same solution to the PDE problem. The details of implementing that idea are left as another open problem for further study.

Proof for Theorem \[theorem-identification\] {#proof-identification}
--------------------------------------------

The identifiability of the model and is equivalent to that as long as $\Omega\not=\Omega_0$, $$\label{unique-maximum}
{\scalebox{0.9}{$
\begin{aligned}
&p_{a}\left(g\left(z,t,t\mid a\right),0\right)\mathcal{J}_{z|a}(t)\exp(-\int_{0}^{t}\exp(b^Tg\left(z,t,s\mid a\right))\lambda(t-s)ds)\lambda(t)\exp\left(b^{T}z\right)\\
&\not=
p_{a_0}\left(g\left(z,t,t\mid a_0\right),0\right)\mathcal{J}_{z|a_0}(t)\exp(-\int_{0}^{t}\exp(b_0^{T}g\left(z,t,s\mid a_0\right))\lambda_0(t-s)ds)\lambda_0(t)\exp\left(b_0^{T}z\right)
\end{aligned}$}}$$ within a positive measure set $M\in\mathbb{R}^{p}\times\mathbb{R}_{+}$ (with respect to the standard Lebesgue Measure).

\(1) Suppose does not hold for some $\Omega$ with $b\not=b_0$. The right-continuity condition in $C1$ requires that $$\label{initial identity}
p_a(z,0)\exp(b^Tz)\equiv p(z,0)\exp(b_0^Tz)$$ where $p(.,0)$ is the true initial pdf. The assumption that the true joint pdf has the full support implies that $p(.,0)$ has full support as well. Therefore, Eq. leads to $$p_a(z,0)\equiv p(z,0)\exp\left(\left(b_0-b\right)^Tz\right),$$ both $p_a(.0)$ and $p(.,0)$ are probability density function, which yield that $$\int_{\mathbb{R}^{p}}p(z,0)\exp(c\cdot z)dz=1$$ for some $c\not=0$ that contradicts to the requirement in $C6$. Consequently, every $\Omega$ that could potentially break down the condition must have $b=b_0$.

\(2) On the other hand, if $b=b_0$, we have for all $t\geq 0$: $$\label{condition-for-lambda}
S_{(a,b_0,\lambda)}(t)\cdot \lambda(t)=S_{(a_0,b_0,\lambda_0)}(t)\cdot \lambda_0(t)$$ where $S_{(a,b,\lambda)}$ is the survival function of the failure time that follows the joint pdf associated with $\Omega=(a,b,\lambda)$, by definition it has the following form: $${\scalebox{0.9}{$
S_{(a,b,\lambda)}(t):\int_{\mathbb{R}^p}p_{a}\left(g\left(z,t,t\mid a\right),0\right)\mathcal{J}_{z|a}(t)\exp(-\int_{0}^{t}\exp(b_0^Tg\left(z,t,s\mid a\right))\lambda(t-s)ds)dz\lambda(t)$}}.$$ Because the survival function is uniquely determined by the pdf of the event time which is furthermore uniquely determined by the joint pdf. Under the assumption that $\Omega$ and $\Omega_0$ corresponds to exactly the same joint pdf, the equation enforces that $\lambda=\lambda_0$ for all $\Omega$ that breaks the condition .

\(3) Suppose there exists $\Omega=(a,b_0,\lambda_0)$ with $a\not=a_0$ for which the condition doesn’t hold. Then, the following identity holds: $$\label{weak-identity}
{\scalebox{0.9}{$
\begin{aligned}
&p_{a}\left(z,0\right)\mathcal{J}_{g^{-1}\left(z,0,t\mid a\right)|a}(t)\exp(-\int_{0}^{t}\exp(b^Tg^{-1}\left(z,0,s\mid a\right))\lambda(s)ds)\\
&=
p_{a_0}\left(z_0,0\right)\mathcal{J}_{g^{-1}\left(z_0,0,t\mid a_0\right)|a_0}(t)\exp(-\int_{0}^{t}\exp(b_0^{T}g^{-1}\left(z_0,0,s\mid a_0\right))\lambda_0(s)ds).
\end{aligned}$}}$$ for all pairs $(z,z_0)$ such that $z_0=g\left(g^{-1}\left(z,0,t\mid a\right),t,t\mid a_0\right)$ where $g^{-1}\left(z,s,t\mid a\right)$ is the inverse trajectories of $g$ and defined through the relation $$\label{g-inverse}
g\left(g^{-1}\left(z,s,t\mid a\right),s+t,t\mid a\right)=z.$$ Factor out Eq. by $\mathcal{J}_{g^{-1}\left(z_0,0,t\mid a_0\right)|a_0}(t)\exp(-\int_{0}^{t}\exp(b_0^{T}g^{-1}\left(z_0,0,s\mid a_0\right))\lambda_0(s)ds)$ and take the limit as $t\rightarrow 0$ yielding the following identity: $$\label{invariant-measure}
p_a(\mathcal{T}_{r}(z_0),0)\cdots\mathcal{J}_{\mathcal{T}_r}(z_0)=p(z,0)$$ where for every $r\in \mathbb{R}$, the map $\mathcal{T}_{r}:\mathbb{R}^{p}\rightarrow \mathbb{R}^{p}$ is the diffeomorphism obtained from solving the ODE system: $$\label{ode2}
z'=q\left(z,0\mid a\right)-q\left(z,0\mid a_0\right),$$ $\mathcal{T}_r(z_0)$ is just the point reached at the time $r$ by the trajectory starting at $z_0$ that solves Eq. . $\mathcal{J}_{\mathcal{T}_r}$ is the Jacobian associated with $\mathcal{T}_r$. By the language of ergodic theory, Eq. implies that the probability measure $p(.,)$ is invariant under the $\mathbb{R}$-action on the space $\mathbb{R}^{p}$ induced by the solutions $\mathcal{T}$. However, under the condition $C5$, the action $\mathcal{T}$ associated with the pair of $a$ and $a_0$ does not allow any invariant probability measure fully supported on $\mathbb{R}^{p}$ unless $a=a_0$. This contradiction guarantees the condition .

The uniqueness of the maximal point $\Omega_{0}$ of function $E_{\Omega_{0}}\left(\log\left(V_{\Omega}\right)\right)$ is simply the consequence of the standard proof the consistency of the full information maximum likelihood estimator, and can be found in every advanced textbook of econometrics, like [@amemiya1985advanced].

The continuity of $E_{\Omega_{0}}\left(\log\left(V_{\Omega}\right)\right)$ with respect to $\Omega$ and its differentiability with respect to the component $b$ comes from $C4$ by the the dominant convergent theorem. This completes the proof for Theorem \[theorem-identification\].

Proof for Theorem \[theorem-consistency\] {#proof-consistency}
-----------------------------------------

The relation defines a map, denoted as $\hat{\iota}_n$, that assigns every $(a,b)\in A\times B$ a step-wise function with the step heights specified by . We can define the asymptotic version of $\hat{\iota}_n$ as below: $$\iota\left(a,b\right)(t):=
\frac{\int_{\mathbb{R}^{p}}pdf_{\Omega_{0}}
\left(z,t\right)dz}{\int_{t}^{\infty}\int_{\mathbb{R}^{p}}pdf_{\Omega_{0}}
\left(z,\tau\right)\cdot\exp\left(b^T g\left(z,\tau,\tau-t\mid a\right)\right)dzd\tau},$$ where $pdf_{\Omega_0}$ is the joint pdf associated with the true parameter $\Omega_0$. It turns out that $\iota$ is continuous with respect to the weak$-\ast$ topology on the $L^{\infty}$ space and has the compact domain $A\times B$. In addition, for every pair $(a,b)$, $\hat{\iota}_n(a,b)\rightarrow \iota(a,b)$ in the weak$-\ast$ topology.

Then, the consistency of $\hat{a}$, $\hat{b}$ and $\hat{\lambda}^{s}$ follows immediately from $C2$ and the facts: (1) the function $l_n\left(a,b,\hat{\iota}_n\left(a,b\right)\right)\rightarrow_p E_{\Omega_{0}}\left(\log\left(V_{\left(a,b,\iota\left(a,b\right)\right)}\right)\right)$ (by the strong law of large number); (2) the function $E_{\Omega_{0}}\left(\log\left(V_{\left(a,b,\iota\left(a,b\right)\right)}\right)\right)$ is continuous and has a unique maximal point, $\Omega_{0}$, within a compact domain (by the theorem \[theorem-identification\]).

The asymptotic normality of $\hat{b}$ can be verified by the standard argument for the asymptotic normality of a maximum likelihood estimator.

To verify the asymptotic normality of $\hat{\lambda}^{s}$, firstly notice that let $\left\{ Z_{t}:t\in[0,\infty)\right\} $ be a process satisfying model associated with the true parameter $a_{0}$ and $\left(Z_{T},T\right)$ be a random vector following the distribution associated with $\Omega_{0}$, denote $N\left(t\right):=I\left(T\leq t\right)$ being the counting process determined by $T$ and $\tilde{N}\left(t\right):=I\left(T> t\right)$. The processes $N\left(t\right)$ and $\tilde{N}(t)$ determines a martingale process as below:

$$\label{estimator-process}
M\left(t\right):=\int_{0}^{t}\frac{dN\left(s\right)}{E\left(\exp\left(b_{0}^{T}Z_{s}\right)\cdot \tilde{N}\left(s\right)\right)}-\frac{
\tilde{N}\left(t\right)}{E\left(\tilde{N}(t)\right)}\int_{0}^{t}\lambda_{0}\left(\tau\right)d\tau$$

it turns out that $E\left(M\left(t\right)\right)\equiv0$, $$Var\left(M\left(t\right)\right)=\int_{0}^{t}
\frac{\lambda_{0}\left(s\right)}{E\left(\exp\left(b_{0}^{T}Z_{s}\right)\cdot \tilde{N}\left(s\right)\right)}ds+\frac{\left(\int_{0}^{t}
\lambda_{0}\left(s\right)ds\right)^{2}}{E\left(\tilde{N}(t)\right)}$$ for $t>s$. On the other hand, by the relation in Eq. , we have:

$${\scalebox{0.8}{$
\begin{aligned}
\sqrt{n}\left(\int_{0}^{t}\hat{\lambda}^{s}\left(\tau\right)-\lambda_{0}\left(\tau\right)d\tau\right)  &=  \frac{1}{\sqrt{n}}\left(
\begin{aligned}
&\sum_{i=1}^{j_{t}}\frac{n\cdot\left(t_{i}-t_{i-1}\right)}{\sum_{j=i+1}^{n}\int_{0}^{t_{i}-t_{i-1}}\exp\left(b^{T}g\left(z_{t_{j}},
t_{j},t_{j}-t_{i}+\tau|a\right)\right)d\tau}\\
&+\frac{n\cdot\left(t-t_{j_{t}}\right)}{\sum_{j=j_t+1}^{n}\int_{0}^{t_{j_{t}}-t_{j_{t}-1}}\exp\left(b^{T}g\left(z_{t_{j}}
,t_{j},t_{j}-t_{j_t}+\tau|a\right)\right)d\tau}-\frac{n-j_t}{\frac{n-j_t}{n}}\int_{0}^{t}\lambda_{0}\left(\tau\right)d\tau
\end{aligned}
\right)\\
 & =  \frac{1}{\sqrt{n}}\cdot\sum_{i=1}^{n}\left(\frac{I\left(t_{i}<t\right)}{\frac{\sum_{j=i+1}^{n}\exp\left(b^{T}g\left(z_{t_{j}},
t_{j},t_{j}-t_{i}|a\right)\right)}{n}}-\frac{I\left(t_{i}\geq t\right)}{\frac{n-j_t}{n}}\int_{0}^{t}\lambda_{0}\left(\tau\right)d\tau\right)+
 O\left(h_n\right)\\
 & =  \frac{\sum_{i}^{n}M_{i}\left(t\right)}{\sqrt{n}}+O\left( h_n\right)
\end{aligned}$}}$$

where $j_{t}=\max\left\{ i\in\left\{ 1,\dots,n\right\} :t_{i}<t\right\} $ and $h_n=\max\left(h_n^1,h_n^2,h_n^3\right)$, with $h_n^1$, $h_n^2$ and $h_n^3$ given as below: $$\begin{aligned}
h_n^1 &:=\sup\left\{ t_{i}-t_{i-1}:i=1,\dots,n\right\}\\
h_n^2 &:=\sup\left\{\left|\frac{\left(n-i\right)}{n}-E\left(\tilde{N}\left(t\right)\right)\right|:\, i\leq j_{t}\right\}\\
h_n^3 &:={\scalebox{0.9}{$\sup\left\{\left|\frac{\sum_{j=i+1}^{n}\exp\left(b^{T}g\left(z_{t_{j}},
t_{j},t_{j}-t_{i}|a\right)\right)}{n}-E\left(\exp\left(b_{0}^{T}Z_{s}\right)\cdot \tilde{N}\left(s\right)\right)\right|:i\leq j_t\right\}$}}
\end{aligned}$$ By the assumption that $\lambda_{0}$ is strictly positive, the fact that $\sup\left\{ t_{i}-t_{i-1}:i=1,\dots,n\right\} \rightarrow_{p}0$, and that for every fixed $t<\infty$, $\left|\frac{\left(n-i\right)}{n}-E\left(1-N\left(t\right)\right)\right|\rightarrow_p 0$, ${\scalebox{0.75}{$\left|\frac{\sum_{j=i+1}^{n}\exp\left(b^{T}g\left(z_{t_{j}},
t_{j},t_{j}-t_{i}|a\right)\right)}{n}-E\left(\exp\left(b_{0}^{T}Z_{s}\right)\cdot \tilde{N}\left(s\right)\right)\right|\rightarrow_p 0$}}$ uniformly by uniform law of large number. Therefore, by central limit theorem, we have:

$$\label{limit-eq1}
\lim_{n\rightarrow\infty}\sqrt{n}\left(\int_{0}^{t}\hat{\lambda}^{s}\left(\tau\right)-\lambda_{0}\left(\tau\right)d\tau\right)=\lim_{n\rightarrow\infty}\frac{\sum_{i}^{n}\left(M_{i}\left(t\right)-\int_{0}^{t}\lambda_{0}\left(\tau\right)d\tau\right)}{\sqrt{n}}=N\left(0,Var\left(M\left(t\right)\right)\right)$$

Applying a vector version of the central limit theorem as well as Eq. , we can prove the weak convergence of $\sqrt{n}\left(\int_{0}^{t}\hat{\lambda}^{s}\left(\tau\right)-\lambda_{0}\left(\tau\right)d\tau\right)$ to the Gaussian Process $B\left(Var\left(M\left(t\right)\right)\right)$.

Finally, the consistency and the weak convergence of the estimator $_{t}\hat{S}_{t'}$ is direct from the consistency and asymptotic normality of $\hat{\lambda}^{s}$.

Proof for Theorem \[theorem-ALASSO\] {#proof-alasso}
------------------------------------

Construct a function $\Psi_{n}\left(a,u,\lambda^{s}\right)$ with $u\in\mathbb{R}^{p}$ as below:

$$\label{Psi-n}
\Psi_{n}\left(a,u,\lambda^{s}\right):=l_{n}\left(a,\hat{b}+u,\lambda^{s}\right)-\Lambda_{n}\sum_{i=1}^{p}\frac{\left|\hat{b}_{j}+\frac{u_{j}}{\sqrt{n}}\right|}{\left|\hat{b}_{j}\right|^{2}}$$

Define $W_{n}\left(a,u,\lambda^{s}\right):=\Psi_{n}\left(a,u,\lambda^{s}\right)-\Psi_{n}\left(a,0,\lambda^{s}\right)$, then $W_{n}$ has the following form: $$\label{W-n}
W_{n}\left(a,u,\lambda^{s}\right) = l_{n}\left(a,\hat{b}+u,\lambda^{s}\right)-l_{n}\left(a,\hat{b},\lambda^{s}\right)-\frac{\Lambda_{n}}{\sqrt{n}}\sum_{i=1}^{p}\frac{\sqrt{n}\left(\left|\hat{b}_{j}+\frac{u_{j}}{\sqrt{n}}\right|-\left|\hat{b}{}_{j}\right|\right)}{\left|\hat{b}{}_{j}\right|^{2}}.$$

If $b_{0,j}\not=0$, then $\hat{b}_{j}\rightarrow_{p}b_{0,j}$ and $\sqrt{n}\left(\left|\hat{b}_{j}+\frac{u_{j}}{\sqrt{n}}\right|-\left|\hat{b}_{j}\right|\right)\rightarrow u_{j}\cdot\frac{b_{0,j}}{\left|b_{0,j}\right|}$. Therefore we have $\frac{\Lambda_{n}}{\sqrt{n}}\frac{\sqrt{n}\left(\left|\hat{b}_{j}+\frac{u_{j}}{\sqrt{n}}\right|-\left|\hat{b}_{j}\right|\right)}{\left|\hat{b}_{j}\right|^{2}}\rightarrow_{p}0$ by the assumption $\frac{\Lambda_{n}}{\sqrt{n}}\rightarrow0$. If $b_{0,j}=0$, $\sqrt{n}\left(\left|\hat{b}_{j}+\frac{u_{j}}{\sqrt{n}}\right|-\left|\hat{b}_{j}\right|\right)=\left|u_{j}\right|$ and $\frac{\Lambda_{n}}{\sqrt{n}\left|\hat{b}_{j}\right|}=\sqrt{n}\Lambda_{n}\frac{1}{\left|\sqrt{n}\hat{b}_{j}\right|^{2}}$, as $\sqrt{n}\hat{b}_{j}=O\left(1\right)$ as $n\rightarrow\infty$ almost surely, we have $\frac{\Lambda_{n}}{\sqrt{n}\left|\hat{b}_{j}\right|}\rightarrow_{p}\infty$ by the assumption $\sqrt{n}\Lambda_{n}\rightarrow\infty$. Consequently,

$$\label{W}
{\scalebox{0.8}{$
W_{n}\left(a,u,\lambda^{s}\right)\rightarrow_{p}W\left(a,u,\lambda^{s}\right)=\begin{cases}
E_{\Omega_{0}}\left(\log\left(V_{\left(a,b_{0}+u,\lambda^{s}\right)}\right)\right)-E_{\Omega_{0}}\left(\log\left(V_{\Omega_{0}}\right)\right) & \textrm{if }u_{j}=0\textrm{ for all }j\in\mathcal{A}\\
-\infty & \textrm{else}
\end{cases}$}}$$

which implies all the estimators, $\hat{a}_{AL}$, $\hat{b}_{AL}$ and $\hat{\lambda}_{AL}^{s}$, converges to the unique maximal points of the function $E_{\Omega_{0}}\left(\log\left(V_{\left(a,b_{0}+u,\lambda^{s}\right)}\right)\right)-E_{\Omega_{0}}\left(\log\left(V_{\Omega_{0}}\right)\right)$, which is the true value $a_{0}$, $b_{0}$ and $\lambda_{0}$. The asymptotic normality of $\sqrt{n}\left(\hat{\lambda}_{AL}^{s}-\lambda_{0}\right)$ is established in the same way as in the proof of Theorem \[theorem-consistency\]. The oracle property of $\hat{b}_{AL}$ can be verified as below:

First, it is obvious that for each $j\in\mathcal{A}$, $\lim_{n\rightarrow\infty}Prob\left(j\in\hat{\mathcal{A}}\right)=1$. And for $j\not\in\mathcal{A}$, if $j\in\hat{\mathcal{A}}$, by the first order condition of the maximization problem , we have the identity that:

$$\label{identity-of-l-n}
\left|\frac{\partial l_{n}}{\partial_{u_{j}}}\right|=\frac{\Lambda_{n}}{\sqrt{n}\left|\hat{b}_{j}\right|^{2}}=\frac{\sqrt{n}\Lambda_{n}}{\left|\sqrt{n}\hat{b}_{j}\right|^{2}}$$

where the right-hand side $\rightarrow_{p}\infty$ while the left-hand side $\rightarrow_{p}\frac{\partial E_{\Omega_{0}}\left(\log\left(V_{\Omega_{0}}\right)\right)}{\partial_{b_{j}}}=0$ , while implies that $Prob\left(j\not\in\mathcal{A},j\in\hat{\mathcal{A}}\right)\leq Prob\left(\left|\frac{\partial l_{n}}{\partial u_{j}}\right|=\frac{\Lambda_{n}}{\sqrt{n}\left|\hat{b}_{j}\right|^{2}}\right)\rightarrow0$. It verifies that $\lim_{n\rightarrow\infty}Prob\left(\hat{\mathcal{A}}=\mathcal{A}\right)=1$.

Second, the asymptotic normality of $\hat{b}_{AL,\mathcal{A}}$. Consider

$$\label{derivative-identity-for-Psi-n}
\nabla_{u_{\mathcal{A}}}\Psi_{n}\left(a,0,\lambda^{s}\right) =  -\nabla_{u_{\mathcal{A}}}\nabla_{u_{\mathcal{A}}}\Psi_{n}\left(a,u',\lambda^{s}\right)u_{\mathcal{A}}$$

where the left-hand side equals to $\nabla_{u_{\mathcal{A}}}l_{n}\left(a,\hat{b},\lambda^{s}\right)-\left(\frac{\Lambda_{n}}{\sqrt{n}\left|\hat{b}_{j}\right|^2}\right)_{j\in \mathcal{A}}$ with $\sqrt{n}\left(\frac{\Lambda_{n}}{\sqrt{n}\left|\hat{b}_{j}\right|^2}\right)_{j\in \mathcal{A}}\rightarrow0$ and the right-hand side equals to $-\nabla_{u_{\mathcal{A}}}\nabla_{u_{\mathcal{A}}}l_{n}\left(a,u',\lambda^{s}\right)u_{\mathcal{A}}$ with $u'$ being some intermediate point between $0$ and $\hat{b}_{AL}-\hat{b}$. By the condition $C2$ and the consistency property of estimator $\hat{b}$, we have

$$\label{normality-conclusion-ALASSO}
\sqrt{n}\left(\hat{b}_{AL,\mathcal{A}}-\hat{b}\right)=\sqrt{n}u_{\mathcal{A}}\sim
\nabla_{u_{\mathcal{A}}}\nabla_{u_{\mathcal{A}}}
l_{n}\left(a,\hat{b}+u',\lambda^{s}\right)^{-1}\cdot\sqrt{n}
\nabla_{u_{\mathcal{A}}}l_{n}\left(a,\hat{b},\lambda^{s}\right)\rightarrow_{d}N\left(0,I_{\mathcal{A}}^{-1}\right).$$

Proof for the Validity of the Two-Step Procedure {#proof-two-step}
------------------------------------------------

It to show the estimator $\hat{a}_E$ from the first step is consistent, which is equivalent to show that the following function has a unique minimum: $$\label{mean_dis}
m(a):=E\left(E\left(\int_{0}^{t}\left(Z(s)-g\left(z,t,t-s\mid a\right)\right)^{2}ds\mid T=t,Z(t)=z\right)\right).$$ In fact, the unique minimal point of the function must be $a_0$ as long as $g\left(z,t,t-s\mid a\right)$ equals to the conditional mean of $E\left(Z(s)\mid Z(t)=z,T=t\right)$ for all $z$ and $s\leq t$, which is implied by the condition .

Tables & Figures
================

Tables
------

\[table: 3.1\]

\

\[table: non-zero param\]

\

\[table: variable selection\]

\

\[table: real\]

\

Figures
-------

![image](Simulation_Study2){width="15cm" height="7cm"} \[fig: fitting\]

\

![image](evaluating_fitting){width="7cm" height="4.5cm"} \[fig: real data study\]

\

[20]{} \[1\][\#1]{} \[1\][`#1`]{} urlstyle \[1\][doi: \#1]{}

Takeshi Amemiya. *Advanced econometrics*. Harvard university press, 1985.

Per Kragh Andersen. Repeated assessment of risk factors in survival analysis. *Statistical Methods in Medical Research*, 10 (3):0 297–315, 1992.

Per Kragh Andersen and Richard David Gill. Cox’s regression model for counting processes: a large sample study. *The annals of statistics*, pages 1100–1120, 1982.

David R Cox. Regression models and life-tables (with discussion). In *J. R. Statist. Soc.*, volume B, pages 187–220. 1972.

Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties. *Journal of the American statistical Association*, 960 (456):0 1348–1360, 2001.

Jianqing Fan and Runze Li. Variable selection for cox’s proportional hazards model and frailty model. *Annals of Statistics*, pages 74–99, 2002.

Robin Henderson, Peter Diggle, and Angela Dobson. Joint modelling of longitudinal measurements and event time data. *Biostatistics*, 10 (4):0 465–480, 2000.

Fushing Hsieh, Yi-Kuan Tseng, and Jane-Ling Wang. Joint modeling of survival and longitudinal data: likelihood approach revisited. *Biometrics*, 620 (4):0 1037–1043, 2006.

Joseph G Ibrahim, Haitao Chu, and Liddy M Chen. Basic concepts and methods for joint models of longitudinal and survival data. *Journal of Clinical Oncology*, 280 (16):0 2796–2801, 2010.

Sehee Kim, Donglin Zeng, Yi Li, and Donna Spiegelman. Joint modeling of longitudinal and cure-survival data. *Journal of statistical theory and practice*, 70 (2):0 324–344, 2013.

A Lawrence Gould, Mark Ernest Boye, Michael J Crowther, Joseph G Ibrahim, George Quartey, Sandrine Micallef, and Frederic Y Bois. Joint modeling of survival and longitudinal non-survival data: current methods and issues. report of the dia bayesian joint modeling working group. *Statistics in medicine*, 340 (14):0 2181–2195, 2015.

Dimitris Rizopoulos. Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data. *Biometrics*, 670 (3):0 819–829, 2011.

Xiao Song, Marie Davidian, and Anastasios A Tsiatis. A semiparametric likelihood approach to joint modeling of longitudinal and time-to-event data. *Biometrics*, 580 (4):0 742–753, 2002.

In[ê]{}s Sousa. A review on joint modelling of longitudinal measurements and time-to-event. *Revstat Stat J*, 9:0 57–81, 2011.

Robert Tibshirani. Regression shrinkage and selection via the lasso. *Journal of the Royal Statistical Society. Series B (Methodological)*, pages 267–288, 1996.

Anastasios A Tsiatis and Marie Davidian. Joint modeling of longitudinal and time-to-event data: an overview. *Statistica Sinica*, pages 809–834, 2004.

Wen Ye, Xihong Lin, and Jeremy MG Taylor. Semiparametric modeling of longitudinal measurements and time-to-event data–a two-stage regression calibration approach. *Biometrics*, 640 (4):0 1238–1246, 2008.

D Zeng and DY Lin. Maximum likelihood estimation in semiparametric regression models with censored data. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 690 (4):0 507–564, 2007.

Xiaoqi Zhang and John Ringland. A note on the joint probability density function of a stochastic growth process with random stopping time. 2017. URL <https://arxiv.org/abs/1701.04423>.

Hui Zou. The adaptive lasso and its oracle properties. *Journal of the American statistical association*, 1010 (476):0 1418–1429, 2006.

[^1]: Email: `xiaoqizh@buffalo.edu`

[^2]: Email: `maxbzhao@126.com`

[^3]: Email: `yanqiaoz@buffalo.edu`; Corresponding Author. Address: No. 18, Xueyuan Street, Xiasha Higher Education Park, Hangzhou, Zhejiang, 310018, China
