Non-Smooth Backfitting for Excess Risk Additive Regression Model with Two Survival Time-Scales

TL;DR

This paper introduces a novel backfitting algorithm for a complex non-parametric survival model with two time-scales, avoiding smoothing and providing large sample properties and confidence bands.

Contribution

It proposes a smoothing-free backfitting method for a two time-scale survival model with additive covariate effects, extending existing estimation techniques.

Findings

Successfully estimates cumulative intensities without smoothing.

Provides large sample properties and confidence bands.

Applied to myocardial infarction data to distinguish effects.

Abstract

We present a new backfitting algorithm estimating the complex structured non-parametric survival model of Scheike (2001) without having to use smoothing. The considered model is a non-parametric survival model with two time-scales that are equivalent up to a constant that varies over the subjects. Covariate effects are modelled linearly on each time scale by additive Aalen models. Estimators of the cumulative intensities on the two time-scales are suggested by solving local estimating equations jointly on the two time-scales. We are able to estimate the cumulative intensities solving backfitting estimating equations without using smoothing methods and we provide large sample properties and simultaneous confidence bands. The model is applied to data on myocardial infarction providing a separation of the two effects stemming from time since diagnosis and age.

Figures

Tables

Table 1. Table 1: Bias of backfitting algorithm for sample sizes n = 100 , 200 , 400 𝑛 100 200 400 n=100,200,400 for the age and time component for selected ages and time points. Based on 1000 realisations.

age	n=100	n=200	n=400
$6.717$	$-0.001$	$0.006$	$-0.004$
$13.788$	$0.009$	$0.003$	$-0.006$
$20.859$	$0.018$	$0.001$	$0.002$
$27.929$	$0.027$	$0.004$	$0.010$
$35$	$0.078$	$0.006$	$0.013$
time	n=100	n=200	n=400
$0.96$	$0.018$	$0.009$	$0.006$
$1.97$	$0.015$	$0.007$	$0.005$
$2.98$	$0.009$	$0.005$	$0.003$
$3.99$	$0.005$	$0.002$	$0.002$
$5$	$0$	$0$	$0$

Table 2. Table 2: Uncertainty estimated from bootstrap for sample sizes n = 100 , 200 , 400 𝑛 100 200 400 n=100,200,400 for the age and time component for selected ages and time points. Based on 1000 realisations and a bootstrap with 100 repetitions. mean of estimated standard errors (mean se), standard deviation of estimates (sd) and 95 % pointwise coverage (cov).

n	age	mean se	sd	cov	time	mean se	sd	cov
$100$	$6.717$	$0.224$	$0.231$	$0.912$	$0.96$	$0.044$	$0.045$	$0.954$
$100$	$13.788$	$0.297$	$0.298$	$0.935$	$1.97$	$0.039$	$0.04$	$0.946$
$100$	$20.859$	$0.351$	$0.357$	$0.943$	$2.98$	$0.032$	$0.034$	$0.951$
$100$	$27.929$	$0.391$	$0.402$	$0.938$	$3.99$	$0.024$	$0.024$	$0.966$
$100$	$35$	$0.460$	$0.464$	$0.932$	$5$	$0.016$	$0.017$	$0.874$
$200$	$6.717$	$0.158$	$0.155$	$0.94$	$0.96$	$0.031$	$0.031$	$0.951$
$200$	$13.788$	$0.207$	$0.206$	$0.942$	$1.97$	$0.027$	$0.027$	$0.960$
$200$	$20.859$	$0.243$	$0.237$	$0.948$	$2.98$	$0.022$	$0.022$	$0.966$
$200$	$27.929$	$0.271$	$0.262$	$0.945$	$3.99$	$0.017$	$0.017$	$0.972$
$200$	$35$	$0.328$	$0.329$	$0.933$	$5$	$0.011$	$0.012$	$0.933$
$400$	$6.717$	$0.114$	$0.118$	$0.948$	$0.96$	$0.022$	$0.022$	$0.951$
$400$	$13.788$	$0.148$	$0.153$	$0.946$	$1.97$	$0.019$	$0.019$	$0.957$
$400$	$20.859$	$0.173$	$0.18$	$0.937$	$2.98$	$0.015$	$0.015$	$0.960$
$400$	$27.929$	$0.192$	$0.196$	$0.943$	$3.99$	$0.012$	$0.012$	$0.970$
$400$	$35$	$0.235$	$0.245$	$0.934$	$5$	$0.008$	$0.008$	$0.950$

Table 3. Table 3: Coverage of confidence bands estimated from bootstrap for sample sizes n = 100 , 200 , 400 𝑛 100 200 400 n=100,200,400 for the age and time component. Based on 1000 realisations and a boostrap with 100 repetitions.

n	coverage (age)	coverage (time)
$100$	$0.797$	$0.792$
$200$	$0.912$	$0.915$
$400$	$0.952$	$0.939$

Equations

$${\lambda }_i(t)$$

$$=X_i(t)\alpha (t)+Z_i(t)\beta (t+a_i),(0\le t\le t_{max}),$$

$$X_{il}=Z_{il},l\le d.$$

$$f_l(s,u)={\left(0,\cdots ,0,c_{l}s,0,\cdots ,0,-c_l(u-a_0),0,\cdots ,0)\right)}^T,(l=1,\dots ,d).$$

$$G=\left\{g:\left(\begin{array}{c}A\\ B\end{array}\right)\mapsto \left(\begin{array}{c}A\\ B\end{array}\right)+h|h\in Lin(f_1,\dots f_d)\right\}.$$

