\section{Data Example}
\label{s: real data}


In this section, we apply the proposed method to a kidney post-transplantation mortality study. End-stage renal disease (ESRD) is one of the most deadly and costly diseases in the United States. From 2004 to 2016, ESRD incident cases increased from 345.6 to 373.4 per million people, with Medicare expenditures escalating from 18 to 35 billion dollars \citep[]{saran2019usrds}.
Kidney transplantation is the renal replacement therapy for  the majority of patients with ESRD. Successful kidney transplantation is associated with improved survival, improved quality of life, and health care cost savings when compared to dialysis.
However, despite aggressive effort to increase the number of donor kidneys, the demand far exceeds the supply of donor kidneys for transplantation and hence, the donor waiting list is very long. 
Currently about 130,000 patients are waiting for lifesaving organ transplants in the U.S., among whom 100,000 await kidney transplants and fewer than 15\% of patients will receive transplants in their lifetime. To optimize the organ allocation, further research is essential to
determine the risk factor associated with post-transplant mortality.


To better understand this problem, we considered the data obtained from the Organ Procurement and Transplantation Network (OPTN).
There were 146,248 patients who received transplants between 1990 and 2008. Failure time (recorded in years) was defined as the time from transplantation to graft failure
or death, whichever occurred first, where graft failure was considered to occur when the transplanted kidney ceased to function. Patient survival was censored 6 year post-transplant or at the end of study (2008). 
The median follow-up time was around 6 years and the censoring rate was 62\%. Covariates included in this study were age at transplantation, race, gender, cold ischemic time, donation after cardiac death (DCD), BMI, expanded criteria donor (ECD), dialysis time, comorbidity conditions such as glomerulonephritis, polycystic kidney disease, diabetes, and hypertension.
Detecting and accounting for time-varying effects are particularly important in the context of kidney transplantation, as
non-proportional hazards have already been reported in the literature \citep[]{Wolfe1999,He2017}. Also, analyses with time-varying effects provide valuable clinical information that could be obscured otherwise.

However, existing statistical software become computationally infeasible when fitting a time-varying effects model on a data set as large as what we have here. 
Thus, to estimate the potential time-varying effects, we fit the time-varying Cox model using the proposed sieve MLE, which is computationally scalable.
Specifically, based on previous studies, DCD, Polycystic, Diabetes and Hypertension are modeled with time-independent effects, and the remaining variables are estimated with time-varying effects. The time-varying effects are all implemented by cubic B-splines with 5 interior knots. Figure \ref{fig: baseline hazard} shows the estimated baseline hazard function. We can see that the post-transplant mortality is high in the short term after surgery, with a weakening association over time. 
Table \ref{tab:application} summarizes the estimated time-independent effects, and Figure \ref{fig: time-varying effects} shows examples of fitted time-varying effects with 95\% pointwise confidence intervals, where the standard error estimators were obtained by inverting the estimated information matrix of all parameters including time-independent coefficients and the coefficients of spline bases. As one can see, the effects of baseline age varied over time, resulting in an eventually strengthened association. Specifically, compared with the reference group (age at transplantation between 19-39), patients 40 to 49 years of age had a protective effect in the short term after transplantation.
We can also see that the high cold ischemic time is a risk factor for mortality in the short run, with a weakening association over time. Thus, special care should be dedicated to improve the short-term outcome.
As expected, longer waiting times on dialysis (greater than 5 years) negatively impact post-transplant survival, especially in the short run. Male gender was not significantly associated with mortality immediately after the renal transplantation but became a risk factor in the long run.  
As can be seen in Figure \ref{fig: time-varying effects}, underweight shows a protective effect in the short run, and then a slightly weakening association over time, which confirms the previous finding of \cite{Lafranca2015}.
The results regarding high BMI  should be interpreted with caution. Although higher levels of BMI in the general population are typically associated with high mortality, in chronic kidney diseases, such as patients with kidney dialysis and kidney transplantation, higher BMI has been associated with better survival, which has been labeled as reverse epidemiology \citep[]{Dekker2008, Kovesdy2010}.
Our results show that both overweight and obesity  improved survival in the short term after kidney transplantation, but obesity became a risk factor after long-term exposure.  One possible explanation is that BMI is a complex marker of visceral and nonvisceral adiposity and also of nutritional status including muscle mass \citep[]{Kovesdy2010}, and the improved short-term outcome associated with higher BMI may be related to differential benefits by one or more of these components.  Our findings indicate a need to critically reassess the role of BMI in the risk stratification of kidney transplantation. 
A further assessment (such as sub-group analysis) of high BMI that differentiates between visceral adiposity, nonvisceral adiposity and higher muscle mass may improve risk stratification in kidney transplant recipients. In addition, our
results show that graft survival for patients with Glomerulonephritis is better than patients with other primary diseases. 
Regarding racial disparities, the long-term survival outcomes for African Americans continue to lag behind non-African Americans.
Finally, as expected,  the effect
of expanded criteria donor (ECD) is not as good as optimal donor.  When a sub-optimal organ becomes available, patients and physicians must decide whether to accept the offer and special care  must be dedicated to improve the survival benefit. 



\begin{table}[!h]
   \begin{center}
   \caption{Summary of estimates for time-independent effects in kidney post-transplantation mortality study} 
   \label{tab:application}
   \begin{threeparttable}
   \begin{tabular}{ccccc}
  \toprule
       Variables & DCD & Polycystic & Diabetes & Hypertension\\
      \midrule
        EST & $-0.081$ & $-0.511$ & $0.333$ & $-0.146$ \\
        ESE & $0.038$ & $0.021$ & $0.012$ & $0.014$  \\
        95\% CI & $[-0.156, -0.007]$ & $[-0.553, -0.469]$ & $[\ \ 0.310,\ \ 0.357]$ & $[-0.172, -0.119]$ \\
        p-value & $0.033$ & $<0.001$ & $<0.001$ & $<0.001$ \\
        \bottomrule
    \end{tabular}
    \begin{tablenotes}
  \item[*]  \footnotesize{EST is the estimated time-independent effect, ESE is the estimated standard error by inverting the estimated information matrix of all parameters including the coefficients of spline basis, and CI is the confidence interval.
  }
  \end{tablenotes}
  \end{threeparttable}
  \end{center}
\end{table}

 \begin{figure}[!h]
     \centering
     \scalebox{0.55}
     {\includegraphics{fig/unos_baseline_summary.eps}}
     \caption{Estimated baseline hazard $\hat{\alpha}(t)$ using the the proposed sieve MLE method for the kidney transplantation data.}
     \label{fig: baseline hazard}
 \end{figure}
 


\begin{figure}
    \centering
    \scalebox{0.65}
    {\includegraphics{fig/unos_cox_selected_tv_summary.eps}}
    \caption{Estimated time-varying effects using the proposed sieve MLE method for the kidney transplantation data.}
    \label{fig: time-varying effects}
\end{figure}


\section{Theoretical Properties}
\label{s: asymptotics}

In this section, we study the theoretical properties of the proposed sieve MLE. Although many works have investigated  asymptotic distributional theories for M-estimation with bundled parameters \citep{ai2003efficient, chen2003estimation, ding2011}, their results cannot be directly applied to our setting. In particular,  the nuisance parameters in existing works often take the form of an unknown function of only some finite-dimensional Euclidian parameters of interest. However, our work focuses on a more general scenario, where the nuisance parameter is an unknown function of not only the Euclidian parameters but also some other infinite-dimensional nuisance parameters. To deal with theoretical challenges due to the additional functional nuisance parameters, we develop a new sieve M-theorem for the asymptotic theory of a general family of semi-parametric M-estimators. 
Moreover, we apply the proposed general theorem to establish the asymptotic normality and semi-parametric efficiency of the proposed sieve MLE $\hat{\beta}_n$ when the convergence rate of the sieve estimator of the nuisance parameter can be slower than $\sqrt{n}$.  We present regularity conditions and main theorems in this section and give all the proofs in the Supplemental Material.

