Assessing the effect of a partly unobserved, exogenous, binary time-dependent covariate on survival probabilities using generalised pseudo-values

The compartment representation for transplant data

By focussing on cumulative treatment effects and prediction of survival probabilities as primary aim, the most interesting outcome of such studies is survival at a pre-specified long-term time point \( {t}^{\ast } \), e.g., in order to investigate 5-years survival probabilities, \( {t}^{\ast } \) would be set equal to 5 years. In line with the approach of genetic randomisation, two separate populations have to be distinguished conditional on their donor availability status, and survival probabilities in patients with and without a donor available, \( {S}_1\left({t}^{\ast}\right) \) and \( {S}_0\left({t}^{\ast}\right) \), have to be compared, respectively. Donor availability is defined as the identification of a donor in the time interval from study entry to maximum donor search time \( {t}^{search} \), with \( {t}^{search}\le {t}^{\ast } \).

Survival in patients with donor available

At first, survival \( {S}_1\left({t}^{\ast}\right) \) in patients with an available donor is investigated. Let \( W \) denote a random variable representing the (waiting) time from a given origin (time 0) to the time when a donor is identified, \( 0\le W\le {t}^{search} \). The distribution of \( W \) is characterised by the density \( {f}_{01}(t) \) with corresponding hazard function \( {\lambda}_{01}(t) dt=P\left(W\in \left[t,t+ dt\right)|W\ge t\right) \). In theory, \( W \) can always be observed for this population independently of the survival or censoring status of the patients by prolonging donor search until \( {t}^{search} \). In practice, as outlined above, donor availability remains unobserved in patients that die or become censored before donor identification.

Survival in the population with available donor corresponds to a stochastic process with three states (Fig. 1a), where all patients start at the initial state 0.

Let \( T \) denote the failure time to reach state 2. The absorbing state 2 can be reached either directly (0→2) if \( T<W \) or through the intermediate state 1 (0→1→2) if \( T\ge W \). Let \( {\lambda}_{02}(t) \) denote the hazard function for a transition at time \( t \) from the initial state 0 direct to the absorbing state 2. Accordingly

Note that in patients with a donor (Fig. 1a), no patients remain in state 0 at \( {t}^{search} \) as all patients have either moved to state 1 or to state 2 until \( {t}^{search} \).

Survival in patients without donor available

Since donor availability is stochastically independent of patients’ prognosis under standard treatment, \( {\lambda}_{02}(v)={\lambda}_{02}^{\hbox{'}}(v) \) until \( {t}^{search} \). This independence has two important consequences. Firstly, \( {S}_0\left({t}^{\ast}\right) \) is identical to the counterfactual survival probabilities in a population where there is no change in therapeutic strategies after a donor is identified, i.e. \( {\lambda}_{12}(v)={\lambda}_{02}^{\hbox{'}}(v),v\ge w \); secondly, in order to estimate \( {S}_0\left({t}^{\ast}\right) \) the survival information of the donor available group can be exploited up to time \( w \).

Pseudo-values

In the following the original pseudo-value approach is briefly outlined. Subsequently, generalised pseudo-values for the estimation of \( {S}_0\left({t}^{\ast}\right) \) and \( {S}_1\left({t}^{\ast }|w\right) \) are introduced. Furthermore, appropriate weights for \( {S}_1\left({t}^{\ast }|w\right) \) are defined to compensate for partly unobserved waiting times when estimating \( {S}_1\left({t}^{\ast}\right) \).

Common pseudo-values

The pseudo-value regression technique [25–30] for censored survival data provides an attractive alternative to methods commonly applied in survival analysis. The approach specifically allows the modelling of survival probabilities at pre-specified interesting time points. Without relying on proportional hazards, pseudo-values can be flexibly regressed on common baseline covariates within the framework of a generalised linear model [15].

Let \( {T}_i \) denote the time to failure and \( {C}_i \) the time to censoring for the \( i \)-th patient, \( i=1,\dots, n \). Usually \( {\tilde{T}}_i=\min \left({T}_i,{C}_i\right) \) and \( {D}_i=I\left({T}_i\le {C}_i\right) \) can be observed. Let \( \widehat{S}(t) \) denote the ordinary Kaplan-Meier estimate and let \( {\widehat{S}}^{-i}(t) \) denote the Kaplan-Meier estimate where the \( i \)-th observation is excluded (jacknife statistic). The pseudo-value at time \( {t}^{\ast } \) [25] for the \( i \)-th patient is now defined as