$$A_l(t_{max})={\int }_0^{t_{max}}{\alpha }_l(s)\mathrm{d}s=0,(l=1,\dots ,d),$$

$$\arg \underset{\overline{A},\overline{B}}{\min }\sum _i\int {\left\{{\int }_0^t\mathrm{d}{N}_i(s)-\sum _j{\int }_0^t{X}_{ij}(s)d{\overline{A}}_j(s)-\sum _k{\int }_0^t{Z}_{ik}(s)d{\overline{B}}_k^{a_i}(s)\right\}}^2\mathrm{d}t,$$

$$\arg \underset{\overline{A},\overline{B}}{\min }\sum _i\int {\left\{{\int }_0^t\mathrm{d}{N}_i(s)-{\int }_0^t{X}_i(s)d\overline{A}(s)-{\int }_0^t{Z}_i(s)d{\overline{B}}^{a_i}(s)\right\}}^2\mathrm{d}t.$$

$$\sum _i{X}_i(t{)}^T\left\{d{N}_i(t)-X_i(t)d\hat{A}(t)-Z_i(t)\mathrm{d}{\hat{B}}^{a_i}(t)\mathrm{d}t\right\}=0,$$

$$\sum _i{Z}_i^{-a_i}(a{)}^T\left\{d{N}_i^{-a_i}(a)-Z_i^{-a_i}(a)d\hat{B}(a)-X_i^{-a_i}(a)\mathrm{d}{\hat{A}}^{-a_i}(a)\right\}=0.$$

$$\sum _i{X}_i(t{)}^{T}d{N}_i(t)-\sum _i{X}_i(t{)}^T{Z}_i(t)\mathrm{d}{\hat{B}}^{a_i}(t)=X(t{)}^{T}X(t)\mathrm{d}\hat{A}(t),$$

$$\sum _i{Z}_i^{-a_i}(a{)}^{T}d{N}_i^{-a_i}(a)-\sum _i{Z}_i^{-a_i}(a{)}^T{X}_i^{-a_i}(a)\mathrm{d}{\hat{A}}^{-a_i}(a)=Z^{-a_{\bullet }}(a{)}^T{Z}^{-a_{\bullet }}(a)d\hat{B}(a).$$

$$\hat{A}(t)$$

$$\hat{B}(a)$$

$$E_1(s|u)$$

$$E_2(u|s)$$

$${\lambda }_i(t)=Y_i(t)\{\alpha (t)+\beta (a_i+t)\},$$

$$E_1(s|u)$$

$$E_2(u|s)$$

$$\left(\begin{array}{c}\hat{A}\\ \hat{B}\end{array}\right)=\left(\begin{array}{c}{\int }_0^{t}X(s{)}^{-}dN(s)\\ {\int }_{a_0}^a{Z}^{-a_{\bullet }}(u{)}^{-}d{N}^{-a_{\bullet }}(u)\end{array}\right)+\left(\begin{array}{cc}0& -E_1\\ -E_2& 0\end{array}\right)\times \left(\begin{array}{c}\hat{A}\\ \hat{B}\end{array}\right),$$

$$\hat{\theta }=\hat{m}+E\hat{\theta },$$

$$\hat{\theta }=\left(\begin{array}{c}\hat{A}\\ \hat{B}\end{array}\right),\hat{m}=\left(\begin{array}{c}{\int }_0^{t}X(s{)}^{-}dN(s)\\ {\int }_{a_0}^a{Z}^{-a_{\bullet }}(u{)}^{-}d{N}^{-a_{\bullet }}(u)\end{array}\right),E=\left(\begin{array}{cc}0& -E_1\\ -E_2& 0\end{array}\right).$$

$$\hat{\theta }=(I-E{)}^{-1}\hat{m}.$$

$$E_2{f}_{11}(u)$$

$$E_1{f}_{12}(s)$$

$$\int E_1(s|u)\mathrm{d}u=\sum _i{\int }_{a_i}^{a_i+s}\frac{1}{{\sum }_{i^{\prime }}{Y}_{i^{\prime }}(u-a_i)}{Y}_i^{-a_i}(u)\mathrm{d}u={\int }_0^s\frac{{\sum }_i{Y}_i(t)}{{\sum }_{i^{\prime }}{Y}_{i^{\prime }}(t)}\mathrm{d}t=s.$$

$$E\left(\begin{array}{c}{f}_{11}\\ f_{12}\end{array}\right)=\left(\begin{array}{c}-E_1{f}_{12}\\ -E_2{f}_{12}\end{array}\right)=\left(\begin{array}{c}{f}_{11}\\ f_{12}\end{array}\right).$$

$$\hat{A}(t)$$

$$\hat{B}(a)$$

$$\hat{\theta }=\hat{m}+\overline{E}\hat{\theta },\overline{E}=\left(\begin{array}{cc}0& -E_1\\ -{\overline{E}}_2& 0\end{array}\right),$$

$${\hat{\theta }}_0=(I-\tilde{\Pi })\hat{\theta },$$

$$\tilde{\Pi }\left(\begin{array}{c}{h}_1(t)\\ h_2(a)\end{array}\right)=\left(\begin{array}{c}t{h}_1^{d_p}(t_{max})t_{max}^{-1}\\ -(a-a_0)h_1^{d_q}(t_{max})t_{max}^{-1}\end{array}\right).$$

$${\hat{\theta }}_0+Lin(f_1,\dots f_d),$$

$$K=\{h=(h_1,\dots ,h_d,0,\dots ,0)|h_l:\Re \to \Re ,x\mapsto c_{l}x,c_l\in \Re ,l=1,\dots ,d\}.$$

Full text