For the simplicity of notation, we focus on model (\ref{model: ltm + cox}) without covariates $Z$, i.e. the linear transformation model (\ref{transformation ode}),  and the results can be similarly extended to the general case with additional regularity conditions on $Z$. Recall that we have set $\gamma(\cdot) = \log \alpha(\cdot)$ and $g(\cdot) = \log q(\cdot)$ to ensure the positivity of $\alpha(\cdot)$ and $q(\cdot)$ in (\ref{transformation ode}). Then we reformulate the ODE model as follows,
\begin{equation}
\label{ode: log double}
\left\{
\begin{array}{lr}
\Lambda'(t) = \exp(x^T\beta + \gamma(t)+g(\Lambda(t))) \\
\Lambda(0) =0
\end{array}
\right..
\end{equation}
Note that the parameter $\beta$ is identifiable when time-varying effects are considered, but in (\ref{ode: log double}) it is identifiable only up to a scaling factor when both $\gamma$ and $g$ are unknown as shown in Proposition~\ref{identi}. To guarantee the identifiability, we constrain the first element of $\beta$ to be $1$ and $\gamma(0)=c$ with some constant $c$ for simplicity in optimization. Let $X=(X_{(1)}, X_{(-1)})$, $\beta=(1, \bar{\beta}^T)^T$, $\bar{\gamma}(\cdot)=\gamma(\cdot)-\gamma(0)$ with $\bar{\gamma}(0)\equiv 0$, and $\bar{X}_{(1)}=X_{(1)}+\gamma(0)$, then we have $X^T\beta + \gamma(t)=\bar{X}_{(1)}+ X_{(-1)}^T \bar{\beta}+ \bar{\gamma}(t)$. We substitute $\bar{\beta}$, $\bar{\gamma}$, and $\bar{X}_{(1)}$ by $\beta$, $\gamma$, and $X_{(1)}$ respectively for notational simplicity hereafter, and the ODE (\ref{ode: log double}) is then equivalent to  
\begin{equation}
\label{ode: log double iden}
\left\{
\begin{array}{lr}
\Lambda'(t) = \exp(x_{(1)}+ x_{(-1)}^T \beta + \gamma(t)+g(\Lambda(t))) \\
\Lambda(0) =0
\end{array}
\right.,
\end{equation}
with $\gamma(0)\equiv 0$. Before stating the regularity conditions, we first introduce some notations. 
We denote the solution of (\ref{ode: log double iden}) by $\Lambda(t, x, \beta, \gamma, g)$ to explicitly indicate that the solution of (\ref{ode: log double iden}) depends on covariates $x$ and parameters $(\beta, \gamma, g)$. We denote the true parameters by $(\beta_0, \gamma_0, g_0)$ and simplify $\Lambda(t, x, \beta_0, \gamma_0, g_0)$ as $\Lambda_0(t,x)$. In addition, some commonly used notations in the empirical process literature will be used in this section as well. Let $Pf = \int f(x) Pr(\textit{d} x)$, where $Pr$ is a probability measure, and denote the empirical probability measure as $\mathbb{P}_n$.  

Then we assume the following regularity conditions. 
\begin{enumerate}[label=(C\arabic*), ref=(C\arabic*)]
	\item\label{c: bounded b} The true parameter $\beta_0$ is an interior point of a compact set $\mc{B}\subset {\mbox{\textbf{R}}}^d$.
	\item\label{c: bounded x} The density of $X$ is bounded below by a constant $c>0$ over its domain $\mc{X}$, which is a compact subset of ${\mbox{\textbf{R}}}^{d+1}$, and $P (X_{(-1)}X_{(-1)}^T)$ is nonsingular.
	\item\label{c: bounded t} There exists a truncation time $\tau <\infty$ such that, for some positive constant $\delta_0$, $Pr(Y>\tau|X)\geq \delta_0$ almost surely with respect to the probability measure of $X$. Then there is a constant $\mu=\sup_{x\in\mc{X}}\Lambda_0(\tau, x) \leq - \log \delta_0 $ such that $\Lambda_0(\tau, X)= -\log Pr(T>\tau|X)\leq \mu$ almost surely with respect to the probability measure of $X$. 
	\item\label{c: smoothness} Let $S^p([a, b])$ be the collection of bounded functions $f$ on $[a, b]$ with bounded derivatives $f^{(j)}$, $j =1, \dots, k$, where the $k$th derivative $f^{(k)}$ satisfies the $m$-H\"older 
	continuity condition: 
	\[|f^{(k)}(s)- f^{(k)}(t)|\leq L|s-t|^m\text{\ \ \ \ \ \ for }s, t\in [a, b],\]
	where $k$ is a positive integer and $m\in (0,1]$ with $p=m+k$, 
	and $L<\infty$ is a constant. The true function $\gamma_0(\cdot)$ belongs to $\Gamma^{p_1}= \{\gamma \in S^{p_1}([0, \tau]):\gamma(0)=0\}$ with $p_1\geq 2$ and the true function $g_0(\cdot)$ belongs to $S^{p_2}([0, \mu+\delta_1]) = \mc{G}^{p_2}$ with some positive constant $\delta_1$ and $p_2 \geq 3$. 
	\item\label{c: lower bound1} Denote $R(t)=\int^t_0\exp(\gamma_0(s))\textit{d} s$, $V=X_{(1)}+X_{(-1)}^T\beta_0$, and $U = e^{V}R(Y)$. There exists $\eta_1\in (0,1)$ such that for all $u\in {\mbox{\textbf{R}}}^d$ with $\|u\|=1$, 
	\[u^T Var(X_{(-1)}\mid U,V)u\geq\eta_1 u^T P(X_{(-1)}X_{(-1)}^T\mid U,V)u \ \ \text{ almost surely.}\]
	\item\label{c: lower bound2} Let $\psi(t, x, \beta, \gamma, g)= x_{(1)}+x_{(-1)}^T\beta + \gamma(t)+g(\Lambda(t, x, \beta, \gamma, g))$ and denote its functional derivatives with respect to $\gamma(\cdot)$ and $g(\cdot)$ along the direction $v(\cdot)$ and $w(\cdot)$ at the true parameter by $\psi'_{0\gamma}(t,x)[v]$ and $\psi'_{0g}(t,x)[w]$ respectively, whose rigorous definitions are given by (S17)-(S18) in the Supplemental Material. For any $v(\cdot)\in \Gamma^{p_1}$ and $w(\cdot)\in \mc{G}^{p_2}$, there exists $ \eta_2 \in (0, 1)$ such that 
	\[(P\{\psi'_{0\gamma}(Y,X)[v]\psi'_{0g}(Y,X)[w]\,|\,\Delta=1\})^2\leq \eta_2 P\{(\psi'_{0\gamma}(Y,X)[v])^2\,|\,\Delta=1\} P\{(\psi'_{0g}(Y,X)[w])^2\,|\,\Delta=1\}\]
	almost surely.
\end{enumerate}

Conditions \ref{c: bounded b}-\ref{c: bounded t} are common regularity assumptions in survival analysis. 
Condition \ref{c: smoothness} requires $p_2\geq 3$ to control the error rates of the spline approximation for the true function $g_0$ and its first and second derivatives. Moreover, together with $p_1\geq 2$, \ref{c: smoothness} will also be used to verify the assumptions (A4)-(A6) for the general M-theorem (Theorem \ref{thm: m-theorem}) when we apply it to derive the asymptotic normality of the proposed sieve MLE (Theorem \ref{thm: asym normal}).  
A similar condition to \ref{c: lower bound1} was  imposed by \citet{wellner2007} for the panel count data, by \citet{ding2011} for the linear transformation model with a known transformation, and by \citet{Zhao_2017} for the accelerated hazards model. When the transformation function is known, condition \ref{c: lower bound1} is equivalent to the assumption C7 in \citet{ding2011} and can be verified in many applications as shown in \citet{wellner2007}. For the general case where both the transformation function and the error distribution are unspecified, condition \ref{c: lower bound2} is assumed to avoid strong collinearity between $\psi'_{0\gamma}(Y,X)[v]$ and $\psi'_{0g}(Y,X)[w]$.

Note that the parameter $g(\cdot)$ takes $\Lambda(t, x, \beta, \gamma, g)$ as its argument in (\ref{ode: log double iden}), which involves the other parameters $\beta$ and $\gamma(\cdot)$. Thus, $\beta$, $\gamma(\cdot)$ and $g(\cdot)$ are bundled parameters. For any $g(\cdot) \in \mc{G}^{p_2}$, we directly consider the composite function $g(\Lambda(t, x, \beta, \gamma, g))$ as a function from $\mc{T}\times \mc{X}\times \mc{B} \times \Gamma^{p_1}$ to ${\mbox{\textbf{R}}}$. And we define the collection of functions
\begin{align*}
	\mc{H}^{p_2}=\{\zeta(\cdot, \beta, \gamma):  \zeta(t, x, \beta, \gamma)=& g(\Lambda(t, x, \beta, \gamma, g)), t\in [0, \tau], x\in \mc{X}, \beta\in \mc{B}, \gamma \in \Gamma^{p_1}, \\
	& g\in \mc{G}^{p_2} \text{ such that } \sup_{t\in[0,\tau],x\in\mc{X}}|\Lambda(t,x,\beta, \gamma, g)| \leq \mu+\delta_1\},
\end{align*}
with $\delta_1$ given in condition \ref{c: smoothness}. 
For any $\zeta(\cdot, \beta, \gamma) \in \mc{H}^{p_2}$, we define its norm as \[\|\zeta(\cdot, \beta, \gamma)\|_2 = \left[ \int_{\mc{X}} \int_{0}^{\tau} [\zeta(t, x, \beta, \gamma)]^2 \textit{d} \Lambda_0(t, x) \textit{d} F_X(x)\right]^{1/2},\] where $F_X(x)$ is the cumulative distribution function of $X$. 
Denote the parameter $\theta = (\beta, \gamma(\cdot), \zeta(\cdot, \beta, \gamma))$ and the true parameter $\theta_0 = (\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0))$ with $\zeta_0(t,x, \beta_0, \gamma_0) =g_0(\Lambda(t,x, \beta_0, \gamma_0, g_0 ))$. 
Denote the parameter space by $\Theta=\mc{B} \times \Gamma^{p_1}\times \mc{H}^{p_2}$. For any $\theta_1$ and $\theta_2$ in $\Theta$, we define the distance \[d(\theta_1, \theta_2)= \left( \|\beta_1-\beta_2\|^2+\|\gamma_1-\gamma_2\|_{2}^2+\|\zeta_1(\cdot, \beta_1, \gamma_1) - \zeta_2(\cdot, \beta_2, \gamma_2)\|_2^2 \right)^{1/2},\] where $\|\cdot\|$ is the Euclidean norm and $\|\gamma\|_{2}=(\int^\tau_0(\gamma(t))^2\textit{d} t)^{1/2}$ is the $L_2$ norm. 

Next, we construct the sieve space as follows. Let $0 = t_0 < t_1 < \cdots < t_{K_n^1}<t_{K_n^1+1}= \tau$ be a partition of $[0, \tau]$ with $K_n^1= O(n^{\nu_1})$ and $\max_{1\leq j\leq K_n^1+1} |t_j - t_{j-1}|=O(n^{-\nu_1})$ for some $\nu_1\in (0,0.5)$. Let $T_{K_n^1}=\{t_1, \cdots, t_{K_n^1}\}$ denote the set of partition points and $S_{n}(T_{K_n^1}, K_n^1, p_1)$ be the space of polynomial splines of order $p_1$ as defined in \citet[page 108, Definition 4.1]{schumaker_2007}. Similarly, let $T_{K_n^2}$ be a set of partition points of $[0,\mu]$ with $K_n^2= O(n^{\nu_2})$ and $\max_{1\leq j\leq K_n^2+1} |t_j - t_{j-1}|=O(n^{-\nu_2})$ for some $\nu_2\in (0,0.5)$, and $S_{n}(T_{K_n^2}, K_n^2, p_2)$ be the space of polynomial splines of order $p_2$. According to \citet[page~117, Corollary~4.10]{schumaker_2007}, there exist two sets of B-spline bases $\{B_j^1, 1\leq j \leq q_n^1\}$ with $q_n^1 = K_n^1+p_1$ and $\{B_j^2, 1\leq j \leq q_n^2\}$ with $q_n^2 = K_n^2+p_2$ such that for any $s_1 \in S_{n}(T_{K_n^1}, K_n^1, p_1)$ and $s_2 \in S_{n}(T_{K_n^2}, K_n^2, p_2)$, we can write $s_1(t) = \sum_{j=1}^{q_n^1}a_j B_j^1(t)$ and $s_2(t) = \sum_{j=1}^{q_n^2}b_j B_j^2(t)$.  Let $\Gamma^{p_1}_n = \{\gamma \in S_{n}(T_{K_{n}^1}, K_{n}^1, p_1): \gamma(0)=0\}$, $\mc{G}^{p_2}_n =S_{n}(T_{K_{n}^2}, K_{n}^2, p_2)$, and \[\mc{H}^{p_2}_n=\{\zeta(\cdot, \beta, \gamma): \zeta(t, x, \beta, \gamma)=g(\Lambda(t, x, \beta, \gamma, g)), g\in \mc{G}^{p_2}_n, t\in [0, \tau], x\in \mc{X}, \beta\in \mc{B}, \gamma \in \Gamma^{p_1}_n\}.\] Let $\Theta_n=\mc{B} \times \Gamma^{p_1}_n\times \mc{H}^{p_2}_n$ be the sieve space. It is not difficult to see that $\Theta_n \subset \Theta_{n+1}\subset \cdots \subset \Theta$. We consider the sieve estimator $\hat{\theta}_n = (\hat{\beta}_n, \hat{\gamma}_n(\cdot), \hat{\zeta}_n(\cdot, \hat{\beta}_n, \hat{\gamma}_n))$, where $ \hat{\zeta}_n(t, x, \hat{\beta}_n, \hat{\gamma}_n) = \hat{g}_n(\Lambda(t, x, \hat{\beta}_n, \hat{\gamma}_n, \hat{g}_n))$, that maximizes the log-likelihood (\ref{obj}) (without covariates $Z$ and parameter $\boldsymbol{\eta}$) over the sieve space $\Theta_n$. The consistency and convergence rate of the sieve MLE $\hat{\theta}_n$ are then established in the following theorem.

\begin{theorem}(Convergence rate of $\hat{\theta}_n$.)
	\label{thm: conv rate}
	Let $\nu_1$ and $\nu_2$ satisfy the restrictions $\max\{\frac{1}{2(2+p_1)},\frac{1}{2p_1}-\frac{\nu_2}{p_1} \}<\nu_1<\frac{1}{2p_1}$, $\max\{\frac{1}{2(1+p_2)},\frac{1}{2(p_2-1)}-\frac{2\nu_1}{p_2-1}\}<\nu_2<\frac{1}{2p_2}$,  and $2\min\{2\nu_1,\nu_2\} > \max\{\nu_1,\nu_2\}$. Suppose conditions \ref{c: bounded b}-\ref{c: lower bound2} hold, then we have
	\[d(\hat{\theta}_n, \theta_0)=O_p(n^{-\min\{p_1\nu_1, p_2\nu_2, \frac{1-\max\{\nu_1, \nu_2\}}{2}\}}).\] 
\end{theorem}

Theorem \ref{thm: conv rate} gives the convergence rate of the proposed estimator $\hat{\theta}_n$ to the true parameter $\theta_0$, and its proof is provided in the Supplemental Material 
by verifying the conditions in~\citet[Theorem 1]{shen1994}. Note the subscripts 1 and 2 correspond to the space of the spline approximation for two infinite-dimensional parameters $\gamma$ and $g$, respectively. 
The restrictions on $\nu_1$ and $\nu_2$ are feasible for $p_1$ and $p_2$ not far away from each other. For example, if $p_1=p_2=p$ and $\nu_1=\nu_2=\nu$, the restriction on~$\nu$ is equivalent to $\frac{1}{2(1+p)}< v<\frac{1}{2p}$, and the convergence rate becomes  $d(\hat{\theta}_n, \theta_0) = O_p(n^{-\min\{p\nu, \frac{1-\nu}{2}\}})$, which is the same as the case when there is only one infinite-dimensional parameter~\citep{ding2011, Zhao_2017}. Further, if $\nu=\frac{1}{1+2p}$, we have $d(\hat{\theta}_n, \theta_0) = O_p(n^{-\frac{p}{1+2p}})$, which achieves the optimal convergence rate in the nonparametric regression setting. 

Although the convergence rate for the nuisance parameter is slower than the typical rate $n^{1/2}$, we will show that the sieve MLE of the regression parameter, i.e. $\hat{\beta}_n$, is still asymptotically normal and achieves the semi-parametric efficiency~bound.  First, we introduce two additional regularity conditions which are stated below.

\begin{enumerate}[label=(C\arabic*), ref=(C\arabic*),resume]
	\item \label{c: pos info exist}
	There exist $\mb{v}^*=(v^*_1, \cdots, v^*_d)^T$ and $\mb{w}^*=(w^*_1, \cdots, w^*_d)^T$, where $v^*_j \in  \Gamma^{2}$ and $w^*_j \in \mc{G}^{2}$ for $j=1, \cdots, d$, such that $P\{\Delta \mb{A}^*(U, X)\psi'_{0\gamma}(Y,X)[v]\}=0$ and $P\{\Delta \mb{A}^*(U, X)\psi'_{0g}(Y,X)[w]\}=0$ 
	hold for any $v\in \Gamma^{p_1}$ and $w\in \mc{G}^{p_2}$. Here $U$ and $V$ are defined the same as in condition \ref{c: lower bound1} and
	\begin{align*}
		\mb{A}^*(t, X)=& -\left(g_0'(\tilde{\Lambda}_{0}(t))\exp(g_0(\tilde{\Lambda}_{0}(t)))t +1\right)X_{(-1)}\\
		& + g_0'(\tilde{\Lambda}_{0}(t))\exp(g_0(\tilde{\Lambda}_{0}(t))) \int^{t}_0 \mb{v}^*(R^{-1}(se^{-V}))\textit{d} s + \mb{v}^*(R^{-1}(te^{-V}))\\
		& + g_0'(\tilde{\Lambda}_{0}(t))\exp(g_0(\tilde{\Lambda}_{0}(t)))\int^{\tilde{\Lambda}_{0}(t)}_0 \exp(-g_0(s))\mb{w}^*(s)\textit{d} s + \mb{w}^*(\tilde{\Lambda}_{0}(t)),
	\end{align*}
	where $\tilde{\Lambda}_0(t)$ is the solution of $\tilde{\Lambda}_0'(t)=\exp(g_0(\tilde{\Lambda}_0))$ with $\tilde{\Lambda}_0(0)=0$.
	\item\label{c: pos fish} Let $\boldsymbol{l}^*(\beta_0, \gamma_0, \zeta_0;W) = \int \mb{A}^*(t, x)\textit{d} M(t)$, where 
$M(t)=\Delta \mathbbm{1}(U \leq t) - \int_{-\infty}^t \mathbbm{1}(U \geq s)\textit{d} \tilde{\Lambda}_0(s)$ is the event counting process martingale. The information matrix $I(\beta_0) = P(\boldsymbol{l}^*(\beta_0, \gamma_0, \zeta_0;W)^{\otimes 2})$ is nonsingular. Here for a vector $a$, $a^{\otimes 2} = aa^T$.
\end{enumerate}

The additional condition \ref{c: pos info exist} essentially requires the existence of the least favorable direction that is used to establish the semi-parametric efficiency bound. The directions $\mb{v}^*$ and $\mb{w}^*$ may be found through the equations in \ref{c: pos info exist}. We illustrate how to construct $\mb{v}^*$ and $\mb{w}^*$ for the Cox model and   the linear transformation model with a known transformation respectively in Remark~\ref{rmk: least favorable direction}. Condition \ref{c: pos fish} is a natural assumption that requires the information matrix to be invertible. The following theorem establishes the asymptotic normality and semi-parametric efficiency of the sieve MLE $\hat{\beta}_n$ of the regression parameter for the general linear transformation model.

\begin{theorem}(Asymptotic normality of $\hat{\beta}_n$)
	\label{thm: asym normal}	 Suppose the conditions in Theorem \ref{thm: conv rate} and  \ref{c: pos info exist}-\ref{c: pos fish} hold, then we have
\begin{align*}
		\sqrt{n}(\hat{\beta}_n - \beta_0) & = \sqrt{n}I^{-1}(\beta_0)\mathbb{P}_n \boldsymbol{l}^*(\beta_0, \gamma_0, \zeta_0;W) + o_p(1) \rightarrow_d N(0, I^{-1}(\beta_0))
	\end{align*}
	with $I(\beta_0)$ given in condition \ref{c: pos fish} and $\rightarrow_d$ denoting convergence in distribution.
\end{theorem}

Theorem \ref{thm: asym normal} states that $\hat{\beta}_n$ is asymptotically normal with variance as the inverse of the information matrix. 
In practice, the information matrix can be approximated by the estimated information matrix of all parameters including the coefficients of spline bases. 

We note that the existing sieve M-theorem for bundled parameters \citep{ding2011,Zhao_2017} cannot be directly applied to prove Theorem \ref{thm: asym normal}. Because it does not allow the infinite-dimensional nuance parameter 
to be a function of another infinite-dimensional nuance parameter. 
Therefore, to study the asymptotic distribution of $\hat{\beta}_n$, we first establish a new general M-theorem for bundled parameters where the infinite-dimensional nuisance parameter is a function of not only the Euclidean parameter of interest but also another infinite-dimensional nuisance parameter. The established M-theorem under such a general scenario then enables us to prove Theorem \ref{thm: asym normal} by verifying its assumptions for the linear transformation model.  The details are provided in the Supplemental Material.  Since the new M-theorem can be useful for developing the asymptotic normality of sieve estimators for other ODE models, we state it below for readers of interest.

We first introduce the general setting and notation for the proposed sieve M-theorem. 
Let $m(\theta; W)$ be an objective function of unknown parameters $\theta = (\beta, \gamma(\cdot), \zeta(\cdot, \beta, \gamma))$ given a single observation $W$, where $\beta$ is the finite-dimensional parameter of interest, $\gamma(\cdot)$ is an infinite-dimensional nuisance parameter, and $\zeta(\cdot,\beta,\gamma)$ is another infinite-dimensional nuisance parameter that can be a function of $\beta$ and $\gamma$. Here ``$\cdot$'' represents some components of $W$. Given i.i.d. observations $\{W_i\}^{n}_{i=1}$, the sieve estimator $\hat{\theta}_n = (\hat{\beta}_n, \hat{\gamma}_n(\cdot), \hat{\zeta}_n(\cdot, \hat{\beta}_n, \hat{\gamma}_n))$ maximizes the objective function, $\mathbb{P}_n m(\theta; W)$, over certain sieve space. For example, $\hat{\theta}_n$ becomes the sieve MLE if $m$ is the log-likelihood function. We denote the derivative of $m$ with respect to $\beta$ as $m'_{\beta}$, the functional derivative of $m$ with respect to $\gamma$ along the direction $v(\cdot)$ as $m'_{\gamma}[v]$, and the functional derivative of $m$ with respect to $\zeta$ along the direction $h(\cdot)$ as $m'_{\zeta}[h]$, whose rigorous definitions are given in the Supplemental Material. The following theorem then establishes the asymptotic normality of the sieve estimator, $\hat{\beta}_n$, under the above general setting.

\begin{theorem}(A general M-theorem for bundled parameters.)
\label{thm: m-theorem}
Under assumptions (A1)-(A6) in the Supplemental Material,  we have
	\begin{align*}
		\sqrt{n}(\hat{\beta}_n - \beta_0) & = A^{-1}\sqrt{n}\mathbb{P}_n \mb{m}^*(\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0);W) + o_p(1) \\
		& \rightarrow_d N(0, A^{-1}B( A^{-1})^T),
	\end{align*}
	where 
	\begin{align*}
		\mb{m}^*(\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0);W) & = m'_{\beta}(\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0);W) -m'_{\gamma}(\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0);W)[\mb{v}^*] \\
		& \ \ \ - m'_{\zeta}(\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0);W)[\mb{h}^*(\cdot, \beta_0, \gamma_0)],\\
		B & = P\{\mb{m}^*(\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0);W)\mb{m}^*(\beta_0, \gamma_0(\cdot), \zeta_0(\cdot, \beta_0, \gamma_0);W)^T\},
	\end{align*}
