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

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

Munir Hiabu munir.hiabu@sydney.edu.au School of Mathematics and Statistics, University of Sydney, Camperdown NSW 2006, Australia Jens P. Nielsen Jens.Nielsen.1@city.ac.uk Cass Business School, City, University of London, 106 Bunhill Row, London, EC1Y 8TZ, United Kingdom Thomas H. Scheike thsc@sund.ku.dk (Corresponding author) Department of Public Health, University of Copenhagen, Øster Farimagsgade 5B, 1014 Copenhagen K, Denmark
###### 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.

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

\NewDocumentCommand\normL

s O m \IfBooleanTF#1∥#3∥#2∥#3#2∥_L_2(Ω)

## 1 Introduction

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 be 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

 λi(t) =p∑j=1Xij(t)αj(t)+q∑k=1Zik(t)βk(t+ai) =Xi(t)α(t)+Zi(t)β(t+ai),(0≤t≤tmax), (1)

where and are tupels of one dimensional deterministic functions, and are predictable cadlag covariate vectors with and having almost surely full rank, and is a real-valued random variable observed at time . If for all , 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 and where the latter time-scale is specific to each individual and 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, might be time since diagnosis (duration) among subjects that have entered the illness stage of the model and could be the age when the transition to the illness stage occurred, such that is the age of the subject.

After introducing some notation we present an estimation procedure that leads to explicit estimators of and . The cumulative effects have the advantage compared to and 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 one may estimate the intensity by smoothing techniques.

2.1  Notation Let such that are martingales. Let further be the n-dimensional counting process, is its compensator, such that is an n-dimensional martingale, and define matrices and with dimensions and , respectively. The individual entry times are summarised in one vector . A superscript denotes a shift in the argument, i.e, for a generic function , . For a generic matrix , with rows , and a n-dimensional vector , is defined through shifting the rows: . For a generic matrix , a minus superscript, , denotes the Moore-Penrose inverse. An integral, , 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 and the design. If this is the case, then the functions and 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 and share the first columns, i.e., for all ,

 Xil=Zil,l≤d.

We formulate the problem using group-theoretic arguments, see also Carstensen (2007); Kuang et al. (2008). Fix constants and define as valued function having all entries but the and the equal zero:

 fl(s,u)=(0,⋯,0,cls,0,⋯,0,−cl(u−a0),0,⋯,0))T, (l=1,…,d).

We define the group by

 G={g:(AB)↦(AB)+h |h∈Lin(f1,…fd)}.

The identification problem can be rephrased as that the intensity defined in (2) is a function of , which is invariant to transformations . In the sequel we circumvent the identification issue by adding the following constraint

 Al(tmax)=∫tmax0αl(s)ds=0,(l=1,…,d), (2)

noting that for any solution of model (2), there exists a unique solution 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 via the following least squares criteria.

 argmin¯¯¯¯A,¯¯¯¯B∑i∫{∫t0dNi(s)−∑j∫t0Xij(s)d¯¯¯¯Aj(s)−∑k∫t0Zik(s)d¯¯¯¯Baik(s)}2dt,

where the integrals can be understood as Stieltjes integrals, noting that and 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 is a step function. To simplify notation we will generally work in matrix notation so that above minimisation criteria can also be written as

 argmin¯¯¯¯A,¯¯¯¯B∑i∫{∫t0dNi(s)−∫t0Xi(s)d¯¯¯¯A(s)−∫t0Zi(s)d¯¯¯¯Bai(s)}2dt.

Straight forward computations utilzing calculus of variations lead to solving the following first order conditions for all , :

 ∑iXi(t)T{dNi(t)−Xi(t)dˆA(t)−Zi(t)dˆBai(t)dt}=0, ∑iZ−aii(a)T{dN−aii(a)−Z−aii(a)dˆB(a)−X−aii(a)dˆA−ai(a) }=0.

Rearranging yields

 ∑iZ−aii(a)TdN−aii(a)−∑iZ−aii(a)TX−aii(a)dˆA−ai(a)=Z−a∙(a)TZ−a∙(a)dˆB(a).

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

 ˆA(t) =∫t0X(s)−dN(s)−∫E1(t|u)dˆB(u) (3) ˆB(a) =∫aa0Z−a∙(u)−dN−a∙(u)−∫E2(a|s)dˆA(s), (4)

where

 E1(s|u) =∑i{XT(u−ai)X(u−ai)}−1X−ai,Ti(u)Z−aii(u)I(ai≤u≤ai+s), E2(u|s) =∑i{Z−a∙,T(s+ai)Z−a∙(s+ai)}−1ZTi(s)Xi(s)I(a0−ai≤s≤u−ai).