Key words: Aalen model, counting process, disability model, illness-death model, generalized additive models, multiple time-scales, non-parametric estimation, varying-coefficient models.

1 Introduction

In many bio-medical applications in survival analysis it is of interest and needed to use multiple time-scales. A medical study will often have a follow-up time (for example time since diagnosis) for patients of different ages, and here both time-scales will contain important but different information about how the risk of, for example, dying is changing. We therefore consider the situation with two time-scales that are equivalent up to a constant for each individual, such as for example follow-up time and age. One may see this as arising from the the illness-death model, or the disability model, where the additional time-scale may be duration in the illness state of the model; see Keiding (1991) for a general discussion of these models. There is rather limited work on how to deal with multiple time-scales in a biomedical context, see for example Oakes (1995); Iacobelli & Carstensen (2013) and Duchesne & Lawless (2000) and references therein. We present a non-parametric regression approach with two time-scales where each time-scale contribute additively to the mortality. The regression setting models the effect of covariates by additive Aalen models on each time-scale (Aalen, 1989; Huffer & McKeague, 1991; Andersen et al., 1993; Martinussen & Scheike, 2006). This allows covariates to have effects that vary on two different time-scales. In a motivating example we consider patients that experience myocardial infarction, and aim at predicting the intensity considering the two time-scales age and time since myocardial infarction. As a consequence, we can make survival predictions for patients given their age at diagnosis. This model was considered previously by Scheike (2001) where estimation was based on smoothing for one of the time-scales. A study closely related to ours is Kauermann & Khomski (2006) who studied the two most common time scales: age and duration. The underlying technical setting of Kauermann & Khomski (2006) was a multiplicative hazard model without covariates that is estimated via splines. In contrast our approach is an additive hazard model including covariates and estimating without smoothing. Alternative smoothing methodologies to multiplicative hazard estimation includes Linton et al. (2003); Huang (1999); Hastie & Tibshirani (1986); Lin et al. (2016). None of the known multiplicative hazard approaches including the ones mentioned above are able to estimate without smoothing, include time varying covariate-effects, or are able to provide simultaneous confidence bands as the additive approach of this paper does provide. We do know that smoothing improves efficiencies of cumulatively estimated quantities, see Guillen et al. (2007) for the simplest possible case. However, smoothing is also a complexity and experts applying survival analysis have developed a practical way of smoothing by eye the underlying rough non-parametric estimators of Kaplan & Meier (1958); Nelson (1972). The advantage of providing estimators without smoothing is that there can be no confusion from the complicated process of picking the smoothing procedure first and the amount of smoothing after that. Even if a smoothing approach is eventually used, then the smoothing free procedure would always count as a benchmark approach to check whether something went wrong during the smoothing. Our backfitting approach is different from standard backfitting in regression, see for example the smooth additive backtiffing approach of Mammen et al. (1999), where data is projected down via a smoothing kernel onto an additive subspace. In the backfitting approach of this paper, the non-parametric dynamics is only taking place in the two time directions, and the end result is therefore closer to the classical approach of Nelson (1972) with a non-smooth estimator of the dynamics in the one-dimensional time axis. What is obtained through Aalen’s additive hazard regression model on two time axis is that the dynamics of the two time effects are adjusted for covariaties in a way that keep the one-dimensional structure of the non-parametric dynamics. The expert user of survival methodology can therefore use the well developed intuition from looking at Nelson-Aalen estimators and Kaplan-Meier estimators when interpreting the empirical results based on the new methodology of this paper. Another advantage of estimating directly the cumulative hazards is that we are able to obtain a simple uniform asymptotic description of our estimators. We are thus able to construct confidence bands and intervals, that are based on bootstrapping the underlying martingales.

The paper is organised as follows. Section 2 presents the model via counting processes. Section 3 gives some least squares based local estimating equations that are solved to give simple explicit estimators of the non-parametric effects of the model. Based on these explicit estimators we are able to derive asymptotic results and provide the estimators with asymptotic standard errors. Sections 4-6 discusses how to solve the equations and compute the estimator practically and how deal with identifiability issues. Section 7 shows how the large sample properties may be derived and in Section 8 we construct confidence bands. Section 9 demonstrates the finite sample properties supporting Section 10 where we use our proposed methods in a worked example. Finally, Section 10 discusses some possible extensions.

2 Aalen’s Additive Hazard Model for Two Time-Scales

Let $N_{i}(t)\;\;i=1,...,n$ be $n$ independent counting processes that do not have common jumps and are adapted to a filtration that satisfy the usual conditions (Andersen et al., 1993). We assume that the counting processes have intensities given by

[TABLE]

where $\alpha=(\alpha_{1},\dots,\alpha_{p})$ and $\beta=(\beta_{1},\dots,\beta_{q})$ are tupels of one dimensional deterministic functions, $X_{i}^{T}(t)\in\Re^{p}$ and $Z_{i}^{T}(t)\in\Re^{q}$ are predictable cadlag covariate vectors with $X(t)$ and $Z(t)$ having almost surely full rank, and $a_{i}$ is a real-valued random variable observed at time $t=0$. If $Z_{i}(t)=0$ for all $t$, $a_{i}$ does not need to be observed.

The model is the sum of two Additive Alalen Models running on two different time scales, see also Scheike(2001). The two time-scales are $t$ and $a=t+a_{i}\in[a_{0},a_{max}]$ where the latter time-scale is specific to each individual and $a_{0}$ is some lower-limit that depends on the observed range of the second time-scale. Note, that no indicator variables are introduced but are absorbed in the covariates. In the illness-death model, say, $t$ might be time since diagnosis (duration) among subjects that have entered the illness stage of the model and $a_{i}$ could be the age when the transition to the illness stage occurred, such that $t+a_{i}$ is the age of the subject.