with $\mb{v}^*=(v^*_1, \dots, v^*_{d})^T$, $\mb{h}^*=(h^*_1,\dots, h^*_{d})^T$ and $A$ given in the assumption (A3). 
\end{theorem}


\begin{remark}
The assumptions needed in Theorem \ref{thm: m-theorem} are similar to those in \citet{ding2011} (see the Supplemental Material for details).  However, our proposed theorem significantly differs from the main theorem in \citet{ding2011}, because the latter considers $\zeta(\cdot, \beta)$ to be a function of only the finite-dimensional parameter $\beta$, while we consider a more general scenario of bundled parameters, where the nuisance parameter $\zeta(\cdot, \beta,\gamma)$ can be a function of both the finite-dimensional parameter $\beta$ and another infinite-dimensional nuisance parameter $\gamma$. The proposed theorem nontrivially extends the asymptotic distributional theories for M-estimation under this general scenario. 
\end{remark}
 
 
\begin{remark}
\label{rmk: least favorable direction}
Finally, we note that to find the least favorable directions $\mb{v}^*$ and $\mb{w}^*$ required in \ref{c: pos info exist}, we may solve the equations in \ref{c: pos info exist}, which can be simplified to equations (S35) and (S37) provided in the Supplemental Material. For illustration, we provide explicit constructions of the least favorable directions for the Cox model and for the linear transformation model with a known transformation respectively. Specifically, for the Cox model, we have $g_0\equiv 0$ and $\mb{v}^*$ can be derived as 
\[\mb{v}^*(t) = \frac{P\{\mathbbm{1}(Y\geq t)e^{X_{(1)}+X_{(-1)}^T\beta_0}X_{(-1)}\}}{P\{\mathbbm{1}(Y\geq t)e^{X_{(1)}+X_{(-1)}^T\beta_0}\}};\]
for the linear transformation model where $\gamma_0$ is known, $\mb{w}^*$ can be obtained as  
\[\mb{w}^*(t) = \pmb{\phi}(t) - g_0'(t)\int^{t}_0\pmb{\phi}(s)\textit{d} s,\]
where 
\[  \pmb{\phi}(t) = \left(g_0'(t)\exp(g_0(t))\tilde{\Lambda}_0^{-1}(t) +1\right) \frac{P\{\mathbbm{1}(\Lambda_0(Y,X)\geq t)X_{(-1)}\}}{P\{\mathbbm{1}(\Lambda_0(Y,X)\geq t)\}}\]
with $\tilde{\Lambda}_0$ defined in \ref{c: pos info exist}.
\end{remark}
\section{Conclusion}
\label{s:discuss}
In this paper, we have proposed a novel ODE framework for survival analysis, which unifies the current literature, along with a general estimation procedure which is scalable and easy to implement. The ODE framework provides a new perspective for modeling censored data, which further allows us to utilize well-developed numerical solvers and local sensitivity analysis tools for ODEs in parameter estimation. Although we have only focused on one class of ODE models in this paper, the ODE framework and the estimation method offer new opportunities for investigating more flexible model structures.

We have also developed a general theory to derive the asymptotic distribution of bundled parameters, where the nuisance parameter is a function of not only the regression parameters of interest but also other infinite-dimensional nuisance parameters. Though we have only illustrated the efficient estimation in the linear transformation model as an example to motivate such a theoretical development, the proposed general theory can be extended to other models.


\subsection{Cox model and time-varying Cox model}
\label{subs: cox}
The Cox proportional hazard model assumes that the covariates have a multiplicative effect on the hazard function, i.e. $\lambda_{x}(t)=\alpha(t)\exp(x^T \beta)$,  
where $\alpha(t)$ is a baseline hazard function and $\exp(x^T \beta)$ is the relative risk, and extensions of the Cox model allow for time-varying coefficients~\citep{Zucker_1990, Gray_1994}. Here we write the Cox model with both time-independent and time-varying effects as a simple ODE, whose right-hand side does not depend on the cumulative function, i.e.
\begin{equation}
\label{cox ode}
	 \Lambda'_{x, z}(t) = \alpha(t) \exp(x^T \beta + z^T \boldsymbol{\eta}(t)),
\end{equation}
which allows covariates $x$ to have time-independent effects and covariates $z$ to have time-varying effects on the hazard function. The baseline hazard function $\alpha(t)$ and time-varying effects $\boldsymbol{\eta}(t)$ can be specified in a parametric model or left unspecified in a semi-parametric model.



\subsection{Accelerated failure time model}
\label{subs: aft}
The AFT model assumes that the log transformation of $T$ is linearly correlated with covariates, i.e. $\log T = - X^T \beta + \epsilon$. 
In the proposed ODE framework, the AFT model can be written as 
\begin{equation}
\label{eq: aft ode}
	 \Lambda'_{x}(t) = q(\Lambda_{x}(t)) \exp(x^T \beta),
\end{equation}
where the function $q(\cdot)$ uniquely determines the distribution of error $\epsilon$ in the following way. Let $H_q(u) = \int_0^{-\ln u} {q^{-1}(v)}\textit{d} v$ and $ G_{q}(u) = H_q^{-1}(u)$, then $G_{q}$ is the survival function of $\delta =\exp(\epsilon)$ as shown in \citet{bagdonavicius2001accelerated}. For example, if $q(t)=v k^{\frac{1}{v}} t^{1-\frac{1}{v}}$, then $\delta$ follows a Weibull distribution with $G_q(t)=\exp(-kt^v)$. When the error distribution is unknown (as in a semi-parametric AFT model), we can leave the function $q(\cdot)$ unspecified.

The ODE (\ref{eq: aft ode}) provides a new and clear interpretation on how covariates affect the hazard for the AFT model. Specifically, it implies that given an individual's own cumulative hazard, covariates $x$ have a multiplicative constant effect on the hazard function. Further, besides the direct effects of covariates, if $q(\cdot)$ is a monotonic increasing function, then an individual with a higher cumulative hazard at a particular time would have a higher ``baseline" hazard. Note that although we can also present the hazard directly as a function of covariates and time, i.e. $\lambda_x(t) = \lambda_{\delta}(t \exp(x^T\beta))\exp(x^T\beta)$, the covariate effects are entangled with the baseline hazard $\lambda_\delta$ in this representation, which is more difficult to interpret. 

\subsection{Linear transformation model}
\label{subs: ltm}
As an extension of the AFT model, the linear transformation model assumes that, after a monotonic increasing transformation $\varphi(\cdot)$, the event time $T$ is linearly correlated with covariates, i.e. $\varphi(T) = -X^T\beta+\epsilon$. In the proposed ODE framework, it can be written as 
\begin{equation}
    \label{transformation ode}
   \Lambda'_{x}(t)  = q(\Lambda_{x}(t)) \exp(x^T \beta) \alpha(t),
\end{equation}
where $q(\cdot)$ corresponds to the distribution of $\epsilon$ in the same way as in the AFT model, and $\alpha(\cdot)$ is uniquely determined by the equation $\varphi(t) = \log \int_0^t \alpha(s)\textit{d} s$.
In comparison to model (\ref{eq: aft ode}), the hazard function at time $t$ depends not only on the current cumulative hazard and covariates, but also on the current time $t$ directly.


Different specifications of $\varphi(\cdot)$ and $\epsilon$ have been proposed in the literature for the linear transformation model. We consider the case where both the transformation and the error distribution are unknown. This specification is especially preferred when parametric assumptions on the transformation function or the error distribution cannot be properly justified. However, when both $q(\cdot)$ and $\alpha(\cdot)$ are unknown, they may not be identifiable. The equivalent linear regression representation, $\varphi(T) = - x^T \beta+\epsilon$, allows us to see the identifiability issue clearly. For example, for any constants $c_1 > 0$ and $c_2 \neq 0$, we have $c_1\varphi(T)+c_2 = -x^T (c_1\beta)+(c_1\epsilon + c_2)$. Therefore, $(\beta, \varphi, \epsilon)$ and $(c_1 \beta, c_1 \varphi+c_2 , c_1 \epsilon+c_2)$ are not distinguishable. Under certain regularity conditions, \citet{horowitz1996semiparametric} showed that the model parameters are identifiable up to a scale and a location normalization. Following that result, we have developed Proposition \ref{identi} that characterizes the identifiability of parameters in (\ref{transformation ode}). Proposition \ref{degeneration} provides necessary and sufficient degeneration conditions for AFT and Cox models.  The proofs are given in the Supplemental Material. 

\begin{proposition}
\label{identi}
Suppose the covariates in $x$ are continuous and at least one of the components has a non-zero $\beta$ coefficient, which without loss of generality is assumed to be positive. Let $(q(\cdot), \beta, \alpha(\cdot))$ specify the survival distribution through (\ref{transformation ode}). Then for any other $(\tilde{q}(\cdot), \tilde{\beta}, \tilde{\alpha}(\cdot))$ that gives the same survival distribution, if and only if there exists positive constants $c_1$ and $c_2$ such that $\tilde{\beta} = c_1\beta,$ $\int_{0}^t\tilde{\alpha}(s)\textit{d} s= c_2 (\int_{0}^t \alpha(s)\textit{d} s)^{c_1}$,
and $\int_{0}^t {\tilde{q}^{-1}(s)}\textit{d} s =  c_2 (\int_{0}^t {q^{-1}(s)}\textit{d} s)^{c_1}$
for any $t$.
\end{proposition}

\begin{proposition}
\label{degeneration}
Suppose the conditions in Proposition~\ref{identi} hold, then the linear transformation model in~(\ref{transformation ode}) coincides with the Cox model if and only if there exist positive constants $c_1$ and $c_2$ such that $q(u) = c_2u^{1-c_1}$, 
and it coincides with the AFT model if and only if there exist positive constants $c_1$ and $c_2$ such that $\alpha(t)=c_2 t^{c_1-1}$. 
\end{proposition}


\begin{remark}
\label{rm: other examples}
The proposed ODE framework is general enough to cover other existing models as well. For example, both the additive hazard model  \citep{aalen1980model, MCKEAGUE_1994} and the additive-multiplicative hazard model \citep{lin1995semiparametric} can be viewed as a specific ODE model, i.e. $\Lambda'_{x, z}(t)=r_1 (x^T \beta) + \alpha (t) r_2(z^T \eta)$, where $r_1(\cdot)$ and $r_2(\cdot)$ are some known link functions. Subsequently, the generalized additive hazards model and the generalized additive-multiplicative hazards model \citep{bagdonavicius2001accelerated} can be written as $\Lambda'_{x}(t)=q(\Lambda_{x}(t))(r_1 (x) + \alpha (t) r_2(x))$. The generalized Sedyakin's model \citep{bagdonavicius2001accelerated}, which was proposed as an extension of the AFT model, can also be viewed as a special case of (\ref{general ode}) with $\Lambda'_{x}(t) = f(\Lambda_{x}(t), x)$.
\end{remark}

\begin{remark}
Further, the proposed ODE framework and the estimation method in Section \ref{s: mle} can also be extended to deal with time-varying covariates. Suppose the covariate is a stochastic process $X(t), t\geq 0$ and $T_{X(\cdot)}$ is the failure time under $X(\cdot)$. Denote the conditional survival, the hazard function, and the cumulative function by
$S_{x(\cdot)}(t) = P\{T_{X(\cdot)}\geq t|X(s) = x(s), 0\leq s\leq t\}$,
$\lambda_{x(\cdot)}(t) =- \frac{S_{x(\cdot)}'(t)}{S_{x(\cdot)}(t)}$, and $\Lambda_{x(\cdot)}(t) =-\log(S_{x(\cdot)}(t))$, respectively.  Then the ODE (\ref{general ode}) can be extended to $\Lambda'_{x(\cdot)}(t) = f(t, \Lambda_{x(\cdot)}(t), x(t))$.
For presentation simplicity, we focus on models in the form of (\ref{model: ltm + cox}) in this paper.
\end{remark}


\subsection{Related estimation methods and their limitations}
\label{s:estimation}

The maximum partial likelihood estimator (MPLE) \citep{cox1975partial} was first proposed for the Cox model, and the asymptotic property of MPLE was established by \citet{Andersen_1982} via the counting process martingale theory. For time-varying Cox model, many different estimation methods have been developed while relying on maximizing the partial likelihood \citep{Zucker_1990, Gray_1994}. 
However, evaluating the partial likelihood for an uncensored individual requires access to all other observations who were in its risk set. This prevents parallel computing for partial likelihood-based methods, which is a drawback when analyzing large scale data.

For the linear transformation model, different specifications of the transformation and the error distribution along with different estimation methods have been proposed. For example, \cite{10.1093/biomet/82.4.835}, \cite{10.1093/biomet/85.4.980}, \cite{chen2002semiparametric}, and \cite{bagdonavicius1999generalized} have considered an unknown transformation with a known error distribution, which includes the Cox model and the proportional odds model \citep{doi:10.1002/sim.4780020223} as special cases.  The corresponding
modified MPLE and NPMLE have also been developed \citep{chen2002semiparametric, bagdonavicius1999generalized, murphy1997maximum, zeng2007maximum}. 
However, due to the large number of nuisance parameters, it is difficult to obtain NPMLE in practice, especially in large-scale applications.  Alternatively, \citet{10.1093/biomet/92.3.619} considered a parametric Box-Cox transformation with an unknown error distribution, which includes the semi-parametric AFT model as a special case, and least square and rank-based methods have been proposed to estimate the regression parameters \citep{10.2307/2335161, lai1991large, tsiatis1990estimating, jin2003rank, jin2006least}. 
Nevertheless, they are not asymptotically efficient and may suffer additional numerical errors resulted from discrete objective functions. Subsequently, under the AFT model, \citet{zeng2007efficient} and \citet{Lin_2012} proposed efficient estimators based on a kernel-smoothed profile likelihood, and \citet{ding2011} developed an efficient sieve MLE. 
When both the transformation function and the error distribution are unknown, a partial rank-based method has been proposed  \citep{khan2007partial, Song_2006}, and its computation is analogous to that of the partial likelihood, where the rank of an uncensored individual is determined by all other individuals in its risk set, and thus the computational challenge for large-scale applications still remains.

As evident from the above discussion, many existing estimation methods suffer from important limitations in practice. In Section \ref{s: mle}, we propose a scalable, easy-to-implement and efficient estimation method that can be applied to a wide range of models. 

\section{The ODE Framework}
\label{s: ode framework}
To characterize the conditional distribution of $T$ given $X$, the conditional hazard function, denoted as $\lambda_x(t)=\Lambda'_x(t)$, 
provides a popular modeling target as it describes the instantaneous rate at which the event occurs given survival. In this paper, we view the hazard function as the dynamic change of the cumulative hazard function and quantify them using an ODE.

In our ODE framework, the hazard function depends not only on the time and covariates but also on the cumulative hazard as shown in (\ref{general ode}), where function $f(\cdot)$ specifies the dynamic change of $\Lambda_{x}(t)$ and covariates $x$ serve as additional parameters in terms of the ODE. 
The initial value in (\ref{general ode}) implies that, for an individual with covariates $x$, the probability for an event to occur after $t_0$ is controlled by $\exp(-c(x))$. For example, it is often the case that time $0$ is defined prior to the occurrence of events, which implies that an event always occurs after time $0$, i.e. the survival function $S_x(0)=1$, and it follows that $\Lambda_x(0)=0$. We use this initial value in the ODE framework hereafter for simplicity, while the estimation method and the theoretical properties established later can be extended to the general case where $c(x)$ can be a function of covariates.
Under certain smoothness conditions~\citep[page 108]{Walter1998}, the initial value problem~(\ref{general ode}) has exactly one solution, which uniquely characterizes the conditional distribution of the event time.

Next, we present a general class of ODE models (other examples beyond this class are included in Remark~\ref{rm: other examples} to illustrate the flexibility of the proposed ODE framework). 
Suppose there are two groups of covariates denoted by $X \in {\mbox{\textbf{R}}}^{d_1}$ and $Z\in {\mbox{\textbf{R}}}^{d_2}$ respectively. 
We consider ODE models in the form of  
\begin{equation}
\label{model: ltm + cox}
    \Lambda'_{x, z}(t) = \alpha(t)\exp(x^T\beta +z^T \boldsymbol{\eta}(t))q(\Lambda_{x, z}(t)),
\end{equation}
where $\alpha(\cdot)$ and $q(\cdot)$ are two unknown positive functions, and 
given an individual's own cumulative hazard, both covariates $x$ and $z$ have multiplicative effects on the hazard, one with time-independent coefficients $\beta \in {\mbox{\textbf{R}}}^{d_1}$ and the other with time-varying coefficients $\boldsymbol{\eta}(t) \in {\mbox{\textbf{R}}}^{d_2}$. Here $\boldsymbol{\eta}(\cdot)=(\eta_1(\cdot), \dots, \eta_{d_2}(\cdot))^T. $\footnote{Through out this paper, we bold vectors only when each element is a function.} 
This general class covers many existing models as special cases. In particular, as shown below, model~(\ref{model: ltm + cox}) reduces to the time-varying Cox model when $q(\cdot)=1$, to the linear transformation model when covariates $z$ are not considered, 
and further reduces to the AFT model if $\alpha(\cdot)=1$. 
In the following subsections, we will also show that by rewriting many existing models under the format (\ref{general ode}), the ODE framework brings them new interpretations in terms of the hazard function.

\section{Introduction}
\label{s:intro}
Survival analysis is an important branch of statistical learning, where the primary outcome of interest is the time to a certain event. In practice, event times may not be observed due to a limited observation time window or missing follow-up during the study, which is referred to as censored data. Many statistical models have been developed to deal with censored data in the literature. For example, the Cox proportional hazard model is probably the most classical semi-parametric model for handling censored data~\citep{cox1975partial}, and it assumes that the covariates have a constant multiplicative effect on the hazard function. Although easy to interpret, the constant hazard ratio assumption is often considered as overly strong for real-world applications. As a result, many other semi-parametric models have been proposed as attractive alternatives, such as accelerated failure time (AFT) models, transformation models, and additive hazards models.  See \citet{aalen1980model}, \citet{10.2307/2335161}, \citet{Gray_1994}, \citet{doi:10.1002/sim.4780020223}, \citet{10.1093/biomet/82.4.835}, \citet{10.1093/biomet/85.4.980}, and \citet{chen2002semiparametric} for a sample of references. Given different assumptions made in these semi-parametric models, different estimation and inference procedures have also been developed accordingly, such as partial likelihood based methods~\citep{Zucker_1990, Gray_1994, bagdonavicius2001accelerated, chen2002semiparametric}, least square and rank-based methods~\citep{10.2307/2335161, lai1991large, tsiatis1990estimating, jin2003rank, jin2006least}, 
non-parametric maximum likelihood estimators (NPMLE)~\citep{murphy1997maximum, zeng2007maximum}, and sieve maximum likelihood estimators (MLE)~\citep{huang1999efficient, shen1994, ding2011, Zhao_2017}. 

Despite successful applications of these existing models and their estimation methods, the growing scale and complexity of modern data have placed new demands that are not entirely met. On one hand, more flexible models are often desired to utilize the rich information in data. On the other hand, most existing estimation methods suffer from at least one of the following limitations\footnote{We refer the readers to Section~\ref{s:estimation} for a detailed discussion.}: 1) not allowing parallel computing for large-scale data; 2) complicated numerical implementations; 3) not achieving semi-parametric efficiency bound. 
Therefore, there is an urgent need for flexible survival models with scalable estimation methods.

To address the aforementioned challenges, we introduce a novel Ordinary Differential Equation (ODE) framework for survival analysis. Specifically, the proposed framework models the dynamic change of the cumulative hazard function through an ODE. Let  $T$ be the event time 
and $X$ be covariates. Denote the conditional cumulative hazard function of $T$ given $X=x$ as $\Lambda_{x}(t)$. Then $\Lambda_x(t)$ is characterized by the following ODE with a fixed initial value
\begin{align}
\label{general ode}
\left\{
\begin{array}{lr}
\Lambda'_{x}(t) = f(t, \Lambda_{x}(t), x) \\
\Lambda_{x}(t_0) = c(x)
\end{array}
\right.,
\end{align}
where the derivative is with respect to $t$, $f(\cdot)$ and $c(\cdot)$ are functions to be specified, and $t_0$ is a predefined initial time point. In particular, function $c(\cdot)$ determines the probability of an event occurring after $t_0$; for instance, $\Lambda_x(0)=0$ corresponds to the case when no event occurs before time $0$. And function $f(\cdot)$ determines how covariates $x$ affect the hazard function at time $t$ given an individual's own cumulative hazard. 

The proposed framework is general enough to unify many aforementioned existing survival models through different specifications of the function $f(\cdot)$. For example, the ODE~(\ref{general ode}) is equivalent to the Cox model when $f(\cdot)$ takes the form $\alpha(t)\exp(x^T\beta)$ for some function $\alpha(\cdot)$, and it is equivalent to the AFT model when $f(\cdot)$ takes the form $q(\Lambda_x(t))\exp(x^T\beta)$ for some function $q(\cdot)$. Similarly, we can obtain many more models such as the time-varying variants of the Cox model, the linear transformation model, and the additive hazards model to name a few (see Section \ref{s: examples} for details). Meanwhile, re-writing the existing models in the ODE form provides new interpretations for these models.  In addition, the proposed new perspective of modeling the hazard offers an opportunity for designing more flexible model structures. 


This unified framework also enables the development of a scalable and easy-to-implement procedure for estimation and inference, which can be applied to a wide range of models. In this paper, we illustrate the estimation procedure by using a general class of ODE models as an example. Specifically, this general class includes 
the most flexible linear transformation model, where both the transformation function and the error distribution are unspecified. Since the $f(\cdot)$ function for the general model contains both finite-dimensional and infinite-dimensional parameters, we propose a spline-based sieve MLE that directly maximizes the likelihood in a sieve space. 
Under the ODE framework, the proposed estimation method allows parallel computing for independent observations, which is especially crucial for large-scale applications. Moreover, to address the optimization challenge due to the ODE model structure, we provide an easy-to-implement gradient-based optimization algorithm. In particular, we utilize local sensitivity analysis~\citep{Dickinson_1976} to derive the gradient and evaluate both the likelihood function and its gradient using numerical ODE solvers. Efficient implementations of both ODE solvers and splines are available in many software. As a result, the proposed estimation procedure is easy to carry out in practice.

Further, we establish a general sieve M-theorem for bundled parameters to analyze the theoretical properties of the proposed sieve estimator. While many asymptotic distributional theories for M-estimation in semi-parametric models have been developed (see \citet{huang1999efficient}, \citet{Shen_1997}, \citet{ai2003efficient}, \citet{wellner2007}, \citet{ZHANG_2010}, \citet{He_2010}, \citet{ding2011} for a sample of references), they cannot be directly applied to our setting. Among them, the proposed theory in \citet{ding2011} considers bundled parameters where the infinite-dimensional parameter is an unknown function of the finite-dimensional Euclidean parameter and has been applied to the AFT model, and recently, to the accelerated hazards model in \citet{Zhao_2017}. However, for the general class of ODE models, the estimation criterion is parameterized with more general bundled parameters where the nuisance parameter is an unknown function of not only finite-dimensional regression parameters of interest but also other infinite-dimensional nuisance parameters. To accommodate this different and challenging scenario, we develop a new general sieve M-theorem for bundled parameters. By applying it to the general class of ODE models along with ODE related methodologies, we show consistency, asymptotic normality, and semi-parametric efficiency for the estimated regression parameters. The proposed theory can also be extended to develop the asymptotic normality of estimators for other ODE models.

The rest of the paper is organized as follows. We introduce the ODE framework and present a general class of ODE models as special cases in Section \ref{s: ode framework}. We provide the estimation procedure in Section \ref{s: mle} and establish theoretical properties in Section \ref{s: asymptotics}. Simulation studies and a real-world data example are presented in Sections \ref{s: simulation} and \ref{s: real data} respectively.



\section{Maximum Likelihood Estimation}
\label{s: mle}

In this section, we propose a general estimation procedure that can be applied to a wide range of ODE models. Here we use the ODE model in~(\ref{model: ltm + cox}) as an illustrative example, and the proposed estimation method can also be applied to other models such as those mentioned in Remark \ref{rm: other examples}. 

We denote the event time as $T$, the censoring time as $C$. Let $Y = \min\{T, C\}$ and $\Delta=\mathbbm{1}(T\leq C)$, where $\mathbbm{1}(\cdot)$ denotes the indicator function. Our  data consist of $n$  independent and identically distributed observations $\{Y_i, \Delta_i, X_i, Z_i\},~ i = 1, \dots, n$.  Since $\alpha(\cdot)$ and $q(\cdot)$ in~(\ref{model: ltm + cox}) are positive, we set $\gamma(\cdot) = \log \alpha(\cdot)$ and $g(\cdot) = \log q(\cdot)$. Under the conditional independence between $T$ and $C$ given covariates~$(X, Z)$, the log-likelihood function of the parameters $(\beta, \gamma(\cdot), \boldsymbol{\eta}(\cdot), g(\cdot))$ is given by
\begin{align}
\label{obj}
	l_n(\beta, \gamma(\cdot), g(\cdot), \boldsymbol{\eta}(\cdot))=\frac{1}{n}\sum_{i=1}^n [ & \Delta_i\{ \gamma(Y_i)+ X_i^T\beta +Z_i^T\boldsymbol{\eta}(Y_i)+g(\Lambda_{i}(Y_i;\beta, \gamma, g, \boldsymbol{\eta}))\}\\
	& - \Lambda_{i}(Y_i;\beta, \gamma, \boldsymbol{\eta}, g) ], \nonumber 
\end{align}
where $\Lambda_{i}(t;\beta, \gamma, \boldsymbol{\eta}, g)$ denotes the solution of ODE (\ref{model: ltm + cox}) parameterized by $(\beta, \gamma, \boldsymbol{\eta}, g)$ given covariates $X=X_i$ and $Z=Z_i$. The log-likelihood function (\ref{obj}) includes both finite-dimensional parameter $\beta$ and infinite-dimensional parameters $\gamma, \boldsymbol{\eta}, g$.

We propose a sieve MLE that maximizes the log-likelihood over a sequence of finite-dimensional parameter spaces that are dense in the original parameter space as the sample size increases. The sieve space can be chosen as linear spans of many types of basis functions with desired properties \citep{Chen_2007}. In particular, we construct the sieve space using polynomial splines due to their capacity in approximating complex functions and the simplicity of their construction. Under suitable smoothness conditions, the true functions $\gamma_0(\cdot)$, $\boldsymbol{\eta}_0(\cdot)$, and $g_0(\cdot)$ can be well approximated by some functions in the space of polynomial splines as defined in \citet[page 108, Definition 4.1]{schumaker_2007}. And there exists a group of spline bases such that functions in the space of polynomial splines can be written as linear combinations of the spline bases \citep[page 117, Corollary~4.10]{schumaker_2007}. Different groups of spline bases may be used for the estimation of different parameters $(\gamma, \boldsymbol{\eta})$ and $g$ because of their different domains. 

Specifically, the proposed sieve estimator is constructed as follows. Let $\mc{B}\subset {\mbox{\textbf{R}}}^{d_1}$ be the parameter space of $\beta$.
Let $\{B_j^1, 1\leq j \leq q_n^1\}$ and $\{B_j^2, 1\leq j \leq q_n^2\}$ be two groups of spline bases that are used for the estimation of parameters $(\gamma, \boldsymbol{\eta})$ and $ g$ respectively (see Section~\ref{s: asymptotics} for rigorous definitions). Overall, we wish to find $d_2+1$ members $(\gamma, \eta_1, \cdots, \eta_{d_2})$ from the space of polynomial splines associated with $\{B_j^1\}$, one member $g$ from that associated with $\{B_j^2\}$, along with $\beta \in \mc{B}$ to maximize the log-likelihood function~(\ref{obj}). Let $Z_{i0}= 1$, $Z_i =(Z_{i1},\cdots, Z_{id_2})^T$. Then the objective function can be written as
\begin{equation}
\label{eq: obj}
	l_n(\beta, a, b)=\frac{1}{n}\sum_{i=1}^n\left[\Delta_i\{X_i^T\beta + \sum_{l=0}^{d_2}\sum_{j=1}^{q_n^1}a_j^l B_j^1(Y_i)Z_{il}+ \sum_{j=1}^{q_n^2}b_j B_j^2(\Lambda_{i}(Y_i;\beta, a, b))\} - \Lambda_{i}(Y_i;\beta, a, b) \right], 
\end{equation}
where $a = \left(a_j^l \right)_{j=1,\cdots,q_n^1,l=0,\cdots,d_2}$ and $b = \left(b_j \right)_{j=1,\cdots,q_n^2}$ are the coefficients of the spline bases, and $\Lambda_{i}(t;\beta, a,b)$ is the solution of 
\begin{equation}
\label{ode: spline}
\left\{
\begin{array}{lr}
\Lambda'_{i}(t) = \exp(X_i^T\beta + \sum_{l=0}^{d_2}\sum_{j=1}^{q_n^1}a_j^l B_j^1(t)Z_{il}+\sum_{j=1}^{q_n^2}b_j B_j^2(\Lambda_{i}(t))), \\
\Lambda_{i}(0) =0.
\end{array}
\right.
\end{equation}
The proposed sieve estimators are given by $\hat{\beta}_n=\hat{\beta}$,   $\hat{\boldsymbol{\eta}}_{n}(\cdot)=\left(\sum_{j=1}^{q_n^1}\hat{a}_j^1 B_j^1(\cdot), \dots, \sum_{j=1}^{q_n^1}\hat{a}_j^{d_2} B_j^1(\cdot)\right)$, $\hat{\gamma}_n(\cdot)=\sum_{j=1}^{q_n^1}\hat{a}_j^0 B_j^1(\cdot)$, 
and $\hat{g}_n(\cdot)=\sum_{j=1}^{q_n^2}\hat{b}_j B_j^2(\cdot)$, where $(\hat{\beta}, \hat{a}, \hat{b})$ maximizes the objective function~(\ref{eq: obj}). 

Note that the objective function (\ref{eq: obj}) contains the solution of a parameterized ODE (i.e. (\ref{ode: spline})), and this is different from most traditional optimization problems. In particular, it is nontrivial to evaluate the objective function and its gradient with respect to parameters when there is no closed-form solution for the ODE. To address this optimization challenge, we develop a gradient-based optimization algorithm by taking advantage of local sensitivity analysis \citep{Dickinson_1976, petzold2006sensitivity} and well-implemented ODE solvers. Specifically, we evaluate the objective function and its gradient as follows:
\begin{enumerate}
	\item we numerically calculate $\Lambda_{i}(Y_i;\beta, a,b)$ by solving (\ref{ode: spline}) 
		given the current parameter estimates $\beta$, $a$, $b$ and covariates $X_i$, $Z_i$, the initial value at $t_0=0$, and the evaluating time $t_1=Y_i$; 
	\item we evaluate the derivative of $\Lambda_{i}(Y_i;\beta, a,b)$ with respect to the parameters $\beta$, $a$ and $b$, denoted as $\Lambda'_{i\beta}(Y_i;\beta, a,b)$, $\Lambda'_{ia}(Y_i;\beta, a,b)$, and $\Lambda'_{ib}(Y_i;\beta, a,b)$, through solving another ODE which is derived by local sensitivity analysis, and calculate the gradient of the objective function by the chain rule.
\end{enumerate}

We summarize the results of the local sensitivity analysis in the following, and provide detailed derivations in the Supplemental Material.
The local sensitivity analysis is a technique that studies the rate of change in the solution of an ODE system with respect to the parameters. There are two ways to obtain the sensitivity: forward sensitivity analysis and adjoint sensitivity analysis. Both of them require solving another ODE with some fixed initial value. 
For example, we consider to compute the gradient of $\Lambda(y; \theta)$ with respect to its parameter $\theta$, where $\Lambda(t;\theta)$ is the solution of (\ref{ode: spline}) and $\theta$ consists of parameters $\beta$, $a$, and $b$ in our case. For presentation simplicity, we denote the right-hand side of (\ref{ode: spline}) by the function $f(t,\Lambda;\theta)$, i.e.
\[f(t, \Lambda; \theta) = \exp(X^T\beta + \sum_{l=0}^{d_2}\sum_{j=1}^{q_n^1}a_j^l B_j^1(t)Z_{j}+\sum_{j=1}^{q_n^2}b_j B_j^2(\Lambda)),\]
and its partial derivative with respect to $\theta$ and $\Lambda$ by $f'_{\theta}$ and $f'_{\Lambda}$ respectively. In forward sensitivity analysis, it can be shown that the partial derivative of $\Lambda(y;\theta)$ with respect to $\theta$ is given by the solution of (\ref{I}) at $t=y$, i.e. $\Lambda'_{\theta}(y;\theta)=F_1(y)$ with $F_1$ satisfying
\begin{equation}
\label{I}
\left\{
\begin{array}{lr}
F'_{1}(t) = f'_{\theta}(t, \Lambda; \theta) + f'_{\Lambda}(t, \Lambda;\theta) F_1,\\
F_1(0) =0.
\end{array}
\right.
\end{equation}
In the alternative adjoint sensitivity analysis, we can show that the partial derivative can also be obtained by evaluating the solution of (\ref{II}) at $t=0$, i.e. $\Lambda'_{\theta}(y;\theta) = F_2(0)$ with $F_2$ satisfying 
\begin{equation}
\label{II}
\left\{
\begin{array}{lr}
(\kappa(t); F'_{2}(t))  = (-\kappa\cdot f'_{\Lambda}(t, \Lambda;\theta); -\kappa\cdot f'_{\theta}(t, \Lambda;\theta) ), \\
(\kappa(t); F_{2}(t))|_{t=y} = (1;\mb{0}).
\end{array}
\right.
\end{equation}
Thus, after plugging the form of $f(t,\Lambda;\theta)$ into either (\ref{I}) or (\ref{II}), we can obtain the gradients through solving the corresponding ODE. In Remark~\ref{rmk: ode}, we compare the computational complexity of forward and adjoint sensitivity analyses and provide a general guidance on which sensitivity analysis to use when computing gradients under different survival models. 

It is worth noting that the proposed estimation method can be easily implemented using existing computing packages. For example, the ``Optimization Toolbox'' in MATLAB contains ``fminunc'' for unconstrained optimization and ``fmincon'' for constrained optimization; both require initialization and the objective function. In our implementation, we also provide evaluation of the gradient for faster and more reliable computations. In particular, we compute both the objective function and the gradient by well-implemented ODE solvers in MATLAB. 
In addition, we construct the sieve space using B-splines for its numerical simplicity, whose implementation is available in the ``Curve Fitting Toolbox'', and utilize the ``Parallel Computing Toolbox'' to help parallel the for-loops.

\begin{remark}
\label{rmk: ode}
In general, forward sensitivity analysis is computationally more efficient when the dimension of the ODE system is relatively large and the number of parameters is small, while adjoint sensitivity analysis is best suited in the complementary scenario. See \citet{Dickinson_1976} and \citet{petzold2006sensitivity} for more details. For an ODE model whose right-hand side does not involve the cumulative hazard, such as (\ref{cox ode}), we can combine ODEs for $n$ individuals into a large ODE system with $n$ dimensions, which  
is larger than the number of parameters and then the forward sensitivity analysis is preferred. For a more general ODE model such as (\ref{transformation ode}) where the size of the ODE system is 1 and the number of parameters increases as the sample size $n$ grows, 
we use the adjoint sensitivity analysis along with parallel computing. 
\end{remark}


\begin{remark}
The proposed sieve MLE can also be applied to many existing models. For example, for the time-varying Cox model where $q(\cdot)=1$, we can remove the function $g(\cdot)$ from the objective function~(\ref{obj}). For the semi-parametric AFT model where $Z$ is not considered and $\alpha(\cdot)=1$, we can just keep parameters $\beta$ and $g(\cdot)$ in (\ref{obj}). For the linear transformation model, if either $q(\cdot)$ or $\alpha(\cdot)$ is specified, we can replace the corresponding term in (\ref{obj}) with the specified finite-dimensional parametric form. Also note that in comparison to existing estimation methods in Section \ref{s:estimation}, the proposed estimation method allows parallel computing, which is especially important for large-scale applications. Specifically, since the log-likelihood of each individual only depends on its own observations, the evaluation for independent data points can be carried out simultaneously. Further, comparing with the NPMLE where the number of optimization parameters is linear in $n$ \citep{murphy1997maximum, zeng2007maximum}, the number of optimization parameters used in sieve MLE increases more slowly with the sample size.
\end{remark}

\begin{remark}
The objective function~(\ref{eq: obj}) is convex with respect to $\beta$ and $a$ for the (time-varying) Cox model, where the parameter $b$ is not included, and the global optimum can be achieved quickly. 
For the general case, the objective function is nonconvex and the optimization algorithm may converge to a local optimum. Nevertheless, based on our extensive simulation studies, the algorithm generally performs well with appropriately chosen initialization, such as initializing the algorithm with the estimates from the Cox model. 
\end{remark}

\begin{remark}
Note that different identifiability conditions are required for different survival models. Thus, we need to add corresponding constraints in the optimization algorithm. For the general ODE model~(\ref{model: ltm + cox}) where both covariates $X$ (with time-independent effects) and $Z$ (with time-varying effects) are considered, two groups of parameters $(\beta, \gamma, g, \boldsymbol{\eta})$ and $(\tilde{\beta}, \tilde{\gamma}, \tilde{g}, \tilde{\boldsymbol{\eta}})$ give the same survival distribution if only if $\beta = \tilde{\beta}$, $\gamma=\tilde{\gamma}+c$, $g=\tilde{g}-c$, and $\boldsymbol{\eta}=\tilde{\boldsymbol{\eta}}$ for some constant $c$. To guarantee the identifiability, we can constrain either the value of $\gamma(\cdot)$ at a fixed time point or the norm of $\gamma(\cdot)$, in which the former leads to a linear constraint on the coefficients of spline bases.  For the linear transformation model where the time-varying effects are not considered, parameters $(\beta, \gamma, g)$ are identifiable up to two scaling factors as shown in Proposition~\ref{identi}.  Then we can either constrain the first element of $\beta$ to be $1$ \citep{khan2007partial, Song_2006}, which can be naturally achieved by arranging covariates if we know which covariate has a non-zero effect, or set $\|\beta\|=1$.  The constraint for $\gamma$ can be similarly added as for the general ODE model (\ref{model: ltm + cox}).  In our implementation, we choose to use two linear constraints, i.e. set the first element of $\beta$ to 1 and $\gamma(0)=0$ for simplicity in optimization.
\end{remark}



\subsection{Time-varying Cox model}
\label{subs: tinv cox}

We generate event times from the model \[\Lambda'_{x}(t) = \alpha(t) \exp(\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \eta(t) x_5 ),\]
where $(x_1, x_2, x_3, x_4, x_5)$ follows a multivariate normal distribution with mean $0$ and autoregressive covariance, $\beta_1=\beta_4=1$, and $\beta_2=\beta_3=-1$. Let $\eta(t) = \sin(\frac{3}{4}\pi t)$ be a time-varying coefficient for $x_5$ and the coefficients of all other covariates be time-independent. The baseline hazard $\alpha(t)$ is set to 0.5. The censoring times are generated from an independent uniform distribution $U(0, 3)$, which leads to a censoring rate around 50\%. The sample size $N$ varies from 500 to 8000. We fit both the baseline hazard function $\alpha(t)$ and time-varying coefficient $\eta(t)$ by cubic B-splines with 6 interior knots. The interior knots equally separate the observed event times. We compare the estimation accuracy and the computing time of the proposed sieve MLE with those of the partial likelihood-based estimator implemented in the ``coxph'' function in R with the ``tt'' argument set as the same cubic B-spline transformation of time.




\begin{table}[!t]
    \begin{center}
    \caption{Simulation results under time-varying Cox model.}
    \label{tab:time-varying cox1}
    \begin{threeparttable}
    \begin{tabular}{cc|cccc|cccc|cc}
    \toprule
        N & Method & \multicolumn{4}{c}{$\beta_1=1$} & \multicolumn{4}{c}{$\beta_2=-1$}& \multicolumn{2}{c}{IMSE($\eta(t)$)}\\
        & & Bias & SE & ESE & CP & Bias & SE & ESE & CP & Mean & SD \\
        \midrule
        \multirow{2}{*}{500} &  ODE & .014& .092& .092& .99& -.026& .104&  .111& .95&  .113 & .083\\
        
        & Cox(R) & .009& .092& .091 & .97 & -.022& .104& .103 & .96 & .155 &.136 \\ 
        \hline
       \multirow{2}{*}{1000} &  ODE & .004& .055& .066 & .99 & -.009& .079& .075 & .94 & .051 & .031\\
        
        & Cox(R) & .001& .055& .063 & .99 & -.007& .078& .072 & .91 &  .089 & .072  \\ 
		\hline
       \multirow{2}{*}{2000} &  ODE & .000& .049& .046 & .92 & -.004& .047& .051 & .98 &  .025 & .015\\
        & Cox(R) & -.001& .049& .045 & .92 & -.003& .046& .050 & .97 & .039 & .024  \\ 
		\hline
       \multirow{2}{*}{4000} &  ODE & .002& .030& .032 & .96 & .008& .031& .036 & .97 &  .013 & .007\\
        & Cox(R) & .002& .030& .031 & .94 & .009& .031& .035 & .97 & .021 & .016  \\ 
        \hline
        \multirow{2}{*}{8000} &  ODE & .000& .023& .022 & .96 & -.001& .025&  .025& .94 &  .008 & .004 \\
        & Cox(R) & .000& .023& .022 & .96 & .000& .025& .025 & .94 & .013 & .009 \\
        \bottomrule
    \end{tabular}
    \begin{tablenotes}
  \item  \footnotesize{Bias is the difference between the mean of estimates and the true value, SE is the sample standard error of the estimates, Mean is the mean of IMSE, and SD is the standard deviation of IMSE. ESE is the mean of the standard error estimators by inverting the estimated information matrix of all parameters, including the coefficients of spline basis, and CP is the corresponding coverage proportion of 95\% confidence intervals.
  }
  \end{tablenotes}
  \end{threeparttable}
  \end{center}
\end{table}



\begin{figure}[!ht]
    \centering
    \scalebox{0.48}
    {\includegraphics{fig/cox_tv_N8000_summary.eps}}
    \caption{True $\alpha_0(t)$ and mean of $\hat{\alpha}(t)$ (left); true $\eta(t)$ and mean of $\hat{\eta}(t)$ (middle) with the sample size $N=8000$; mean computing time (right) with respect to various sample size.} 
    \label{fig:cox_tv_N8000_ode}
\end{figure}

Table~\ref{tab:time-varying cox1} summarizes the estimates of regression coefficients $\beta_1$ and $\beta_2$ based on 100 replications. The estimates of the other two regression coefficients $\beta_3$ and $\beta_4$ perform similarly, and the results are included in the Supplemental Material. For the time-varying coefficient $\eta(t)$, we report the integrated mean square error (IMSE), which is the weighted sum of mean square error (MSE) of pointwise estimates over simulated time points from $0$ to $2$. As one can see, the mean and standard deviation of IMSE of the proposed sieve estimator decrease as the sample size increases. Remarkably, they are consistently smaller than those of the partial likelihood-based estimator. For time-independent coefficients, the proposed sieve estimator performs as well as the partial likelihood-based estimator. The mean of the standard error estimator, which is obtained by inverting the estimated information matrix of all parameters including the coefficients of spline bases, is approximate to the sample standard error, and the corresponding 95\% confidence interval achieves a proper coverage proportion.  From the left and middle panels of Figure~\ref{fig:cox_tv_N8000_ode}, we can see that the means of the estimated $\alpha(t)$ and $\eta(t)$ are close to the true functions, and the 95\% pointwise confidence bands cover the true functions well.

It is also worth noting that, in comparison to the partial likelihood-based estimation method whose computing time increases quickly as the sample size grows, the proposed estimation method is computationally more efficient, especially when the sample size is large (see the right panel of Figure~\ref{fig:cox_tv_N8000_ode}).


\section{Simulation Studies}
\label{s: simulation}
In this section, we use simulation studies to show the finite sample performance of the sieve MLE under the time-varying Cox model and the general linear transformation model.  


\input{simulation2.tex}

\subsection{Linear transformation model}
\label{subs: simu ltm}
We generate event times from the model $ \Lambda'_{x}(t) = q(\Lambda_x(t)) \exp(\beta_1x_1+\beta_2x_2+\beta_3x_3)\alpha(t)$. 
The covariates $(x_1, x_2)$ are bivariate normal with mean $0$, variance $0.5$, and covariance $0.2$; and $x_3$ is independent Bernoulli with probability $0.5$. We consider four different settings for $q(\cdot)$ and $\alpha(\cdot)$: 1) a constant $q(t)=1$ and a monotonic increasing $\alpha(t)=t^3$, in which case the Cox model is correctly specified; 2) a monotonic decreasing $q(t)=2/(1+t)$ and a constant $\alpha(t)=1$, where the AFT model is correctly specified; 3) an increasing $q(t)=\log(1+t)+2$ and an increasing $\alpha(t)=log(1+t)$; 4) the same $q(t)$ as setting 3) and a periodic function $\alpha(t)=1+\cos(\pi t+1)$ that decreases first and increases later. Note the Cox model and the AFT model are mis-specified in settings 3) and 4). Several other settings have also been considered but are included in the Supplemental Material for the sake of space.  In each setting, we generate the censoring time from an independent uniform distribution $U(0, c)$, where $c$ is chosen to achieve approximately 20-30\% censoring rates. The sample size $N$ is set to 4000. 