Individual pseudo-values can attain values below zero and above one. The approach relies on the fact that the Kaplan-Meier estimate is (approximately) unbiased for the marginal survival function [25]. Accordingly, i.i.d. observations, independent right censoring and a sufficiently large risk set at time \( {t}^{\ast } \) have to be assumed. Asymptotic conditional unbiasedness of the pseudo-values given covariates can be shown [25, 32–34], which allows the use of regression models with pseudo-values as outcome (further details can be found in [25, 29]).

Generalised 0➔2 pseudo-values for the estimation of \( {S}_0\left({t}^{\ast}\right) \)

For notational convenience let the first \( m \) out of \( n \) patients be those with an observed transition to state 1 at times \( {w}_1,\dots, {w}_m \).

Let \( {\widehat{S}}_0(t) \) denote the Kaplan-Meier estimate for the risk of a direct transition from state 0 to state 2 with censoring at \( {w}_1,\dots, {w}_m \) of those patients with observed transitions to state 1. Hence, only direct transitions from state 0 to state 2 are considered as event. Given the assumed independence between the two stochastic processes (transitions 0→1 and 0→2), \( {\widehat{S}}_0(t) \) is an asymptotically unbiased estimate of \( {S}_0(t) \).

As all patients start in state 0 and are at risk for a direct transition to state 2, it is of importance that a 0→2 pseudo-value has to be calculated for each patient. That is, patients with a transition to state 1 are censored at \( {w}_1,\dots, {w}_m \) but not ignored, as the latter would lead to selection bias due to over-representing direct transitions from state 0 to state 2. Since \( {\widehat{V}}_{i,0}\left({t}^{\ast}\right) \) is an asymptotically unbiasedly estimate for \( {S}_0\left({t}^{\ast}\right) \), its mean \( {\overline{V}}_0=\frac{1}{n}{\sum}_{i=1}^n{\widehat{V}}_{i,0}\left({t}^{\ast}\right) \) is asymptotically unbiased as well [25, 32].

Generalised 0➔1➔2 pseudo-values for the estimation of \( {S}_1\left(\left.{t}^{\ast}\right|w\right) \)

Utilising the independence between \( {\widehat{S}}_0\left({w}_i\right) \) and \( {\widehat{U}}_{i,1}\left({t}^{\ast}\left|{w}_i\right.\right) \), \( {\widehat{V}}_{i,1}\left({t}^{\ast}\left|{w}_i\right.\right) \) is asymptotically unbiased too,

\( E\left[{\widehat{V}}_{i,1}\left({t}^{\ast}\left|{w}_i\right.\right)\right]={S}_0\left({w}_i\right)S\left({t}^{\ast}\left|T\ge {w}_i,W={w}_i\right.\right)+{o}_p(1)={S}_1\left({t}^{\ast}\left|{w}_i\right.\right)+{o}_p(1). \)

Estimation of \( {S}_1\left({t}^{\ast}\right) \)

\( {S}_1\left({t}^{\ast}\right)=\underset{0}{\overset{t^{\ast }}{\int }}{f}_{01}(w)\;{S}_1\left(t|w\right) dw. \)

Note that, unless \( {\lambda}_{02}(v)=0 \) and \( {\lambda}_C(v)=0 \) for \( v<{t}^{search},q(w)\ne {f}_{01}(w) \) and consequently \( {\overline{S}}_1\left({t}^{\ast}\right)\ne {S}_1\left({t}^{\ast}\right) \). Since for some patients donor availability is unknown due to censoring or 0→2 transitions before time \( w \), long waiting times are under-represented in \( q(w) \) compared to \( {f}_{01}(w) \).

The weights are then \( {\widehat{\gamma}}_{i,1}=\frac{1}{\widehat{G}\left({w}_i\right)}{\widehat{p}}_m \), \( i=1,\dots, m \), and \( {\widehat{p}}_m=\frac{m}{\sum \limits_{j=1}^m\frac{1}{\widehat{G}\left({w}_j\right)}} \) so that \( \sum \limits_{i=1}^m{\widehat{\gamma}}_{i,1}=m \).