After introducing some notation we present an estimation procedure that leads to explicit estimators of $A(t)=\int_{0}^{t}\alpha(s)ds=(\int_{0}^{t}\alpha_{1}(s)ds,\dots,\int_{0}^{t}\alpha_{p}(s)ds)^{T}$ and $B(a)=\int_{a_{0}}^{a}\beta(u)du=(\int_{a_{0}}^{a}\beta_{1}(u)du,\dots,\int_{a_{0}}^{a}\beta_{q}(u)du)^{T}$. The cumulative effects have the advantage compared to $\alpha(s)$ and $\beta(a)$ that they may be used for inferential purposes since a more satisfactory simultaneous convergence can be established for these processes. We derive the asymptotic distribution for these estimators and a bootstrapping procedure quantifying the estimation uncertainty. Based on the cumulative intensity $A(t)$ one may estimate the intensity $\alpha(t)$ by smoothing techniques.

2.1  Notation Let $\Lambda_{i}(t)=\int_{0}^{t}\lambda_{i}(s)ds$ such that $M_{i}(t)=N_{i}(t)-\Lambda_{i}(t)$ are martingales. Let further $N(t)=(N_{1}(t),...,N_{n}(t))^{T}$ be the n-dimensional counting process, $\Lambda(t)=(\Lambda_{1}(t),...,\Lambda_{n}(t))^{T}$ is its compensator, such that $M(t)=(M_{1}(t),...,M_{n}(t))^{T}$ is an n-dimensional martingale, and define matrices $X(t)=(X_{1}(t),\ldots,X_{n}(t))^{T}$ and $Z(t)=(Z_{1}(t),\ldots,Z_{n}(t))^{T},$ with dimensions $n\times p$ and $n\times q$, respectively. The individual entry times are summarised in one vector $a_{\bullet}=(a_{1},\dots,a_{n})$. A superscript $a>0$ denotes a shift in the argument, i.e, for a generic function $f$, $f^{a}(y)=f(y+a)$. For a generic matrix $C(t)$, with $n$ rows $C_{i}(t)$, and a n-dimensional vector $v$, $C^{v}(t)$ is defined through shifting the rows: $C_{i}^{v}(t)=C_{i}(t+v_{i})$. For a generic matrix $C$, a minus superscript, $C^{-}$, denotes the Moore-Penrose inverse. An integral, $\int$, with no limits denotes integration over the whole range.

3 Identification of the entering nonparametric parameters

In many cases some covariates will enter both the $X$ and the $Z$ design. If this is the case, then the functions $\alpha$ and $\beta$ are not identified in model (2) – constants can be shifted for the components that share the same covariate without altering the intensity. Without loss of generality we assume that $X$ and $Z$ share the first $d$ $(0\leq d\leq\min(p,q))$ columns, i.e., for all $i=1,\dots,n$,

[TABLE]

We formulate the problem using group-theoretic arguments, see also Carstensen (2007); Kuang et al. (2008). Fix constants $c_{1},\dots,c_{d}$ and define $f_{l}$ as $\Re^{p+q}$ valued function having all entries but the $l^{\prime}th$ and the $(d+l)^{\prime}th$ equal zero:

[TABLE]

We define the group $G$ by

[TABLE]

The identification problem can be rephrased as that the intensity defined in (2) is a function of $(A,B)^{T}$, which is invariant to transformations $g\in G$. In the sequel we circumvent the identification issue by adding the following constraint

[TABLE]

noting that for any solution $(A_{0},B_{0})$ of model (2), there exists a unique solution $(A,B)=g(A_{0},B_{0})$ that fulfills (2). Clearly other choices are also possible.

4 Least squares minimisation ignoring the identificating of the nonparametric parameters

We split the identification challenge in two. First we estimate ignoring identification of the parameters, and then we show in next section how to identify the estimated parameters. In this section we therefore ignore the identification problem keeping in mind that the solutions below are hence not unique. We motivate our estimator $(\widehat{A},\widehat{B})$ via the following least squares criteria.

[TABLE]

where the integrals can be understood as Stieltjes integrals, noting that $X_{i}$ and $Z_{i}$ are left continuous. Minimisation runs over all possible integrators. One can already see that the minimiser, if it exists, will be a step-function, since $\int_{0}^{t}\mathrm{d}N_{i}(s)$ is a step function. To simplify notation we will generally work in matrix notation so that above minimisation criteria can also be written as

[TABLE]

Straight forward computations utilzing calculus of variations lead to $(\widehat{A},\widehat{B})$ solving the following first order conditions for all $t\in[0,t_{max}]$, $a\in[a_{0},a_{max}]$:

[TABLE]

Rearranging yields

[TABLE]

The last set of equations can be further rewritten to the backfitting equations

[TABLE]

where

[TABLE]

Remark 1

In the case with no covariates, i.e.,

[TABLE]

with $X_{i}(s)=Z_{i}(s)=Y_{i}(s)\in\Re$, the risk indicators are

[TABLE]

5 Establishing existence, identification and uniqueness of the estimator

In section 3 we outlined the identification problem but ignored it when establishing the estimator in the previous section. In this section we provide a fully identified estimator of our problem. When aiming to solve equations (3) and (4) the identification problem can no longer be ignored. In order to get a better grip of the situation we will now rewrite the backfitting equations as a linear operator equation. We can compress equations (3) and (4) into one matrix equation:

[TABLE]

