Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC5116003

Formats

Article sections

- Summary
- 1. Introduction
- 2. Proposed Methods
- 3. Asymptotic Properties
- 4. Simulations
- 5. Application to Liver Transplant Data
- 6. Discussion
- Supplementary Material
- References

Authors

Related links

Biometrics. Author manuscript; available in PMC 2017 March 24.

Published in final edited form as:

PMCID: PMC5116003

NIHMSID: NIHMS794870

The publisher's final edited version of this article is available at Biometrics

Treatments are frequently evaluated in terms of their effect on patient survival. In settings where randomization of treatment is not feasible, observational data are employed, necessitating correction for covariate imbalances. Treatments are usually compared using a hazard ratio. Most existing methods which quantify the treatment effect through the survival function are applicable to treatments assigned at time 0. In the data structure of our interest, subjects typically begin follow-up untreated; time-until-treatment and the pre-treatment death hazard are both heavily influenced by longitudinal covariates; and subjects may experience periods of treatment ineligibility. We propose semiparametric methods for estimating the average difference in restricted mean survival time attributable to a time-dependent treatment, the average effect of treatment among the treated, under current treatment assignment patterns. The pre- and post-treatment models are partly conditional, in that they use the covariate history up to the time of treatment. The pre-treatment model is estimated through recently developed landmark analysis methods. For each treated patient, fitted pre- and post-treatment survival curves are projected out, then averaged in a manner which accounts for the censoring of treatment times. Asymptotic properties are derived and evaluated through simulation. The proposed methods are applied to liver transplant data in order to estimate the effect of liver transplantation on survival among transplant recipients under current practice patterns.

It is often of interest in biomedical settings to evaluate the benefit of a treatment on survival. In many clinical settings, treatment is not administered right at the time of diagnosis, such that a period of waiting time occurs for some (or perhaps all) patients. In cases where treatment is not randomized, it is often useful to assess the benefit of treatment under current treatment assignment patterns. Through the average effect-of-treatment-on-the-treated (ETT; Pearl, 2009), one can evaluate the benefit of treatment as currently practiced.

Survival probabilities are easily understood by health care professionals, as is the area under the survival curve (restricted mean lifetime). Various authors have proposed using Cox regression with the primary goal not being to estimate hazard ratios, but to compare differences in survival and/or restricted mean lifetime. For example, Zucker (1998) and Chen and Tsiatis (2001) proposed methods that involved averaging over fitted values from Cox models. Zhang and Schaubel (2011) extended the methods of Chen and Tsiatis (2001) to accommodate dependent censoring, then subsequently developed double-robust methods (Zhang and Schaubel, 2012). Each of the afore-described methods applies to treatments assigned at baseline, as opposed to time-varying treatments.

In the data structure of interest in this report, all patients begin follow-up untreated, with some eventually receiving treatment and others dying beforehand. Pre-treatment mortality and treatment assignment rates are dependent on longitudinal covariates (including periods during which a subject is declared treatment-ineligible), such that a patient’s pre-treatment death is dependently censored by the receipt of treatment. Post-treatment survival is dependent on a subject’s condition at the time of treatment, and the duration of pre-treatment follow-up time. Our objective is to contrast two scenarios: (a) treatment is never assigned (b) treatment is assigned according to current practice patterns.

The proposed methods are motivated by the end-stage liver disease (ESLD) setting. The number of available deceased-donor livers is always less than the number of patients in need of liver transplantation. As a result, medically suitable patients are placed on a liver transplant waiting list. Patients typically begin follow-up on the wait list (‘untreated’; i.e., not transplanted), such that transplantation can be viewed as a time-dependent treatment. In the U.S., chronic end-stage liver disease patients are sequenced in decreasing order of Model for End-Stage Liver Disease (MELD) score, a very strong predictor of pre-transplant mortality. Transplantation results in the dependent censoring of pre-transplant death, since MELD scores predict both wait list mortality and transplant rates. Note that patients may be removed from the wait list (or made inactive) and, in such cases, are permanently (or temporarily) ineligible to receive a transplant. In the setting of our interest, the effect of treatment on the treated is of greater interest than the average causal effect, due to the implausibility of all patients receiving treatment.

Our analysis in Section 5 is different from that in Gong and Schaubel (2013) since (i) the former only looked at pre-transplant survival; (ii) did not compare post- versus pre-transplant survival; (iii) reported contrasts only in terms of the hazard ratio; and (iv) did not exclude Status 1 (acute liver failure) patients and, in fact, focused on contrasting them with chronic ESLD patients.

We develop semiparametric methods to estimate the average effect-of-treatment-on-the-treated through partly conditional modeling. The proposed method involves averaging over the observed instances of treatment initiation, with the averaging accounting for the various complexities in data structure; e.g., treatment initiation times are subject to right censoring; patients may die before treatment is received; and patients cannot initiate treatment while ineligible. For each treated patient, we use the accrued history (up to the time of treatment initiation) to project out a survival curve for post-treatment residual lifetime. Based on the same accrued pre-treatment history, we also project out the survival curve that would apply in the absence of treatment. This set-up lends itself well to partly conditional modeling (Zheng and Heagerty, 2005; Gong and Schaubel, 2013); see also the closely related concept of landmark analysis (Feuer et al., 1992; van Houwelingen, 2007; van Houwelingen and Putter, 2012; Parast, Tian and Cai, 2014). Gong and Schaubel (2013) developed methods for fitting partly conditional hazard regression models which apply to the absence-of-treatment setting in our set-up. We extend the ideas in Gong and Schaubel (2013) to estimate the average ETT through residual survival and restricted mean survival time. Although we focus on partly conditional modeling in this report, it should be noted that other pertinent methods exist, as described in Section 6.