Now, the weighted mean \( {\overline{V}}_1=\frac{1}{m}\sum \limits_{i=1}^m{\widehat{\gamma}}_{i,1}{\widehat{V}}_{i,1}\left({t}^{\ast}\left|{w}_i\right.\right) \) is a consistent estimator of \( {S}_1\left({t}^{\ast}\right) \) [35].

Weighted generalised linear model for comparing \( {S}_0\left({t}^{\ast}\right) \) with \( {S}_1\left({t}^{\ast}\right) \)

Analogous to the original pseudo-value approach, a weighted generalised linear model is utilized for inferential purposes with log-log link function, normal response probability distribution, and the parameters \( {\beta}_0 \) and \( {\beta}_1 \). Whereas \( {\beta}_0 \) corresponds to the intercept, \( {\beta}_1 \) corresponds to the binary indicator \( {x}_{1i} \), with \( {x}_{1i}=0 \) for \( i=1,\dots, n \) and \( {x}_{1i}=1 \) for \( i=n+1,\dots, n+m \). The \( \left(n+m\right) \)-dimensional response vector is \( \widehat{\mathbf{V}}={\left({\widehat{V}}_{1,0}\left({t}^{\ast}\right),\dots, {\widehat{V}}_{n,0}\left({t}^{\ast}\right),{\widehat{V}}_{1,1}\left({t}^{\ast}\left|{w}_1\right.\right),\dots, {\widehat{V}}_{m,1}\left({t}^{\ast}\left|{w}_m\right.\right)\right)}^{\hbox{'}} \). The weights \( {\widehat{\gamma}}_{i,1} \) for \( {\widehat{V}}_{i,1}\left({t}^{\ast}\left|{w}_i\right.\right) \) are defined in the previous chapter and the weights \( {\widehat{\gamma}}_{i,0} \) for \( {\widehat{V}}_{i,0}\left({t}^{\ast}\right) \) are set to one.

The parameter estimates are functions of the weighted means of the two types of pseudo-values; \( {\widehat{\beta}}_0=g\left({\overline{V}}_0\right) \) and \( {\widehat{\beta}}_1=g\left({\overline{V}}_1\right)-g\left({\overline{V}}_0\right) \). When the log-log link function is used, \( g(V)=\log \left(-\log (V)\right) \) and \( \exp \left({\widehat{\beta}}_1\right) \) is an asymptotically unbiased estimate of the cumulative hazard ratio at time \( {t}^{\ast } \) comparing patients with and without donor available.

Estimation of standard errors

When using the Kaplan-Meier estimate \( {\widehat{S}}_0\left({w}_i\right) \) in the estimation of \( {\widehat{V}}_{i,1}\left({t}^{\ast}\left|{w}_i\right.\right),i=1,\dots, m \) (eq. 6), the individual patients’ variation is not properly represented. Consequently, the weighted generalised linear model with a sandwich estimator underestimates the standard errors of \( \widehat{\beta^{\hbox{'}}}=\left({\widehat{\beta}}_0,{\widehat{\beta}}_1\right) \) on average. More appropriate standard errors can be attained by replacing \( {\widehat{S}}_0\left({w}_i\right) \) in eq. 6 by random draws from Bernoulli variables \( {\mathrm{B}}_i \) with

\( {\mathrm{B}}_i\sim \mathrm{Bernoulli}\left[\exp \left\{-\exp \left({p}_i\right)\right\}\right] \) where \( {p}_i\sim N\left(\log \left[-\log \left\{{\widehat{S}}_0\left({w}_i\right)\right\}\right],\frac{{\widehat{\sigma}}^2\left({\widehat{S}}_0\left({w}_i\right)\right)}{{\left[{\widehat{S}}_0\left({w}_i\right)\kern0.24em \log \left\{{\widehat{S}}_0\left({w}_i\right)\right\}\right]}^2}\right) \)

and \( {\widehat{\sigma}}^2\left({\widehat{S}}_0\left({w}_i\right)\right) \) denotes the variance estimator of the Kaplan-Meier estimator (Greenwood’s formula). This ad-hoc approach exploits the asymptotic normality of the log-log transformation of the Kaplan-Meier estimator [31]. In practical applications, the imputations should be repeated several times and the obtained standard errors should be averaged for stability reasons.

Software implementation

The proposed method can be straightforwardly implemented using standard routines available in the majority of statistical software packages. Details on the implementation in SAS and R are described in the Additional file 1 (Section D).

Simulation studies

Simulation study 1

This simulation study was performed to examine the properties of individual components of our approach, i.e. (1) the weights \( {\widehat{\gamma}}_{i,1} \) to estimate \( {S}_1\left({t}^{\ast}\right) \), (2) the survival estimates for \( {S}_0\left({t}^{\ast}\right) \), \( {S}_1\left({t}^{\ast}\left|{w}_i\right.\right) \) and \( {S}_1\left({t}^{\ast}\right) \), and (3) the standard error estimates. Details of the simulation setup are provided in scenario I in Table S1 of the Additional file 1.

Seventy-five percent of all patients have a donor available. These patients equally split between discrete waiting times at \( w \) = 0.5, 1 and 3 years and the according probability mass function is \( {f}_{01}(w)=1/3 \). The true survival probabilities and simulation results are given in Table 1.

For a sample size of 1000, the means of the observed distribution of waiting times \( \widehat{q}(w) \) are 0.46, 0.39 and 0.15 for \( w \) = 0.5, 1 and 3, respectively. Hence short waiting times are over- and long-waiting times are underrepresented compared to \( {f}_{01}(w) \). The means of the estimated weights \( {\widehat{\gamma}}_{i,1} \) are 0.72, 0.86 and 2.26 for \( w \) = 0.5, 1 and 3, respectively, and the correct waiting time distribution can be restored for all \( w \) with \( {\widehat{f}}_{01}(w)=0.33 \).

Table 1 also considers estimates from a weighted generalised linear model with a log-log link with and without the ad-hoc correction for standard error estimation. For a sample size of 1000 the bias of the \( \log \left(-\log \left(\mathrm{S}\right)\right) \) estimate for both approaches is negligible with a maximum absolute value of 0.011. SE_est is the mean of the empirical sandwich estimates of the standard error of the \( \log \left(-\log \left(\mathrm{S}\right)\right) \) estimates from the 1000 simulation runs (Table 1). Note that in the analysis of a single data set, several repeated imputations are needed to get stable standard error estimates. The mean empirical standard errors are already replicated within the 1000 simulation runs and so no repeated imputations are needed to get stable results.

Monte Carlo standard deviations are similar to mean SE_est and shown in Table S2 (Additional file 1). In the uncorrected case, SE_est tends to become smaller for increasing \( w \). As expected, SE_est is generally larger in the ad-hoc corrected case. Consequently, with uncorrected SE_est estimation the coverage of the 95% confidence intervals decreases with increasing \( w \) and is only 86.4% for \( w=3 \). When \( {\widehat{V}}_{i,1}\left({t}^{\ast}\right) \) are dervied using the ad-hoc approach, the coverage substantially improves to 93.7% for \( {S}_1\left(5\left|w\right.=3\right) \). The results for a sample size of 400 show a similar good performance with regard to unbiasedness of survival estimates and confidence interval coverage (Table 1).

Stimulation study 2

This simulation study was performed to demonstrate both, approximate unbiasedness of the parameter estimates and approximate confidence interval coverage within the weighted generalised linear model with a log-log link. The standard error ad-hoc correction was used throughout in this section.

Four realistic scenarios from pediatric oncology and four theoretical scenarios are considered. Scenarios A-D are based on published results [5, 37–39] and represent realistic situations that cover all main potential types of departure from proportional hazards: differences in long-term survival only (A), crossing survival curves (B-C) and a situation where only differences in short-term survival (D) are observed. Additionally, a proportional hazards model (E), a null-model (F) and a scenario with unrealistically long-waiting times (G) were investigated. Scenario G is implemented under both, the commonly used and a much more pronounced censoring scheme, that is 0–11 and 0–6 years uniform censoring, respectively. True survival curves of patients with and without donor available are shown in Additional file 1: Figure S1 for the various scenarios together with the corresponding waiting time distributions. Details of the simulation setup are provided in scenarios A-G Table S1 and in Section C in the Additional file 1.

The results are summarised in Table 2. In scenario A, initially similar survival is observed until year 1 and afterwards stem-cell transplantation has lower hazards, leading to an inferior long-term survival probability with chemotherapy alone; \( {S}_0(5)=40.4\% \) compared to \( {S}_1(5)=56.2\% \). Accordingly, using the log-log link, \( {\beta}_0,{\beta}_1 \) and the cumulative hazard ratio (cHR) are −0.098, −0.453 and 0.636, respectively. For a sample size of 1000, the bias is 0.001 and −0.001 for \( {\widehat{\beta}}_0 \) and \( {\widehat{\beta}}_1 \), respectively. Similar small biases of 0.001 and −0.002 are seen for a sample size of 400, respectively. The mean standard errors of \( {\widehat{\beta}}_0 \) (SE_est) obtained from the weighted generalised linear models are 0.053 and 0.085 for sample sizes of 1000 and 400, respectively. Furthermore, the observed coverages of the 95% confidence intervals (CI) for \( {\beta}_0 \) are 96.2% and 95.4%, respectively. Similar satisfying results are seen for the CI coverages of \( {\beta}_1 \) and \( {\beta}_0+{\beta}_1 \).

Likewise, convincing results are seen for Scenarios B-G as well (Table 2). In particular scenario G with uniform censoring between 0 and 6 is challenging - long waiting times coincides with heavy censoring. Even in this difficult situation, the biases of \( {\widehat{\beta}}_0 \) and \( {\widehat{\beta}}_1 \) are close to zero; the corresponding CI coverages are 95.0% and 96.0% for \( n \) = 1000 and 94.9% and 96.5% for \( n \) = 400, respectively.

Monte Carlo standard deviations are similar to mean SE_est and shown in Table S3 (Additional file 1).

In summary, in all scenarios the biases are negligibly small. Note that, the maximum bias on the scale of survival probabilities is smaller than one percentage point over all scenarios and sample sizes. The coverage probabilities of the confidence intervals are convincingly close to the nominal 95% level. As expected, the good performance of the method does not depend on the proportional hazards assumption.

Benefit of SCT in paediatric leukaemia patients

A real dataset from a recently published international study in children with newly diagnosed Philadelphia chromosome-positive (PH+) acute lymphoblastic leukaemia [40] is used for illustrative purposes. The aim of the study is to compare SCT to conventional chemotherapy only with disease free survival (DFS) as primary endpoint. A total of 542 patients were eligible and included in the study. Of these, 217 were treated with chemotherapy only and 325 switched to SCT after a suitable donor was identified. Note that, here we only have information when an SCT was actually performed and we assume that patients received SCT immediately after their donor was identified. SCT was performed after a median waiting time of 5.1 months; the majority (95%) of SCTs were performed within 1 year.

The application of the generalised pseudo-values approach is contrasted to the popular Cox model. The original analysis allowed for non-proportional hazards [40] by including a time-dependent treatment indicator and its interaction with log(time) in a Cox regression model. The hazard ratios at 6 months and 5 years were estimated to be 1.56 and 0.39, respectively (Fig. 2). While there is a clear evidence of lower hazards with SCT at later time points, the average hazard ratio of 0.91 (p = 0.39) from a Cox regression model with a time-dependent treatment indicator shows little evidence of a beneficial SCT effect (Fig. 2). However due to the ignorance of the non-proportional hazards these results may be overly sensitive to the initial disadvantage of SCT and may understate the positive impact of SCT in the long term.

For 5-years DFS, the generalisation of pseudo-value regression for a time-dependent covariate is applied with a log-log link and \( {t}^{search}={t}^{\ast }=5 \) years. The log-cumulative hazard ratio is −0.26 and the estimated 5-year cHR is 0.77. Using the ad-hoc approach as described in subsection ‘Estimation of standard errors’ with 1000 repeated imputations, the mean standard error of the log-cumulative hazard ratio is 0.123. The corresponding 95% confidence interval for cHR is 0.61 to 0.98 (Fig. 2) showing a positive cumulative effect of SCT at 5 years with a Wald-type p-value of 0.042. This cHR is directly related to the 5-year DFS probability estimates of 42% and 32% with and without donor, respectively.