In setting 1), we compare the proposed sieve MLE for the ODE-Cox model, where the function $q$ is set to 1, with the partial-likelihood based estimator implemented using the R package {\it survival}. We fit $\alpha(t)$ by cubic B-splines with 6 interior knots that equally separate the time interval from $0$ to the maximum of observed times. In setting 2), we compare the proposed sieve MLE for the ODE-AFT model, where the function $\alpha$ is set to 1, with the rank-based estimation approach implemented using the R package {\it aftgee}. We fit $q(t)$ by cubic B-splines with 6 interior knots that equally separate the time interval from $0$ to the maximum of estimated cumulative hazards under the Cox model. 

In addition, in all settings, we also fit the general linear transformation model (ODE-Flex) where both $q(\cdot)$ and $\alpha(\cdot)$ are unspecified and compare the sieve MLE with the estimators under the Cox model and the AFT model. We constrain $\beta_1=1$ and $q(0)=1$ to guarantee identifiability for ODE-Flex, which introduces two linear equality constraints to the optimization problem. To make the comparison with ODE-Flex fair, the Cox and AFT models are also fitted with the constraint $\beta_1=1$.


\begin{table}[!t]
    \begin{center}
    \caption{Estimates of regression coefficients under correctly-specified ODE-Cox with $q(\cdot) \equiv 1$ and ODE-AFT with $\alpha(\cdot)\equiv 1$. Bias, SE, ESE and CP contain the same meanings as those in Table \ref{tab:time-varying cox1}.}
    \label{tab:ltm}
    \begin{threeparttable}
    \begin{tabular}{cc|cccc|cccc|cccc}
    \toprule
        \small{Setting} & Method & \multicolumn{4}{c}{$\beta_1=1$} & \multicolumn{4}{c}{$\beta_2=1$} & \multicolumn{4}{c}{$\beta_3=1$}\\
        & & \small{Bias} & \small{SE} & \small{ESE} & \small{CP} & \small{Bias} & \small{SE} & \small{ESE} & \small{CP} &  \small{Bias} & \small{SE} & \small{ESE} & \small{CP} \\
        \midrule
        \multirow{2}{*}{1)} &  \small{Cox(R)} & \small{.005} & \small{.033} & \small{.032} & \small{.93} & \small{.001} & \small{.034} & \small{.032} & \small{.96}  & \small{.001} & \small{.037} & \small{.039} & \small{.97} \\
        & \small{ODE-Cox} & \small{.005} & \small{.034} & \small{.032} & \small{.93} & \small{.001} & \small{.034} & \small{.032} & \small{.96} & \small{.002} & \small{.038} & \small{.039} & \small{.97} \\ 
        \hline
       \multirow{2}{*}{2)} &  \small{AFT(R)} & \small{.006} & \small{.041} & \small{.040} & \small{.95} & \small{.000} & \small{.041} & \small{.040} & \small{.95} & \small{.002} & \small{.050} & \small{.051} & \small{.94}\\
        & \small{ODE-AFT} & \small{.004} & \small{.037} & \small{.039} & \small{.94} & \small{-.001} & \small{.041} & \small{.038} & \small{.95} & \small{-.003} & \small{.049} & \small{.049} & \small{.93} \\	
        \bottomrule
    \end{tabular}
    \begin{tablenotes}
  \item  \footnotesize{Setting 1): the Cox model is correctly specified. Setting 2): the AFT model is correctly specified.
  }
  \end{tablenotes}
  \end{threeparttable}
  \end{center}