###### Remark 1

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

 λi(t)=Yi(t){α(t)+β(ai+t)},

with , the risk indicators are

 E1(s|u) =∑i1∑i′Yi′(u−ai)Y−aii(u)I(ai≤u≤ai+s), E2(u|s) =∑i1∑i′Y−ai′i′(s+ai)Yi(s)I(a0−ai≤s≤u−ai).

## 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:

 (ˆAˆB)=(∫t0X(s)−dN(s)∫aa0Z−a∙(u)−dN−a∙(u))+(0−E1−E20)×(ˆAˆB),

where with some miss-use of notation . Or even simpler

 ˆθ=ˆm+Eˆθ, (5)

with obvious notation, and linear operator :

 ˆθ=(ˆAˆB),ˆm=(∫t0X(s)−dN(s)∫aa0Z−a∙(u)−dN−a∙(u)),E=(0−E1−E20).

Note that is composed of the marginal Aalen estimators of the two time scales, and . Additionally, the operator is compact because it is the composition of an integral operator, which is compact, and a derivative operator, which is bounded. The operator 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 are bounded away from one, then, is invertible and we have

 ˆθ=(I−E)−1ˆm.

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

###### Remark 2

Consider the most simple case , i.e., . Given a constant , consider the pair of linear function with , as defined in Section 3. Assuming that and are bounded away from zero on the whole range , one can easily verify that

 E2f11(u) =c∫E2(u|s)ds=c(u−a0), E1f12(s) =−c∫E1(s|u)du=−cs.

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

 ∫E1(s|u)du=∑i∫ai+sai1∑i′Yi′(u−ai)Y−aii(u)du=∫s0∑iYi(t)∑i′Yi′(t)dt=s.

Hence, we have

 E(f11f12)=(−E1f12−E2f12)=(f11f12).

So that one is clearly an eigenvalue of with corresponding eigenfunction . 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 . Functions are defined in Section 3.

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

 ˆA(t) =∫t0X(s)−dN(s)−∫E1(t|u)dˆB(u), (6) ˆB(a) =∫aa0Z−a∙(u)−dN−a∙(u)−∫E2(a|s)dˆA(s)+ˆAdq(tmax)tmax(a−a0), (7)

where is the q-dimensional vector . This translates to the new operator equation

 ˆθ=ˆm+¯¯¯¯Eˆθ,¯¯¯¯E=(0−E1−¯¯¯¯E20), (8)

where . The next proposition states that the solutions of include all relevant solutions of (5) and that every solution of is a solution of .

###### Proposition 1

For every solution of (5), define

 ˆθ0=(I−˜Π)ˆθ,

where

 ˜Π(h1(t)h2(a))=(thdp1(tmax)t−1max−(a−a0)hdq1(tmax)t−1max).

Then is a solution of and

 ˆθ0+Lin(f1,…fd), (9)

are further solutions of (5). Reversly, for every solution 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 is known and hence one can calculate a numerical approximation of its eigenvalues by working on a grid. Consider the sub-space

 K={h=(h1,…,hd,0,…,0)| hl:R→R, x↦clx, cl∈R, l=1,…,d}.

It holds that , where is a projection into . We have . We can check whether equals . This can be done by calculating the dimension of the eigenspace of corresponding to an eigenvalue equal one. The dimension will be at least . If it is exactly , then .

The next proposition states that if , and , then both and are bijective.

###### Proposition 2

Assume that has Eigenvalue 1 with multiplicity . Then, will be bijective. If furthermore E has Eigenvalue 1 with multiplicity , then 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 and applies it on or something closer to an iterative procedure. For the latter, by iterative application of (8) we derive that

 ˆθ=∞∑r=0¯¯¯¯Er(ˆm)+¯¯¯¯E∞(ˆθ). (10)

If the absolute values of the eigenvalues of are bounded from above by a constant strictly smaller than 1, then (10) is well defined with , and the converging series

 ˆθ=∞∑r=0¯¯¯¯Er(ˆm),

so that the iterative algorithm

 ˆθ(r)=ˆm+¯¯¯¯Eˆθ(r−1) (11)

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

 ˜θ=¯r∑r=0¯¯¯¯Er(ˆm),

with some stopping criteria . We conclude that the proposed estimator can be calculated in a straight forward manner from the compound Aalen estimator and the operator

We now briefly discuss how can be calculated in the simple case . Here, can be approximated by a matrix where are the number of grid points in and , respectively. This is done by first calculating the values and for every grid point; see Remark 1 for the definitions of the the functions. We call the resulting matrices and . Afterwards, is derived from , via

 ¯¯¯¯Emx2=Emx2+⎛⎜ ⎜ ⎜ ⎜ ⎜⎝0⋯0s1/sj0…0s2/sj⋮⋮⋮0…01⎞⎟ ⎟ ⎟ ⎟ ⎟⎠.