where with some miss-use of notation $E_{l}f(\cdot)=\int E_{l}(\cdot,y)f(y)\mathrm{d}x,(l=1,2)$. Or even simpler

[TABLE]

with obvious notation, and linear operator $E$:

[TABLE]

Note that $\widehat{m}$ is composed of the marginal Aalen estimators of the two time scales, $t$ and $a$. Additionally, the operator $E$ is compact because it is the composition of an integral operator, which is compact, and a derivative operator, which is bounded. The operator $E$ being compact means that it can be arbitrarily close approximated by a finite dimensional matrix which simplifies both the numerical and theoretical considerations. If the eigenvalues of $E$ are bounded away from one, then, $(I-E)$ is invertible and we have

[TABLE]

Hence existence and uniqueness of our proposed estimator can be translated to properties of the eigenvalues of $E$. One can for instance easily verify that if some covariates are both in the $X$ and the $Z$ design, then $E$ will have an eigenvalue equal to one - as discussed in the following remark.

Remark 2

Consider the most simple case $1=d=p=q$, i.e., $\lambda_{i}(t)=Y_{i}(t)\{\alpha(t)+\beta(a_{i}+t)\}$. Given a constant $c\in\Re$, consider the pair of linear function $f_{1}=(f_{11},f_{12})^{T}$ with $f_{11}(s)=cs,\ f_{12}(u)=-c(u-a_{0})$, as defined in Section 3. Assuming that $\sum Y_{i}(s)$ and $\sum Y_{i}(u-a_{i})$ are bounded away from zero on the whole range $s\in[0,t_{max}],u\in[a_{0},a_{max}]$, one can easily verify that

[TABLE]

To see this, e.g., for the second equation, note

[TABLE]

Hence, we have

[TABLE]

So that one is clearly an eigenvalue of $E$ with corresponding eigenfunction $f_{1}=(f_{11},f_{12})^{T}$. In other words the identification issue of the model carries over to the estimator. With analogue arguments one can show that in the more general case the eigenspace corresponding to eigenvalue equal one includes the functions in $Lin(f_{1},\dots f_{d})$. Functions $f_{1},\dots,f_{d}$ are defined in Section 3.

We now utilize constraint (2) and incorporate it into new backfitting equations:

[TABLE]

where $\widehat{A}^{d_{q}}$ is the q-dimensional vector $\widehat{A}^{d_{q}}=(A_{1},\dots,A_{d},0,\dots,0)^{T}$. This translates to the new operator equation

[TABLE]

where $\overline{E}_{2}h(a)=\int E_{2}(a|s)dh(s)-(a-a_{0})h^{d_{q}}(t_{max})t_{max}^{-1}$. The next proposition states that the solutions of $\eqref{eq:operator2}$ include all relevant solutions of (5) and that every solution of $\eqref{eq:operator2}$ is a solution of $\eqref{eq:operator1}$.

Proposition 1

For every solution $\widehat{\theta}$ of (5), define

[TABLE]

where

[TABLE]

Then $\widehat{\theta}_{0}$ is a solution of $\eqref{eq:operator2}$ and

[TABLE]

are further solutions of (5). Reversly, for every solution ${\widehat{\theta}}_{0}$ of (8), all functions of the form (9) are solutions of (5).

The proof can be found in the appendix.

With Proposition 1 at hand it is justified to define our estimator as the solution of (8). We will now discuss existence and uniqueness of the solution of (8).

Note that $E$ is known and hence one can calculate a numerical approximation of its eigenvalues by working on a grid. Consider the sub-space

[TABLE]

It holds that $\overline{E}_{2}=E_{2}(I-\Pi)$, where $\Pi$ is a projection into $K$. We have $K\subseteq kern(I-E_{2})$. We can check whether $K$ equals $kern(I-E_{2})$. This can be done by calculating the dimension of the eigenspace of $E_{2}$ corresponding to an eigenvalue equal one. The dimension will be at least $d$. If it is exactly $d$, then $K=kern(I-E_{2})$.

The next proposition states that if $kern(I-E_{2})=K$, and $kern(I-E)=Lin(f_{1},\dots,f_{d})$, then both $I-\overline{E}_{2}$ and $I-\overline{E}$ are bijective.

Proposition 2

Assume that $E_{2}$ has Eigenvalue 1 with multiplicity $d$. Then, $(I-\overline{E}_{2})$ will be bijective. If furthermore E has Eigenvalue 1 with multiplicity $d$, then $(I-\overline{E})$ is bijective and hence invertible. In particular a solution of equations (8) exists and it is unique.

The proof can be found in the Appendix.

6 Calculating the estimator

There are two major ways of calculating the proposed estimator. Either one directly calculates $(I-\overline{E})^{-1}$ and applies it on $\widehat{\theta}$ or something closer to an iterative procedure. For the latter, by iterative application of (8) we derive that

[TABLE]

If the absolute values of the eigenvalues of $\overline{E}$ are bounded from above by a constant strictly smaller than 1, then (10) is well defined with $E^{\infty}=0$, and the converging series

[TABLE]

so that the iterative algorithm

[TABLE]

converges from any starting point. Note that (11) is the usual way the backfiting equations (6),(7) or equivalently (8) are solved. Another way is to calculate the finite sum

[TABLE]

with some stopping criteria $\overline{r}$. We conclude that the proposed estimator can be calculated in a straight forward manner from the compound Aalen estimator $\widehat{m}$ and the operator $\overline{E}.$

We now briefly discuss how $\overline{E}$ can be calculated in the simple case $1=d=p=q$. Here, $\overline{E}$ can be approximated by a $j\times k$ matrix where $j,k$ are the number of grid points in $[0,t_{max}]$ and $[a_{0},a_{max}]$, respectively. This is done by first calculating the values $E_{1}(s_{0},a_{0})$ and $E_{2}(a_{0},s_{0})$ for every grid point; see Remark 1 for the definitions of the the functions. We call the resulting matrices $E_{1}^{mx}$ and $E_{2}^{mx}$. Afterwards, $\overline{E}_{2}^{mx}$ is derived from $E_{2}^{mx}$, via

[TABLE]

The matrices are then transformed to the desired operator via

[TABLE]

Finally,

[TABLE]

So that given a function $h:[0,t_{max}]\times[a_{0},a_{max}]\rightarrow\Re$, one calculates its values on the grid and summarises it in a vector $h^{grid}$. The function $\overline{E}h$ is then approximated via $\overline{E}^{op}h^{grid}$ where the latter is a simple matrix multiplication.

7 Asymptotics

Note that we have

[TABLE]

where $m$ arises from $\widehat{m}$ by replacing $N$ by $\Lambda$. It is hereby quite remarkable that $\overline{E}$ is the observable operator from the previous sections and not some asymptotic limit. We further conclude that the least square solution (6) and (7) is a plug-in estimator of (12). The estimation error is then given as

[TABLE]

As in the last section, If $\overline{E}$ has eigenvalues all bounded away from one, then

[TABLE]

So the asymptotic behaviour of $\widehat{\theta}-\theta$ can be deduced from the asymptotic behaviour of $(I-\overline{E})^{-1}$ and $(\widehat{m}-m)$, with the latter being the compound estimation error of two additive Aalen models on different time-scales.

Theorem 1

Under assumptions (A)–(G), the estimator $\widehat{\theta}$ exists. Furthermore the estimator $\widehat{\theta}$ is $n^{1/2}$ consistent:

[TABLE]

in Skorohod space $D^{p+q}[0,a_{max}]$. Here, $(\widehat{\theta}-\theta)$ is treated as one stochastic process defined on $[0,a_{max}]$ by setting for $j=1,\dots,p$ and $\nu\in[t_{max},a_{max}]$, $(\widehat{\theta}-\theta)_{j}(\nu)=(\widehat{\theta}-\theta)_{j}(t_{max})$. And similarly, for $j=p+1,\dots,p+q$ and $\nu\in[0,a_{0}]$, $(\widehat{\theta}-\theta)_{j}(\nu)=0$. The process $U$ is a $p+q$ dimensional mean-zero Gaussian process with covariation matrix $\Sigma(\nu_{1},\nu_{2})$ described in the Appendix, and $\widetilde{E}$ is the limit of $\overline{E}$.

The proof can be found in the Appendix.

8 Confidence Bands

While we could use the central limit theorem of the previous section to construct confidence bands, it has been suggested that better small sample performance can be achieved by directly bootstrapping the estimation error. We propose a wild bootstrap approach based on the relationship

[TABLE]

Since ($I-\overline{E})^{-1}$ is known, it is enough to to only approximate $\mathcal{M}$. We do this via the wild bootstrap version

[TABLE]

or

[TABLE]

where $G_{i}$ is a mean zero random variable with unit variance. The random variable $G_{i}$ is generated such that for fixed $i$, it is independent to all other variables. It is straight forward to confirm that $\widehat{\mathcal{M}}^{(r)},\ r=1,2$ is a mean zero process that has the same covariance as $\mathcal{M}$ (The covariance of $\mathcal{M}$ is given in the appendix). Hence, we directly derive the following proposition.

Proposition 3

Under assumptions (A)–(G), the bootstrapped estimation error is uniformly consistent, i.e., for $r=1,2$

[TABLE]

in Skorohod space $D^{p+q}[0,T]$, where $U$ is is described in Theorem 1.

The proof can be found in the Appendix.

One useful consequence of this is that we can estimate standard errors of our estimator $\hat{\theta}$ based on the approximation from the bootstrap. We denote these estimators as $\hat{\sigma}_{r}(t)$ for the two components $r=1,2$.

Corollary 1

Under assumptions (A)–(G), the bootstrapped errors lead to confidence bands $CB^{(r)}$ for ${\theta}(\nu)$ over $\nu\in[\nu_{1},\nu_{2}]$ providing an asymptotic coverage probability of $1-\alpha$, where

[TABLE]

and

[TABLE]

We explore the performance of the estimator of the standard error and the uniform bands in the next section.

9 Simulations

We generated data from the simple two-time scale model with age and duration that resemble the data we consider in worked example in the next section. Thus assuming that the hazard for those under risk is given as $\beta(t+a_{i})+\alpha(t)$, where $\beta(a)\equiv 0.067$ and the entry ages where drawn uniformly from $[0,25]$ but making sure that 10 % of the data started in [math] to (to avoid difficulties with left truncation in the estimation). The $\alpha(t)$ component was piecewise constant with rate $0.32$ in the time-interval $[0,0.25]$, then $0.48$ in $(0.25,0.5]$ and then finally to satisfy our constraint $-0.044$ in $(0.5,5]$, so that $\int_{0}^{5}\alpha(s)ds=0$. All subjects were censored after $5$ years of follow up.

In all simulations we used a discrete approximation based on a time-grid of either 100 points in both the age direction $[0,30]$ and on the duration time-scale $[0,5]$.

9.1 Bias of backfitting

We considered sample sizes 100, 200 and 400 and show the bias for the two-components in Table 1 based on 1000 realizations.

We note that the the backfitting algorithm is almost unbiased across all sample size and improves as the sample size increases. This is despite the fact that the simulated component in the time-direction really is quite wild.

9.2 Bootstrap uncertainty

Secondly, we demonstrate that our bootstrap seems to work well to describe the uncertainty of the estimates. We simulated data as before and based on 1000 realisations with 100 bootstrap’s based on $G_{i}dN_{i}$ we estimated: a) the point-wise standard error for the two-components; b) computed the pointwise coverage baed on these; c) and constructed uniform confidence bands, as described in Corollary 1, for the the two components and its coverage.

Table 2 around here

We note that the standard error is well estimated by the bootstrapped standard deviation across all sample sizes and for both components. In addition the pointwise coverage is close to the nominal 95 % level for the larger sample sizes. But even for $n=100$ the coverage is reasonable for most time-points for the two components.

Finally, we also considered the performance of the confidence bands based on our bootstrap approach.

Table 3 around here

When $n$ gets larger these bands are quite close to the nominal 95 % level, but for $n=100$ the asymptotics have not quite set in to make the entire band work well.

10 Application to the TRACE study

The TRACE study group (see e.g. Jensen et al. (1997) ) has collected information on more than 4000 consecutive patients with acute myocardial infarction (AMI) with the aim of studying the prognostic importance of various risk factors on mortality. We here consider a subset of 1878 of these patients that are available in the timereg R package. At the age of entry (age of diagnosis) the patients had various risk factors recorded, but we here just show the simple model with the effects of the two-time-scales age and duration. It is expected that the duration time-scale has a strong initial effect of dying that then disappears when patients survive the first period right after their AMI.

We then estimated the two-time-scale model $\alpha(t)+\beta(t+a_{i})$ under the identifiability condition that $\int_{0}^{5}\alpha(s)ds=0$. Restricting attention to patients more than 40 years of age, and within the first 5 duration years after the diagnosis.

First we estimate the mortality on the two time-scales separately, the two marginal estimates, see Figure 1. Panel (a) shows the cumulative hazard on the age time-scale with the marginal estimate (full line) and the one with adjustment for duration effects (broken line), and panel (b) the mortality on the duration time-scale with the marginal estimate (full line) and with adjustment for age effects (broken line). We note that on the duration time-scale the cumulative hazard is quite steep. In addition we show 95 % confidence bands based on our bootstrap (regions), and the pointwise confidence intervals (dotted line).

Figure 1 about here

Taking out the duration effect slightly alters the estimate of the age-effect. In contrast the duration effect is strongly confounded by age effect estimates, and here the two-time scale model more clearly demonstrates what is going on on the duration time-scale. The duration effect is strong initially and then after surviving the first 220 days we see a protective effect (dotted vertical line).

We stress that the interpretation of the hazards on the two-time scales are difficult, due to, for example, the constraint that needs to be imposed to identify a specific solution. Nevertheless, it very useful to see the components from the two time-scales that jointly make up the hazard for an individual, and can be used for the prediction purposes as we demonstrate further below. Note also that due to the additive structure the duration effect can be interpreted as giving relative survival due to the duration time-scale.

Figure 2 about here

In Figure 2 we show the survival predictions for subjects that are 60, 70, or 80, respectively, using the two-time scale model. Thus computing $\exp(-(\hat{B}(a_{0}+t)-\hat{B}(a_{0}))+\hat{A}(t))$ and constructing the confidence bands using the bootstrap approach for $(\hat{B}(a_{0}+t)-\hat{B}(a_{0}))+\hat{A}(t)$ for $t\in[0,5]$. These curves are a direct consequence of having the two-components and are directly interpretable.

11 Discussion

By utilising the additive structure we have demonstrated that one can estimate the effect of two time-scales directly by a backfitting algorithm that does not involve smoothing. By working on the cumulative this also lead to uniform asymptotic description and a simple bootstrap procedure for getting estimates of the uncertainty and for constructing for example confidence intervals. These cumulative may form the basis for smoothing based estimates when the hazard are of interest, but often the cumulative are the quantities of key interest for example when interest is on survival predictions.

Clearly, the model could also be fitted by a more standard backfitting approach working on the hazard scale as in … for multiplicative hazard models.

Our backfitting approach can be extended for example the age-period-cohort model but here identifiability conditions are more complex to build into the estimation.

Appendix A Proofs

A.1 Proof of Proposition 1

With $f_{k}(s,u)=\left(0,\cdots,0,c_{k}s,0,\cdots,0,-c_{k}(u-a_{0}),0,\cdots,0)\right)^{T},\ (k=1,\dots,d)$, the proposition directly follows from $Lin(f_{1},\dots f_{d})\subseteq Kern(I-E)$, and the fact that $\widetilde{\Pi}$ is a projection into $Lin(f_{1},\dots f_{d})$.

A.2 Proof of Proposition 2

Since the eigenspace of $E_{2}$ corresponding to the eigenvalue equal 1 has dimension $d$, we know its exact form:

[TABLE]