\end{table}

\begin{table}[!t]
    \begin{center}
    \caption{Estimates of regression coefficients under the general linear transformation model ODE-Flex with both $q(\cdot)$ and $\alpha(\cdot)$ unspecified. Bias, SE, ESE and CP contain the same meanings as those in Table \ref{tab:time-varying cox1}.}
    \label{tab:ltm2}
    \begin{threeparttable}
    \begin{tabular}{p{3cm}c|cccc|cccc}
    \toprule
        \multicolumn{1}{c}{Setting} & Method & \multicolumn{4}{c}{$\beta_2=1$} & \multicolumn{4}{c}{$\beta_3=1$}\\
        & & Bias & SE & ESE & CP &  Bias & SE & ESE & CP \\
        \midrule
        \multirow{2}{*}{\parbox{2.5cm}{\small{1) Cox is correctly specified}}} & Cox(R)$^\star$ & .001 & .034 & .031 & .95 & .000 & .036 & .038 & .97 \\
        & ODE-Flex &  -.004 & .047 & .045 & .95 & -.003 & .042 & .045 & .97 \\ 
        \hline
        \multirow{2}{*}{\parbox{2.5cm}{\small{2) AFT is correctly specified}}} & AFT(R)$^\star$ & .003 & .039 & .038 & .97 & .002 & .050 & .051 & .93 \\
        & ODE-Flex &  -.007 & .057 & .056 & .94 & -.005 & .056 & .058 & .94 \\
        \hline
       \multirow{3}{*}{\parbox{2.8cm}{\small{3) Cox and AFT are mis-specified}}} & Cox(R)$^\star$ & .128 & .034 & .032 & .02 & .093 & .035 & .038 & .29 \\
        & AFT(R)$^\star$ & -.618 & .017 & .018 & .00 & -.443 & .023 & .024 & .00 \\
        & ODE-Flex &  -.003 & .042 & .038 & .94 & .000 & .037 & .038 & .93 \\
       \hline
       \multirow{3}{*}{\parbox{2.8cm}{\small{4) Cox and AFT are mis-specified}}} & Cox(R)$^\star$ & .131 & .032 & .032 & .01 & .093 & .035 & .039 & .29 \\
        & AFT(R)$^\star$ & .301 & .035 & .035 & .00 & .221 & .046 & .048 & .01 \\
        & ODE-Flex & -.004 & .039 & .036 &.92 & -.003 & .035 &.037 & .97\\
        \bottomrule
    \end{tabular}
    \begin{tablenotes}
  \item  \footnotesize{$^\star$ Cox and AFT models are fitted with the constraint $\beta_1=1$ to make a fair comparison with ODE-Flex.
  }
  \end{tablenotes}
  \end{threeparttable}
  \end{center}
\end{table}



\begin{figure}[!t]
    \centering
    \scalebox{0.48}
    {\includegraphics{fig/fitted_curves_N4000_setting1.eps}}
    \scalebox{0.48}
    {\includegraphics{fig/fitted_curves_N4000_setting3.eps}}
    \scalebox{0.48}
    {\includegraphics{fig/fitted_curves_N4000_setting8.eps}}
    \scalebox{0.48}
    {\includegraphics{fig/fitted_curves_N4000_setting9.eps}}
    \caption{The solid blue curves are the rue $q(\cdot)$, $\alpha(\cdot)$, and the survival functions respectively (from left to right).  The solid red curves are the corresponding estimated $\hat{q}(\cdot)$, $\hat{\alpha}(\cdot)$ and survival functions under the general linear transformation model.  The dashed yellow curves represent 95\% pointwise confidence bands.  From top to bottom, the four rows correspond to settings 1)-4) respectively.}
    \label{fig:ltm}
\end{figure}


Tables~\ref{tab:ltm} and \ref{tab:ltm2} summarize the estimates of regression coefficients based on $100$ replications. 
Table~\ref{tab:ltm} indicates that when either the Cox model or the AFT model is correctly specified, the sieve estimator for the corresponding correctly specified ODE model achieves similar performance as the partial-likelihood based estimator for the Cox model or the rank-based estimator for the AFT model.  From Table~\ref{tab:ltm2}, we can see that the bias of the estimator for the general linear transformation model, i.e. ODE-Flex, is nearly negligible in all settings. The standard error estimators are close to the sample standard errors, and the corresponding 95\% confidence intervals achieve a reasonable coverage proportion. When the Cox model and the AFT model are correctly specified, their estimators achieve smaller standard errors than those for ODE-Flex, which is expected because both $q(\cdot)$ and $\alpha(\cdot)$ are unspecified in ODE-FLex. However, when the Cox model and the AFT model are mis-specified, such as in settings 3) and 4), the estimates under the Cox model and the AFT model clearly deviate from the true values. 
Figure~\ref{fig:ltm} shows the mean of $\hat{\alpha}(\cdot)$ and $\hat{q}(\cdot)$ respectively, and the mean of estimated survival functions for two representative groups of covariates, where $x_1$ and $x_2$ are set to the mean $0$ and the binary covariate $x_3$ is either $0$ or $1$.  As one can see from  Figure~\ref{fig:ltm}, the means of $\hat{\alpha}(\cdot)$, $\hat{q}(\cdot)$, and the estimated survival functions under the general linear transformation model are all close to the true functions.