The remainder of this article is organized as follows. In Section 2, we describe the proposed methods. Asymptotic properties are provided in Section 3 (for proofs, see Supplementary Materials), with finite-sample properties evaluated through simulation in Section 4. We apply the proposed methods to the motivating data set in Section 5. Concluding remarks are made in Section 6.

We now formalize the ideas introduced in Section 1, in the absence of censoring. We remove subscripting, such that defined variates pertain to any hypothetical subject. We let *T* represent treatment time, with *T* > 0 since subjects begin follow-up untreated. Death time in the absence of treatment is denoted by *D*^{0}. Note that, consistent with the intent-to-treat principle, patients that initiate treatment are considered to be ‘treated’ thereafter. Let (*s*) = 1 if the patient is treatment-eligible (i.e., eligible to initiate treatment) at time *s*, and 0 otherwise. A patient may oscillate between the eligible and ineligible states before time *D*^{0} ∧ *T*, where *a* ∧ *b* = min(*a, b*). In particular, (*s*) = 0 for *s* > *D*^{0} ∧ *T*, since a patient cannot *initiate* treatment more than once, and cannot initiate treatment after death. Note that a patient may only initiate treatment while eligible; i.e., *dI*(*T* ≤ *s*) = (*s*)*dI*(*T* ≤ *s*).

Under the above-listed Scenario (a), *T* = ∞. Under Scenario (b), treatment only occurs when *T* < *D*^{0}, in which case *D*^{0} is considered latent; *D*^{0} serves as a competing risk for *T*. For a patient with treatment time *T* = *s, D*^{1} is the death time, such that (*D*^{1} − *s*)_{+} is the residual post-treatment survival, with *a*_{+} = *a* · *I*(*a* > 0) and *I*(·) being the familiar 0/1 indicator function. The quantity (*D*^{0} − *s*)_{+} then represents the residual survival beyond *s* that would have been observed in the absence of treatment. Note that if *D*^{0} < *T*, then *D*^{1} is undefined.