see also Remark 2. One can then verify that $\textrm{kern}(I-E_{2})=\textrm{kern}(I-E_{2})^{l},l=2,3,\dots$. This is because linear functions cannot be constructed as sum of a linear and non-linear functions. Noting that $E_{2}$ is a compact operator, we conclude that $I-E_{2}$ is an isomorphism from $\textrm{Im}(I-E_{2})$ to $\textrm{Im}(I-E_{2})$. We introduce the operator $\overline{E}_{2}=E_{2}(I-\Pi)$, where $\Pi$ is the projection onto $\textrm{ker}(I-E_{2})$. The condition $\overline{E}(h_{1},h_{2})=(h_{1},h_{2})$ is equivalent to $\overline{E}_{2}E_{1}h_{2}=h_{2}$ and $E_{1}\overline{E}_{2}h_{1}=h_{1}$. Since the eigenspace of $E$ corresponding to an eigenvalue of 1 has dimension $d$, $\overline{E}(h_{1},h_{2})=(h_{1},h_{2})$ is not true for non-linear $h_{1},h_{2}$. This is because $\overline{E}=E$ when restricted on non-linear functions $h_{1},h_{2}$. When considering a linear $h_{1}$, then $\overline{E}_{2}h_{1}=0$. We conclude that the solution of $\overline{E}(h_{1},h_{2})=(h_{1},h_{2})$ is trivial. Hence the kern of $(I-\overline{E})$ is trivial. Since $\overline{E}$ is compact this means $(I-\overline{E})$ is bijective, in particular invertible.

A.3 Assumptions

We first define a few quantities.

For every $\nu$ in, $[0,a_{max}]$, we define the following matrices

[TABLE]

as well as

[TABLE]

We further define

[TABLE]

The limiting operator $\widetilde{E}$ is then defined analogue to $\overline{E}$ by replacing $E_{1}(s|u),E_{2}(u|s)$ by $\widetilde{E}_{1}(s|u),\widetilde{E}_{2}(u|s)$.

We make the following assumptions.

1. (A)

There exist continuous functions $r^{(1)}_{j},r^{(2)}_{jk},v^{(2)}_{jk},v^{(3)}_{jkl}$, $(j,k,l=1,\dots,p+q)$, such that for $n\rightarrow\infty$

[TABLE] 2. (B)

For $j=0,\dots,p$ and $k=0,\dots,q$, and $n\rightarrow\infty$

[TABLE] 3. (C)

For every $\nu$, the matrix $(r^{(2)}_{ij}(\nu))$is non-singular. 4. (D)

$\sup_{\left\lVert m\right\rVert_{2}\leq 1}\left\lVert(I-\widetilde{E})^{-1}m\right\rVert_{\infty}<\infty$, 5. (E)

The random variables $(a_{i})_{i=1,\dots,n}$ are $iid$, independent of (X,Z) and are absolutely continuous with continuous density $h$.

A.4 Proof of Theorem 1

We first prove the central limit theorem:

[TABLE]

We write

[TABLE]

Since $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are square integrable martingales (each with respect to its natural filtration), $\mathcal{M}$ is tight under the condition that its jumps are uniformly bounded. This follows from assumption (B), so $\mathcal{M}$ is indeed tight. Furthermore, under assumption (A), $\mathcal{M}$ is asymptotically uniformly close to

[TABLE]

which is the sum of $n$ $iid$ random processes. So the limit of $\mathcal{M}$, if it exists, must be Gaussian. Hence convergence of $\mathcal{M}$ to $U$ is verified by establishing point-wise convergence of the covariance matrix of $\mathcal{M}$ to the covariance matrix of $U$. For two points $\nu_{1},\nu_{2}$ in $[0,a_{max}]$ with $\nu_{1}\leq\nu_{2}$, $\textrm{Cov}\left(\mathcal{M}(\nu_{1}),\mathcal{M}(\nu_{2})\right)$ is a $(p+q)\times(p+q)$ matrix. We have

[TABLE]

With entry $(j,k)$ given by

[TABLE]

where

[TABLE]

The two processes in the covariance are running in in the same time-interval. This is because we could eliminate the non-intersecting time points due to independence. Under assumption (B), the entries converge to

[TABLE]

so that the two process in the covariance are now even martingales with respect to the same filtration $\mathcal{F}_{i}(\nu_{0})=\sigma\{V_{i}(u),N_{i}(u),u\leq\nu_{0}\}$. We can hence first calculate the conditional covariance, given $\mathcal{F}_{i}$, i.e., the predictable covariation process. Afterwards, the covariance is given as the expectation of predictable covariation process. For the predictable covariation process we get

[TABLE]

From assumptions (A),(C), (E) we conclude that $\textrm{Cov}\left(\mathcal{M}(\nu_{1}),\mathcal{M}(\nu_{2})\right)\rightarrow\Sigma(\nu_{1},\nu_{2})$ with entries

[TABLE]

Since the integral is well defined, we conclude convergence of $\mathcal{M}$ to $U$.

We now need to handle the operator $\overline{E}$. We have

[TABLE]

Equation (14) follows directly from the uniform convergence of the kernel functions $E_{1}(s|u),E_{2}(u|s)$ to $\widetilde{E}_{1}(s|u),\widetilde{E}_{2}(u|s)$ which is ensured via Assumptions (A)-(C), (E). Inequality (15) is ensured, since the kernel functions are bounded using the same assumptions. Together with Assumption (D) it follows that the operator $(I-\overline{E})^{-1}$ converges to the linear and bounded operator $(I-\widetilde{E})^{-1}$ which gives the desired central limit theorem.

A.5 Poof of proposition 3

For two points $\nu_{1},\nu_{2}$ in $[0,a_{max}]$ with $\nu_{1}\leq\nu_{2}$, the covariance of $\widehat{\mathcal{M}}^{(1)}$ is given by

[TABLE]

Under assumption (B) this is uniformly close to

[TABLE]

As in the proof of Theorem 1, the two processes in the covariance are martingales with respect to $\mathcal{F}_{i}$, so we can calculate the covariance as expectation of the predictable covariation process. Hence, $\textrm{Cov}\left(\widehat{\mathcal{M}}^{(1)}(\nu_{1}),\widehat{\mathcal{M}}^{(1)}(\nu_{2})\right)\rightarrow\Sigma(\nu_{1},\nu_{2})$.