The matrices are then transformed to the desired operator via

 Δ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1−10⋯00⋱⋱⋱⋮⋮⋱⋱⋱0⋮⋱⋱⋱−10⋯⋯01⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,Eop1=Emx1×Δ,¯¯¯¯Eop2=¯¯¯¯Emx2×Δ.

Finally,

 ¯¯¯¯Eop=(0−Eop1−¯¯¯¯Eop20).

So that given a function , one calculates its values on the grid and summarises it in a vector . The function is then approximated via where the latter is a simple matrix multiplication.

## 7 Asymptotics

Note that we have

 θ=m+¯¯¯¯Eθ, (12)

where arises from by replacing by . It is hereby quite remarkable that 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

 ˆθ−θ=ˆm−m+¯¯¯¯E(ˆθ−θ). (13)

As in the last section, If has eigenvalues all bounded away from one, then

 ˆθ−θ=(I−¯¯¯¯E)−1(ˆm−m).

So the asymptotic behaviour of can be deduced from the asymptotic behaviour of and , 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 exists. Furthermore the estimator is consistent:

 n−1/2(ˆθ−θ)→(I−˜E)−1U,

in Skorohod space . Here, is treated as one stochastic process defined on by setting for and , . And similarly, for and , . The process is a dimensional mean-zero Gaussian process with covariation matrix described in the Appendix, and is the limit of .

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

 ˆθ−θ=(I−¯¯¯¯E)−1(ˆm−m) =(I−¯¯¯¯E)−1(∫t0X(s)−dM(s)∫aa0Z−a∙(u)−dM−a∙(u)) =(I−¯¯¯¯E)−1(M1M2)

Since ( is known, it is enough to to only approximate . We do this via the wild bootstrap version

or

 ˆM(2) ∫t0˜Mi(s)ds =Gi(∫t0Ni(s)ds−(∫t0(Xi(s)dˆA(s)+∫t0Zi(s)dˆB(s+ai))),

where is a mean zero random variable with unit variance. The random variable is generated such that for fixed , it is independent to all other variables. It is straight forward to confirm that is a mean zero process that has the same covariance as (The covariance of 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

 n−1/2((I−¯¯¯¯E)−1ˆM(r))→(I−˜E)−1U,

in Skorohod space , where 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 based on the approximation from the bootstrap. We denote these estimators as for the two components .

###### Corollary 1

Under assumptions (A)–(G), the bootstrapped errors lead to confidence bands for over providing an asymptotic coverage probability of , where

 CB(r)(ν)=θ(ν)+/−c1−α^σr(ν),

and

 c1−α=(1−α)quantile of L⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩sup[ν1,ν2]n−1/2∣∣∣(I−¯¯¯¯E)−1ˆM(r)∣∣∣^σr|X,Z,N⎫⎪ ⎪ ⎪⎬⎪ ⎪ ⎪⎭

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 , where and the entry ages where drawn uniformly from but making sure that 10 % of the data started in to (to avoid difficulties with left truncation in the estimation). The component was piecewise constant with rate in the time-interval , then in and then finally to satisfy our constraint in , so that . All subjects were censored after 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 and on the duration time-scale .

### 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 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 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 gets larger these bands are quite close to the nominal 95 % level, but for 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 under the identifiability condition that . 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).

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.

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 and constructing the confidence bands using the bootstrap approach for for . 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 , the proposition directly follows from , and the fact that is a projection into .

### a.2 Proof of Proposition 2

Since the eigenspace of corresponding to the eigenvalue equal 1 has dimension , we know its exact form:

 {h=(h1,…,hd,0,…,0)| hl is linear, l=1,…,d};

see also Remark 2. One can then verify that . This is because linear functions cannot be constructed as sum of a linear and non-linear functions. Noting that is a compact operator, we conclude that is an isomorphism from to . We introduce the operator , where is the projection onto . The condition is equivalent to and . Since the eigenspace of corresponding to an eigenvalue of 1 has dimension , is not true for non-linear . This is because when restricted on non-linear functions . When considering a linear , then . We conclude that the solution of is trivial. Hence the kern of is trivial. Since is compact this means is bijective, in particular invertible.

### a.3 Assumptions

We first define a few quantities.
For every in, , we define the following matrices

 R(ν)=(X(ν)00Z−a∙(ν)),V(ν)=(X(ν),Z(ν)),

as well as

 R(1)j(ν) =∑iRij(ν) R(2)jk(ν) =∑iRij(ν)Rik(ν), V(2)jk(ν)