The covariate vector, which contains some time-varying elements, is denoted by * Z**(

For a patient initiating treatment at time *T* = *s*, there are two death times of interest; the post-treatment residual death time, (*D*^{1} − *s*)_{+}, and the residual death time that would have occurred in the absence of treatment, (*D*^{0} − *s*)_{+}. At the time of treatment (e.g., *T* = *s*), we observe (*s*), and (*s*) = 1. Conditional on [(*s*), *T* = *s*], we contrast

$${S}_{1}(t;s|\mathit{\mathscr{H}}(s),T=s)=P\{({D}^{1}-s)>t|\mathit{\mathscr{H}}(s),T=s)$$

(1)

$${S}_{0}(t;s|\mathit{\mathscr{H}}(s),T=s)=P\{({D}^{0}-s)>t|\mathit{\mathscr{H}}(s),T=s)$$

(2)

the survival functions pertaining to the counterfactual variates (*D*^{1} − *s*)_{+} and (*D*^{0} − *s*)_{+}, respectively. Note that, in both *S*_{1}(*t*; *s*|·) and *S*_{0}(*t*; *s*|·), the time index *s* represents conditioning time, while *t* refers to residual survival *t* time units beyond the conditioning time, *s*. That is, *S _{j}*(

For fixed *L* > 0, restricted mean residual survival times are given by

$${\mu}_{1}(L;s|\mathit{\mathscr{H}}(s),T=s)={\int}_{0}^{L}{S}_{1}(t;s|\mathit{\mathscr{H}}(s),T=s)dt$$

(3)

$${\mu}_{0}(L;s|\mathit{\mathscr{H}}(s),T=s)={\int}_{0}^{L}{S}_{0}(t;s|\mathit{\mathscr{H}}(s),T=s)dt.$$

(4)

Conditioning on [(*s*), *T* = *s*], a pertinent contrast in survival functions is then

$$\delta (t;s|\mathit{\mathscr{H}}(s),T=s)={S}_{1}(t;s|\mathit{\mathscr{H}}(s),T=s)-{S}_{0}(t;s|\mathit{\mathscr{H}}(s),T=s),$$

(5)

while a contrast in restricted mean residual lifetime is defined as

$$\mathrm{\Delta}(L;s|\mathit{\mathscr{H}}(s),T=s)={\mu}_{1}(L;s|\mathit{\mathscr{H}}(s),T=s)-{\mu}_{0}(L;s|\mathit{\mathscr{H}}(s),T=s).$$

(6)

Average survival functions are then defined as

$${S}_{1}(t)=E[{S}_{1}(t;T|\mathit{\mathscr{H}}(T),T)]$$

$${S}_{0}(t)=E[{S}_{0}(t;T|\mathit{\mathscr{H}}(T),T)],$$

(7)

where, in each case, the expectation is taken with respect to the joint distribution of [(*T*), *T*)] over the identifiable range of *T* which would in practice be capped by the maximum follow-up time. Correspondingly, average restricted mean residual lifetimes are:

$${\mu}_{1}(L)=E[{\mu}_{1}(L;T|\mathit{\mathscr{H}}(T),T)]={\int}_{0}^{L}{S}_{1}(t)dt$$

$${\mu}_{0}(L)=E[{\mu}_{0}(L;T|\mathit{\mathscr{H}}(T),T)]={\int}_{0}^{L}{S}_{0}(t)dt.$$

(8)

The ETT can then be defined in terms of mean difference in survival probability as

$$\delta (t)=E[\delta (t;T|\mathit{\mathscr{H}}(T),T)]={S}_{1}(t)-{S}_{0}(t)$$

(9)

and in terms of average difference in residual mean survival time, by

$$\mathrm{\Delta}(L)=E[\mathrm{\Delta}(L|\mathit{\mathscr{H}}(T),T)]={\mu}_{1}(L)-{\mu}_{0}(L)={\int}_{0}^{L}\delta (t)dt.$$

(10)

Having specified the treatment effect of interest, the remaining subsections in Section 2 describe the proposed methods for estimating δ(*t*) and Δ(*L*).

We let *D _{i}* denote the death time for subject

We now describe the assumed models for ${({D}_{i}^{1}-{T}_{i})}_{+},{({D}_{i}^{0}-{T}_{i})}_{+}$, *T _{i}* and

We let λ_{1}(*t*; *s*|(*s*), *T* = *s*) denote the hazard function corresponding to *S*_{1}(*t*; *s*|(*s*), *T* = *s*) from (1). We assume the following proportional hazards model,

$${\lambda}_{1}(t;s|{\mathit{\mathscr{H}}}_{i}(s),{T}_{i}=s)={\lambda}_{01}(t)\text{exp}\{{\mathit{\beta}}_{1}^{\prime}{\mathit{Z}}_{i1}(s)\},$$

(11)

where the covariate **Z**_{i1}(*s*) is chosen to summarize the pre-treatment history, {_{i}(*s*), *T _{i}* =

We assume that treatment times are independently censored by *C _{i}*. Assuming that (

We begin by describing the assumed hazard model for survival in the absence of treatment. We then outline the proposed data augmentation, which involves selecting calendar date cross-sections. Next, we detail fitting the model through an inverse weighted and stratified log rank estimating function.

We let λ_{0}(*t*; *s*|(*s*), *T* = *s*) denote the hazard function corresponding to (2). Under strong ignorability, note that λ_{0}(*t*; *s*|_{i}(*s*), *T _{i}* =

$${\lambda}_{0}(t;s|{\mathit{\mathscr{H}}}_{i}(s),{\mathcal{E}}_{i}(s)=1)={\lambda}_{00}(t)\text{exp}\{{\mathit{\beta}}_{0}^{\prime}{\mathit{Z}}_{i0}(s)\},$$

(12)

where **Z**_{i0}(*s*) is chosen such that λ_{0}(*t*; *s*|_{i}(*s*), _{i}(*s*) = 1) = λ_{0}(*t*; *s*|**Z**_{i0}(*s*)), implying that **Z**_{i0}(*s*) contains all elements of the history pertinent to predicting ${({D}_{i}^{0}-s)}_{+}$, including all appropriate functions of time-already-survived, *s*. Model (12) is partly conditional (Zheng and Heagerty, 2005; Gong and Schaubel, 2013) since, although the hazard at time *s* + *t* is of interest, the model conditions on information which is ‘frozen’ at time *s*. In contrast, a typical (fully) conditional or ‘time-dependent’ model would condition on _{i}(*s* + *t*).

The motivation for using a partly conditional model is described at the start of Section 2.4. Generally, fitting a partly conditional model requires some form of data augmentation in which the records corresponding to each subject’s observed data are restructured in order to facilitate fitting the assumed model. After such augmentation, each input record has a prior time survived (e.g., *s _{i}*) and corresponding prior history

In landmark analysis, typically survival from a specific follow-up time point (or set of specific time points) is desired, with survival probability projected out after the chosen landmark time(s). In our case, since treatment can occur at any time point (e.g., *T* = *s*), we need to be able to project conditional survival forward from any conditioning time *s*. This suggests a partly conditional model which includes terms representing previous time survived, *s*. Variation in previous time survived is then required, which means that sampling component of the data augmentation should be based on something other then *s* itself. We choose to sample based on calendar time, since each calendar time cross-section will contain wide variation in previous time survived. As we later describe, we stratify the model by cross-section for computational savings, which is important in large data sets like that we analyze in Section 5. For instance, Gong and Schaubel (2013) developed a partly conditional model which chooses the conditioning times to be the *s _{i}* values observed on a randomly selected calendar date. For example, consider a particular calendar date (e.g., 07/01/2004); input records for fitting the model would consist of

The estimation of **β**_{0} from model (12) was developed by Gong and Schaubel (2013). The essential ideas are presented here for continuity, and because the authors only derived the properties of _{0}, but not those of *Ŝ*_{0}(*t*; *s*|**Z**_{i0}(*s*)), _{0}(*L*; *s*|**Z**_{i0}(*s*)), *Ŝ*_{0}(*t*) or _{0}(*L*).

To begin, we choose a set of *K* calendar dates, {*CS*_{1}, …, *CS _{K}*}. Each cross-section date

We now establish additional notation pertinent to model (12). Since survival time from cross-section is modeled, we define the following times-since-cross-section: *D _{ik}* = (

Examples of the relationship between cross-section time and follow-up time. Four subjects (*i* = 1, …, *i* = 4) and two cross sections (*k* = 1, 2) are shown. The four subjects begin follow-up at different calendar dates. For subject *i* = 1, failure **...**

Following Gong and Schaubel (2013), we estimate **β**_{0} through the stratified model,

$${\lambda}_{0k}(t;s|{\mathit{\mathscr{H}}}_{i}({s}_{ik}),{\mathcal{E}}_{i}({s}_{ik})=1)={\lambda}_{00k}(t)\text{exp}\{{\mathit{\beta}}_{0}^{\prime}{\mathit{Z}}_{i0}({s}_{ik})\},$$

(13)

where **β**_{0} is the same parameter in the unstratified model of interest, (12). Model (13) is quite flexible. Non-proportionality can be accommodated by replacing **β**_{0} with **β**_{0}(*t*), a parametric function on *t*. The parameter vector could also be allowed to be a parametric function of previous time survived; e.g., **β**_{0k}, or **β**_{0}(*s _{ik}*). Moreover, interactions between

Model (13) conditions on _{i}(*s _{ik}*). However, we anticipate that

$${\lambda}_{i}^{T}(t|{\mathit{\mathscr{H}}}_{i}(t),{\mathcal{E}}_{i}(t))=\mathcal{E}({s}_{ik}){\mathcal{E}}_{i}(t){\lambda}_{0}^{T}(t)\text{exp}\{{\mathit{\theta}}_{0}^{\prime}{\mathit{Z}}_{i}(t)\},$$

(14)

$$\mathcal{E}({s}_{ik}{\lambda}_{ik}^{\u2020}(t;{s}_{ik}|{\mathit{Z}}_{i0}({s}_{ik}),\mathcal{E}({s}_{ik})=\mathcal{E}({s}_{ik}){\lambda}_{0k}^{\u2020}(t)\text{exp}\{{\mathit{\theta}}_{1}^{\prime}{\mathit{Z}}_{i}({s}_{ik}))\},$$

(15)

with model (14) assumed to be the correct model; model (15) is expected to be misspecified, but is only used to provide a weight stabilizer. We assume no-unmeasured-confounders for treatment, ${\lambda}_{i}^{T}(t|{\mathit{\mathscr{H}}}_{i}(t))={\lambda}_{i}^{T}(t|{\mathit{\mathscr{H}}}_{i}({D}_{i}),{D}_{i})$, and that ${\lambda}_{i}^{T}(t|{\mathit{\mathscr{H}}}_{i}(t))={\lambda}_{i}^{T}(t|{\mathit{Z}}_{i}(t))$. Note that ${\lambda}_{i}^{T}(t|{\mathit{\mathscr{H}}}_{i}(t),{\mathcal{E}}_{i}(t))$ in (14) uses (total) follow-up time *t* (measured from time 0) as the time axis, conditions on information on [0, *t*), while ${\lambda}_{ik}^{\u2020}(t;{s}_{ik}|{\mathit{Z}}_{i0}({s}_{ik}),\mathcal{E}({s}_{ik})=1))$ in (15) uses (residual) time since *s _{ik}* and conditions on the history over [0,

As derived in Gong and Schaubel (2013), an appropriate weight function is given by

$${W}_{ik}^{A}(t)={Y}_{i0k}(t)\text{exp}\{{\mathrm{\Lambda}}_{i}^{T}({s}_{ik}+t)-{\mathrm{\Lambda}}_{i}^{T}({s}_{ik})\},$$

(16)

where ${\mathrm{\Lambda}}_{i}^{T}(t)={\int}_{0}^{t}{\mathcal{E}}_{i}(u){\lambda}_{0}^{T}(u)\text{exp}\{{\mathit{\theta}}_{0}^{\prime}{\mathit{Z}}_{i}(u)\}du$. The quantity ${W}_{ik}^{A}(t)$ can be thought of as the inverse of the conditional probability of remaining untreated at time (*s _{ik}* +

$${W}_{ik}^{B}(t)={Y}_{i0k}(t)\frac{\text{exp}\{{\mathrm{\Lambda}}_{i}^{T}({s}_{ik}+t)-{\mathrm{\Lambda}}_{i}^{T}({s}_{ik})\}}{\text{exp}\{{\mathrm{\Lambda}}_{ik}^{\u2020}(t)\}}.$$

(17)

Note that artificially censoring subjects at *t* = *L* would be an alternative to the stabilizer.

An estimator for **β**_{0}, denoted by _{0}, is obtained through solving the following inverse-weighted score function,

$${\mathit{U}}_{0}(\mathit{\beta})=\sum _{i=1}^{n}\sum _{k=1}^{K}{\int}_{0}^{{\tau}_{0k}}{\mathcal{E}}_{i}({S}_{ik})\{{\mathit{Z}}_{i0}({s}_{ik})-{\overline{\mathit{Z}}}_{0k}(t;\mathit{\beta},W)\}{W}_{ik}^{B}(t)d{N}_{i0k}(t),$$

(18)

with ${\overline{\mathit{Z}}}_{0k}(t;{\mathit{\beta}}_{0})={\mathit{R}}_{0k}^{(1)}(t;{\mathit{\beta}}_{0})/{R}_{0k}^{(0)}(t;{\mathit{\beta}}_{0})$ and ${\mathit{R}}_{0k}^{(d)}(t;{\mathit{\beta}}_{0})={n}^{-1}{\sum}_{i=1}^{n}{\mathcal{E}}_{i}({s}_{ik}){W}_{ik}(t){\mathit{Z}}_{i0}{({s}_{ik})}^{\otimes d}\text{exp}\{{\mathit{\beta}}_{0}^{\prime}{\mathit{Z}}_{i0}({s}_{ik})\}$ with *d* = 0, 1, 2 and where τ_{0k} satisfies *P*{*Y*_{i0k}(τ_{0k}) = 1} > 0, and can in practice be set to the largest *X _{ik}* among subjects with

$${\hat{\mathrm{\Lambda}}}_{00}(t;{\hat{\mathit{\beta}}}_{0})={n}^{-1}\sum _{i=1}^{n}\sum _{k=1}^{K}{\int}_{0}^{t}{R}_{0}^{(0)}{(u;{\hat{\mathit{\beta}}}_{0})}^{-1}{\mathcal{E}}_{i}({s}_{ik}){W}_{ik}^{B}(u)d{N}_{i0k}(u)$$

(19)

for *t* (0, *L*], where ${R}_{0}^{(0)}(u;{\mathit{\beta}}_{0})={\sum}_{k=1}^{K}{R}_{0k}^{(0)}(u;{\mathit{\beta}}_{0})$.

Consider patient *i*, treated at follow-up time *T _{i}* =

$$\hat{\delta}(t;{T}_{i}|{\mathit{\mathscr{H}}}_{i}({T}_{i}),{T}_{i})={\u015c}_{1}(t;{T}_{i}|{\mathit{\mathscr{H}}}_{i}({T}_{i}),{T}_{i})-{\u015c}_{0}(t;{T}_{i}|{\mathit{\mathscr{H}}}_{i}({T}_{i}),{\mathcal{E}}_{i}({T}_{i})=1,{T}_{i}),$$

(20)

in terms of survival probability, and

$$\hat{\mathrm{\Delta}}(L;{T}_{i}|{\mathit{\mathscr{H}}}_{i}({T}_{i}),{T}_{i})={\hat{\mu}}_{1}(L;{T}_{i}|{\mathit{\mathscr{H}}}_{i}({T}_{i}),{T}_{i})-{\hat{\mu}}_{0}(L;{T}_{i}|{\mathit{\mathscr{H}}}_{i}({T}_{i}),{\mathcal{E}}_{i}({T}_{i})=1,{T}_{i})$$

(21)

in terms of restricted residual mean survival time.

Having established how to estimate the treatment effect for a subject treated at *T _{i}* =

$$E\left[{\int}_{0}^{t}\frac{d{N}_{i}^{T}(u)}{{G}_{i}(u)}\right|{\mathit{\mathscr{H}}}_{i}(u)]={F}_{i}^{T}(t|{\mathit{\mathscr{H}}}_{i}(t)),$$

(22)

where ${F}_{i}^{T}(t|{\mathit{\mathscr{H}}}_{i}(t))=E[{\int}_{0}^{t}dI({T}_{i}\le u)|{\mathit{\mathscr{H}}}_{i}(u)]$ is analogous to the cumulative incidence function for *T _{i}* (with

$${\lambda}_{i}^{C}(t)={\lambda}_{0}^{C}(t)\text{exp}\{{\mathit{\alpha}}_{0}^{\prime}{\mathit{Z}}_{i}(0)\}.$$

(23)

Observed data used to fit model (23) include {*X _{i}, I*(

Finally, estimators of δ(*t*) and Δ(*L*) are given by

$$\hat{\delta}(t)=\frac{{\sum}_{i=1}^{n}{\int}_{0}^{\tau}\hat{\delta}(t;u|{\mathit{\mathscr{H}}}_{i}(u),{T}_{i}=u){\u011c}_{i}{(u)}^{-1}d{N}_{i}^{T}(u)}{{\sum}_{i=1}^{n}{\int}_{0}^{\tau}{\u011c}_{i}{(u)}^{-1}d{N}_{i}^{T}(u)},$$

(24)

$$\hat{\mathrm{\Delta}}(L)={\int}_{0}^{L}\hat{\delta}(t)dt$$

(25)

respectively, where ${\u011c}_{i}(u)=\text{exp}\{-{\hat{\mathrm{\Lambda}}}_{i}^{C}(u)\}$, and with τ satisfying *P*(*X _{i}* ≥ τ) > 0 and typically chosen to be the maximum observed follow-up time.

We assume that the random vectors {*X _{i}, N_{i}*(

*Under certain regularity conditions, n*^{1/2}{(*t*) − δ(*t*)} *and n*^{1/2}{(*L*) − Δ(*L*)} *each converge asymptotically to zero-mean Gaussian processes with covariance functions E*[ξ_{j}(*t*)^{2}] *and*
$E[{\eta}_{j}^{2}]$, *respectively, where* {ξ_{1}(*t*), …, ξ_{n}(*t*)} *and* {η_{1}(*L*), …, η_{n}(*L*)} *are i.i.d. with mean 0 asymptotically. Expressions for* ξ_{i}(*t*) *and*
${\eta}_{i}(L)={\int}_{0}^{L}{\xi}_{i}(t)dt$, *which are quite lengthly, are provided in the*
Supplementary Materials.

Variance estimators for (*t*) and (*L*) are given by ${n}^{-2}{\sum}_{i=1}^{n}{\hat{\xi}}_{i}{(t)}^{2}$ and ${n}^{-2}{\sum}_{i=1}^{n}{\hat{\eta}}_{i}{(L)}^{2}$, respectively; where _{i}(*L*) and _{i}(*t*) are computed by replacing all limiting values by their empirical counterparts. A proof of Theorem 1 is given in the Appendix. The essence of the proof is demonstrating that, asymptotically, ${n}^{1/2}\{\hat{\delta}(t)-\delta (t)\}={n}^{-1/2}{\sum}_{i=1}^{n}{\xi}_{i}(t)+{o}_{p}(1)$ through a sequence of Taylor series expansions and applications of empirical process results.

The proof is provided for the weight, ${\u0174}_{ik}^{A}(t)$. In practice, the stabilized weight, ${\u0174}_{ik}^{B}(t)$ would often be preferred. As implied by Theorem 1, the computation of the variance is quite involved, and such computation becomes more complicated when a stabilizer is incorporated. Such concerns motivate a computationally simpler form for the variance estimator, resulting from taking *Ĝ _{i}*(

We generated treatment-free survival to follow the assumed partly conditional model using methods from Gong and Schaubel (2013). First, subject *i* enters the study on calendar date, *B _{i}*, which is generated from a Uniform(0,

After obtaining the pertinent survival function, transforming the time scale to represent time since cross-section (setting *t _{k}* =

$${\lambda}_{i}({t}_{k}|{Z}_{ia},{Z}_{i}({s}_{ik}),{D}_{i}>{s}_{ik})=\frac{{\lambda}_{0}({t}_{k}+{s}_{ik}){\rho}^{2}{\{{\mathrm{\Lambda}}_{0}({t}_{k}+{s}_{ik})\}}^{({\rho}^{2-1})}}{\text{cos}{(\pi \rho /2)}^{(\rho +1)}}\text{exp}\{{\rho}^{2}{\gamma}_{1}{Z}_{ia}+{\rho}^{2}{\gamma}_{2}{Z}_{i}({s}_{ik})\}.$$

Setting Λ_{0}(*t*) = (*t/a*)^{1/ρ2} and λ_{0}(*t _{k}* +

$${\lambda}_{i}({t}_{k}|{Z}_{ia},{Z}_{i}({s}_{ik}),{D}_{i}>{s}_{ik})=\text{exp}\{{\rho}^{2}{\gamma}_{1}{Z}_{ia}+{\rho}^{2}{\gamma}_{2}{Z}_{i}({s}_{ik})\}/[a\text{cos}{(\pi \rho /2)}^{(\rho +1)}].$$

If we define λ_{i0k}(*t*; *s _{ik}*) = λ

For patients who received treatment prior to dying (*D _{i}* >

The complexity in the data generator is necessary to induce the partly conditional structure of the pre-treatment survival model. The positive stable frailty has become a common choice in the simulation of multivariate survival set-ups due to its preservation of the proportional hazards assumption both conditionally and marginally. Analogous set-ups were used by Zheng and Heagerty (2005) and Gong and Schaubel (2013).

We used *K* = 10 cross section dates, with *CS _{k}* = 100 ×

We present settings where treatment has no effect (δ(*t*) = Δ(*L*) = 0), for which *a*_{1} = [*a* cos(πρ/2)^{(ρ+1)}]^{−1}. We also list results for a setting with a positive treatment effect (δ(*t*) > 0, Δ(*L*) > 0) induced by specifying *a*_{1} = 0.5 × 10^{−4}. In developing appropriate parameter settings, we conceptualized the time scale as representing days. For reporting purposes, time is recorded in years, with results presented for (1), (2), (3) and (3). The weight ${\u0174}_{ik}^{B}(t)$ was used throughout, with the simplified variance estimators applied.

Table 1 presents simulation results for settings with Δ(*L*) = 0 and Δ(*L*) > 0. The quantity Δ(*L*), with *L* = 3, can be interpreted as the difference of 3-year restricted mean survival time due to treatment, among the treated. The proposed estimators appear to be approximately unbiased, with coverage probabilities close to the nominal 95% level. Some degree of under-coverage is observed, which is due to the approximation of the results from Section 3 by treated the (random) weights as fixed. The under-coverage is not in unacceptable amounts, particularly relative to the great reduction in complexity and hence computational burden associated with the approximation.

We examined the performance of the proposed methods under various degrees of model misspecification (see Supplementary Materials). The methods generally perform adequately, although some bias is introduced, and increases with increasing model misfit. The method appears to be most sensitive to misspecification of the treatment initiation hazard.

We applied the proposed methods to estimate the average effect of liver transplantation among the transplanted, by Model for End-stage Liver Disease (MELD) score. This study used data from the Scientific Registry of Transplant Recipients (SRTR). The SRTR data system includes data on all donor, wait-listed candidates, and transplant recipients in the U.S., submitted by the members of the Organ Procurement and Transplantation Network (OPTN), and has been described elsewhere. The Health Resources and Services Administration (HRSA), U.S. Department of Health and Human Services provides oversight to the activities of the OPTN and SRTR contractors.

The study population included patients age ≥ 18 wait listed between 03/01/2002 and 12/31/2009. We excluded patients who were Status 1 (acute liver failure) or previously transplanted. Cross-section dates were chosen every 7 days, 30 days or 90 days from 03/01/2002 to 12/31/2009, which led to *K*=409, 96, or 32 cross sections respectively. The transplant hazard model, ${\lambda}_{ir}^{T}(t)={\mathcal{E}}_{i}(t){\lambda}_{0r}^{T}(t)\text{exp}\{{\mathit{\theta}}_{0}^{\prime}{\mathit{Z}}_{i}(t)\}$, was stratified by United Network for Organ Sharing (UNOS) Region (*r* = 1, …, 11). The covariate, * Z_{i}*(

The pre-transplant death model, ${\lambda}_{i0kr}(t)={\lambda}_{00kr}(t)\text{exp}\{{\mathit{\beta}}_{0}^{\prime}{\mathit{Z}}_{i}({s}_{ik})\}$, was also stratified, where *k* = 1, …, *K* stands for cross section and *r* again denotes UNOS Region. The covariate, * Z_{i}*(

The pre-transplant study sample consisted of *n* = 66, 884 patients, of which 34,539 were observed to receive a deceased-donor liver transplant. For the MELD 30–40 subgroup, weekly cross section dates were chosen. For MELD 18–29 cross sections were drawn monthly. For MELD 6–17, cross sections were drawn every 3 months. Note that, we also tried weekly cross section dates for MELD 6–29 patients, which yielded almost identical results. The analysis was based on the weight, ${W}_{ik}^{B}(t)$.

Figure 3 shows the estimated survival curves for MELD groups 6–8, 15–17, 20–22 and 36–40. Note that the MELD score categories refer to MELD at transplant. Within a MELD category, *Ŝ*_{1}(*t*) can be interpreted as the average survival probability, with *t* representing residual time post-transplant. Analogously, *Ŝ*_{0}(*t*) can be interpreted as the average survival that would have resulted in the absence of liver transplantation, among patients who received a liver transplant. For the MELD 6–8 group, survival in the absence-of-transplantation exceeds post-transplant survival until approximately *t* = 2 years post-transplant. However, *Ŝ*_{1}(*t*) > *Ŝ*_{0}(*t*) for *t* > 2 years, with the distance between the curves widening as *t* increases. The early survival advantage (absence-of-transplant versus with a transplant) for patients in the MELD 6–8 group is the combination of relatively mortality in this subgroup, combined with the risk of surgery-related mortality (not faced unless transplantation occurs). The early survival advantage without transplant is even observed in MELD 15–17 patients, but is much less pronounced and very short-lived. In fact, *Ŝ*_{1}(*t*) > *Ŝ*_{0}(*t*) for *t* > 0.25 years in this subgroup. For MELD 36–40 group, the absence-of-transplant survival curve drops dramatically during the first couple of months, then steadily declines thereafter. Note that *Ŝ*_{1}(*t*) curves are quite similar across MELD subgroups, with *Ŝ*_{0}(*t*) decreasing strongly as MELD increases.

Analysis of SRTR data: Estimated survival curves after with a liver transplant (solid line) and in the absence of liver transplantation (dashed line) among liver transplant recipients. The time axis *t* is years post-transplant.

In Table 2, we list estimates of the difference in survival probability, (*t*) for *t* = 1, 3, 5 years, as well as (5), the difference in 5-year restricted mean residual lifetime. The group that benefits the most from liver transplantation is clearly MELD 36–40, with an average gain in residual survival time of (5) ≈ 2.4 years. The next greatest gain is observed in the MELD 30–35 group, with (5) = 1.4 years. For MELD scores between 15 and 30, there is little difference in the gain in 5-year restricted mean residual survival time, with (5) fluctuating about 1 year across the MELD 26–29, 23–25, 20–22, 18–19 and 15–17 subgroups. Only for the MELD 6–8 group is *H*_{0} : Δ(5) = 0 not rejected.

Analysis of SRTR data: Estimating the effect of liver transplantation on the transplanted (with 95% confidence interval in parentheses), by MELD score at transplant.

In the Supplementary Materials, we provide results based on the Sequential Stratification method (Schaubel, Wolfe and Port, 2006; Schaubel et al., 2009), which features inverse weighted time-dependent stratification to create customized comparisons groups for each subject receiving the time-dependent treatment. Comparing our results in Table 2 to those based on Sequential Stratification, the main difference is in the MELD 6–8 group; the models from Sharma et al (2015) report a hazard ratio of 2.04 (*p* < 10^{−4}), indicating that liver transplant is associated with a doubling of the mortality hazard in this subgroup. In the presence of non-proportionality (which is clear in Figure 3, particularly for this subgroup), the hazard ratio and difference in restricted mean do not have to agree.

Additional analysis is presented in the Supplementary Materials. For each MELD category, multiplying the number of transplants by the (5) yields the number of life-years saved via liver transplantation (considering only the first 5 post-transplant). The largest number of transplants was in the MELD 15–17 category (5,028), but the greatest number of life-years saved (7,649) was in the MELD 36–40 group. We estimate that 34,757 years of life were spared based on the liver transplants observed in this analysis. The Supplementary Materials also present plots of pre-transplant MELD profiles over time, the baseline pre-transplant mortality hazard, the liver transplant baseline hazard, and cumulative incidence of transplantation.

In this report, we develop methods for estimating the average effect on the treated of a time-dependent treatment. The methods can be used to evaluate the benefit, in terms of patient survival, of a treatment under current treatment assignment practices. The methods were applied to quantify the survival benefit of deceased-donor liver transplantation among the transplanted, by Model for End-stage Liver Disease (MELD) score.

The proposed methods are not intended to guide treatment decisions. For example, the fact that we estimate a larger treatment effect for MELD 36–40 than for 30–35 does not imply that a patient with MELD=32 should wait until his/her MELD score increases to ≥ 36 before they agree to be transplanted. The proposed methods cannot generally be used to compare treatment effects, since each treatment effect is averaged differently. For example, the difference in the treatment effect between patients transplanted at MELD 15–17 ((5) = 1.00) and MELD 12–14 ((5) = 0.59) is partly attributable to the former group being transplanted with higher quality donor livers.

There are now many methods available for evaluating a time-dependent treatments. Marginal Structural Models (MSM; e.g., Hernán, Brumback and Robins, 2000; Robins, Hernán and Brumback, 2000) are not well-suited to our set-up due to the potential for treatment to interact with time-varying covariates. Structural Nested Failure Time Models (SNFTMs; e.g., Robins, 1988; Joffe et al., 1998; Keiding et al., 1999; Hernan et al., 2005; Taubman et al., 2009; Vock et al., 2013) are an alternative. In particular, the method of Vock et al. (2013) was motivated by the lung transplant setting. Versions of Sequential Stratification, which involves stratified and inverse weighted Cox regression, have been used to evaluate the benefit of kidney transplantation (Schaubel, Wolfe and Port, 2006) and liver transplantation (Schaubel et al., 2009). An advantage the proposed method over SNFTMs and Sequential Stratification is the avoidance of any parametric assumptions regarding the treatment effect. SNFTMs assume that treatment alters the time scale through a constant, while Sequential Stratification assumes proportionality of the pre- and post-treatment hazard functions. A further advantage of our proposed methods over SNFTMs relates to implementation. Although explicit coding would be required for either approach, the ‘core’ models in our method merely involve Cox regression and, therefore, can be fitted using standard statistical software (SAS, R) after modifying the input data appropriately.

In estimating the ETT, we consider the absence of treatment; i.e., *T _{i}* = ∞. In setting where this is found to be too ambitious a goal (e.g., lack of sufficiently long follow-up, in a setting where treatment is inevitable), one could change [

An alternative to the measures proposed in (9) and (10) would be to redefine *S*_{1}(*t*) to be the population average survival (i.e., averaging over the current treated and untreated experiences), with *S*_{0}(*t*) then representing the average population survival in the absence of treatment. Unless strong or unrealistic assumptions were made, the ‘core’ models for this approach would be quite similar to those in the proposed approach, except for the pre-treatment hazard model. The proposed averaging would be preferred in many practical settings (including the liver transplant setting which motivated our current work) since the absence of a treatment benefit among non-recipients is made explicit.

This work was supported in part by National Institutes of Health Grant R01-DK070869. The data reported here have been supplied by the Minneapolis Medical Research Foundation (MMRF) as the contractor for the Scientific Registry of Transplant Recipients (SRTR). The interpretation and reporting of these data are the responsibility of the authors and in no way should be seen as an official policy of or interpretation by the SRTR or the U.S. Government. The authors wish to thank the Associate Editor and two Reviewers, whose comments and suggestions led to considerable improvement of the manuscript. They also thank Min Zhang for her many thoughtful suggestions.

Supplementary Materials

Supplementary Materials, referenced in Sections 3, 4 and 5, are available with this paper at the Biometrics website on Wiley Online Library.

- Andersen PK, Gill RD. Cox’s regression model for counting processes: A large sample study. The Annals of Statistics. 1982;10:1100–1120.
- Breslow NE. Contribution to the discussion of paper by D.R. Cox. Journal of the Royal Statistical Society Series B. 1972;34:216–217.
- Chen P, Tsiatis AA. Causal inference on the difference of the restricted mean life between two groups. Biometrics. 2001;57:1030–1038. [PubMed]
- Cox DR. Regression models and life tables (with Discussion) Journal of the Royal Statistical Society, Series B. 1972;34:187–200.
- Cox DR. Partial likelihood. Biometrika. 1975;62:269–275.
- Feng S, Goodrich NP, Bragg-Gresham JL, Dykstra DM, Punch JD, DebRoy MA, Greenstein SM, Merion RM. Characteristics associated with liver graft failure: The concept of a donor risk index. American Journal of Transplantation. 2006;6:783–790. [PubMed]
- Feuer EJ, Hankey BF, Gaynor JJ, Wesley MN, Baker SG, Meyer JS. Graphical representation of survival curves associated with a binary non-reversible time dependent covariate. Statistics in Medicine. 1992;11:455–474. [PubMed]
- Gong Q, Schaubel DE. Partly conditional estimation of the effect of a time-dependent factor in the presence of dependent censoring. Biometrics. 2013;69:338–347. [PMC free article] [PubMed]
- Hernán MA, Cole SR, Margolick J, Cohen M, Robins JM. Structural accelerated failure time models for survival analysis in studies with time-varying treatments. Pharmacoepidemiology and drug safety. 2005;14:477–491. [PubMed]
- Joffe MM, Hoover DR, Jacobson LP, Kingsley L, Chmiel JS, Visscher BR, Robins JM. Estimating the effect of Zidovudine on Karposi’s sarcoma form observational data using a rank preserving structural failure time model. Statistics in Medicine. 1998;17:1073–1102. [PubMed]
- Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association. 1958;282:457–481.
- Kalbfleisch JD, Prentice RL. The Statistical Analysis of Failure Time Data: 2nd Edition. New York: Wiley; 2002.
- Keiding N, Filiberti M, Esbjerg S, Robins JM, Jacobsen N. The graft versus leukemia effect after bone marrow transplantation: A case study using structural nested failure time models. Biometrics. 1999;57:23–28. [PubMed]
- Parast L, Tian L, Cai T. Landmark estimation of survival and treatment effect in a randomized clinical trial. Journal of the American Statistical Association. 2014;109:383–394. [PMC free article] [PubMed]
- Pearl J. Causality: Models, Reasoning, and Inference. New York: Cambridge; 2009.
- Robins JM. The control of confounding by intermediate variables. Statistics in Medicine. 1988;8:679–701. [PubMed]
- Robins JM, Rotnitzky A. Recovery of information and adjustment for dependent censoring using surrogate markers. In: Jewell N, Dietz K, Farewell B, editors. AIDS Epidemiology - Methodological Issues. Boston: Birkhäuser; 1992. pp. 297–331.
- Rubin DB. Bayesian inference for causal effect: The role of randomization. Annals of Statistics. 1978;6:34–58.
- Schaubel DE, Wolfe RA, Port FK. A sequential stratification method for estimating the effect of a time-dependent experimental treatment in observational studies. Biometrics. 2006;62:910–917. [PubMed]
- Schaubel DE, Wolfe RA, Sima CS, Merion RM. Estimating the effect of a time-dependent treatment by levels of an internal time-dependent covariate. Journal of the American Statistical Association. 2009;104:49–59.
- Sharma P, Schaubel DE, Goodrich NP, Merion RM. Serum sodium and the survival benefit of liver transplantation. Liver Transplantation. 2015;21:308–313. [PMC free article] [PubMed]
- Van Houwelingen HC. Dynamic prediction by landmarking in event history analysis. Scandinavian Journal of Statistics. 2007;34:70–85.
- Van Houwelingen HC, Putter H. Dynamic prediction in clinical survival analysis. New York: CRC Press; 2012.
- Van Houwelingen HC, Putter H. Comparison of stopped Cox regression with direct methods such as pseudo-values and binomial regression. Lifetime Data Analysis. 2015;21:180–196. [PubMed]
- Vock DM, Tsiatis AA, Davidian M, Laber EB, Tsuang WM, Finlen Copeland CA, Palmer SM. Assessing the causal effect of organ transplantation on the distribution of residual lifetime. Biometrics. 2013;69:820–829. [PMC free article] [PubMed]
- Zhang M, Schaubel DE. Estimating differences in restricted mean lifetime using observational data subject to dependent censoring. Biometrics. 2011;67:740–749. [PMC free article] [PubMed]
- Zhang M, Schaubel DE. Double-robust semiparametric estimator for differences in restricted mean lifetimes in observational studies. Biometrics. 2012;68:999–1009. [PMC free article] [PubMed]
- Zheng YY, Heagerty PJ. Partly conditional survival models for longitudinal data. Biometrics. 2005;61:379–391. [PubMed]
- Zucker DM. Restricted mean life with covariates: Modification and extension of a useful survival analysis method. Journal of the American Statistical Association. 1998;93:702–709.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |