Analysis of survival data with non-proportional hazards: A comparison of propensity score weighted methods

In the setting of survival outcomes with non-proportionality, we must first decide what our estimand will be. It should be readily interpretable by the intended audience, meaningful given the observed data features, and estimable using available software. If we fit a proportional hazards model to non-proportional data, we would obtain an estimate of a single HR, but this would be a misleading measure because it would essentially be a summary of the different hazard ratios at each observed failure time. Alternative hazard ratios have been proposed including the average hazard ratio which is interpreted as the average of the hazard ratios weighted by the number of patients at risk. [24]

The difficulties inherent in reporting hazard ratios in the setting of non-proportionality can be avoided by focusing on measures that can be defined directly from the survival functions for the treated and control arms, $S_{1}(t)$ and $S_{0}(t)$, respectively.

The simplest approach compares differences in survival probabilities at a particular time $t=T$.

Quantiles of the survival functions can also be compared. For instance, if the survival curves are defined for $S(t)=0.5$, median survival times can be compared.

These measures have the benefit of being readily understood by a clinical audience, but the drawback is that they do not comprehensively measure survival effects over the length of follow-up. An alternative measure is restricted mean survival (RMS), the mean survival up through some time $T$ where $S(t)$ is defined (i.e. mean survival up through the maximum observation time). [25] RMS can be found using the area under the survival curve, where the difference in RMS between the treated and control conditions is

Each of these measures are estimated directly from the survival curve, so they are always meaningful regardless of proportionality. However, not every method for the analysis of survival data explicitly estimates the survival curve, and not all methods can accommodate weighting. In the subsequent sections, we focus on methods which can 1) accommodate IPTW, 2) allow the estimation of the measure of interest ($\Delta S_{T}$, $\Delta T_{50}$, and $\Delta\mathrm{RMS}_{T}$), and 3) are readily implemented using published software.

Although the Cox proportional hazards model does rely on the assumption of proportionality, the model tends to be robust to small deviations from this assumption. However, if $n$ is large, as is common is registry studies, small deviations from proportionality may be statistically significant when tested using standard diagnostic tests (e.g. Cox-Snell residuals, Schoenfeld residuals). [26] Therefore, the Cox model may still be useful in some cases were deviations from proportionality are demonstrated. The Cox model is defined as follows:

IPTW can be readily accommodated in the Cox model, by weighting the contribution of each observation to the partial likelihood function:

where $R(t_{i})$ is the risk set at time $t_{i}$. The parameter $\beta$ is the average of the log-hazard ratios at each failure time. Importantly, one can find $\hat{S}(t)$ based on the estimates $\hat{h}_{0}(t)$ and $\hat{\beta}$.

We can also extend the standard Cox model, relaxing the assumption of proportionality, by allowing the hazard ratio to vary as a function of time.

In our simulations, we consider two particular cases, 1) allowing the hazard to vary by log-time.

and 2) using a piecewise constant treatment effect

These models can all be implemented with the widely-used survival package in R. [27]

Accelerated Failure Time (AFT) models assume a linear model for $\mathrm{log}(T)$. [3, 7, 28]

Here, we consider parametric models, where the error term $\epsilon$ is assumed to take on a given distribution. We have many choices for the distribution of $\epsilon$. A common choice is the Gumbel distribution, which yields a Weibull hazard model; however, the Weibull model is also a proportional hazards model. We can allow for non-proportionality by allowing the hazard ratio to change over time

where $\lambda$ is the scale and $\gamma$ is the shape parameter. If we let

then (1) can be re-written as

Therefore, if we include an additional parameter in the Weibull model, where the shape can vary by treatment, the model specification allows the hazard ratio to vary by log-time.

Another flexible alternative is the generalized gamma model, which is also a three parameter model. [29] The distribution of the errors has an additional parameter $Q$.

This model encompasses other distributions as special cases; when $Q=1$ we have a Weibull AFT model, and if $Q=0$ it results in a log-normal model.

For any AFT model, the weights $\hat{w}_{i}$ can be applied during estimation of parameters. Here, we focus on the 3-parameter Weibull and generalized gamma AFT models, but a wide range of AFT models can be fit in R using the flexsurv package. [30]

Another strategy is to avoid parametric assumptions entirely, and use Kaplan-Meier curves with IPTW to directly estimate the survival functions. [31, 32]

Where $\widehat{d_{j}^{w}}=\sum_{i:T_{i}=t_{j}}\hat{w}_{i}\delta_{i}$ is the weighted number of events and $\widehat{n_{j}^{w}}\sum_{i:T_{i}\geq t_{j}}\hat{w}_{i}$ is the weighted number of people at risk at time $t_{j}$. Weighted curves are easily estimated in R via the survival package.

Another non-parametric method which can be used to estimate the survival curve is the recently-developed pseudo-observations approach. [4, 5, 6] This framework uses a missing data approach to account for censoring. If no missing data or confounders are present, the ATE can be found using a simple average over the (possibly transformed) outcomes for each subject

But if $f(Y_{i})$ is unknown for some subjects, we cannot calculate $E(f(y))$ directly. Instead, we can define the pseudo-observation for subject $i$ as

where $\hat{\theta}$ is a consistent estimate for $\theta$, and $\hat{\theta}^{-i}$ is the estimator applied to all observations excluding subject $i$. The pseudo-observations are then used in place of all observations, not only the ones which are missing due to censoring.

In a survival context, our parameter of interest, $\theta$, is the probability of survival at time $t$, and $\theta_{i}$ is a binary indicator for the status (alive/dead) of subject $i$ at time $t$. We can use a non-parametric method, such as Kaplan-Meier, to estimate $S_{Z}(t)$, which is then used to find $\theta_{i}$ for all observations in arm $Z$.

Heuristically, $n_{Z}\hat{S}(t)$ is the expected number of patients alive at time $t$, and $(n_{Z}-1)\hat{S}(t)^{-i}$ is the expected number of patients other than patient $i$ alive at time $t$. $\hat{\theta_{i}}$ is used in place of (possibly unobserved) survival status at time $t$. Therefore, the probability of being alive at time $t$ in arm $Z$ is

Incorporating propensity-score weights to account for confounding is straightforward.

There are several benefits of using the pseudo-observations approach. It has been demonstrated to provide a valid estimate of the ATE in a causal framework, and requires weaker assumptions than regression-based models. It also supports estimation of many outcomes of interest. We can estimate $\hat{\mathrm{P}}(Y|Z,t)$ at each failure time, defining the survival curve, which in turn allows us to estimate the measures of interest described in section 3.1. Pseudo-observations can also be used to estimate RMS. This method can be implemented using the R package pseudo. [33]

To conduct inference for the outcomes of interest, we need to quantify uncertainty. Some methods have closed-form variance estimators for our estimands of interest, but others do not. This motivates the use of a non-parametric bootstrap to estimate variances and find confidence intervals.[34] Bootstrap-based estimates are available for every combination of model and outcome we discuss above. Further, to correctly account for the variance of $\hat{e}$, the estimated propensity scores, we can re-estimate the propensity score within each bootstrap iteration. [35]

We tested the performance of the methods described above using simulation studies. We assessed the bias, coverage, and variance of each method’s estimates of our outcomes of interest: median survival, restricted mean survival, and survival probability at time $T$. We based our covariate and outcome distributions, and effect sizes on our NCDB cancer studies to mimic real world scenarios.

First, we simulated a set of covariates based on the observed multivariate distribution of the NCDB covariates in the renal dataset. We used the R package genOrd [36] to draw variables representing gender, age, stage, histology, tumor grade, Charlson comorbidity score, race, ethnicity, insurance status, facility type, income, and education. As our goal here is to evaluate the survival analysis methods, and not the propensity score estimation methods, we held this set of covariates constant for each simulation; we employ simulated covariates instead of the actual covariates so that our data and simulation code can be shared.
(Available online at https://github.com/BethHandorf/NonPH_IPTW)

We used a simple logistic model with linear covariate effects to determine each subjects’ probability of being assigned to the treatment condition (versus control). The effect size (log-scale) for each covariate ranged from 0.8 (effect of high vs low stage) to 0.03 (for each additional year of age). (See Appendix Table 1 for full specification.) This corresponds roughly to the effect sizes estimated when modeling actual treatment allocation in the NCDB sample. This model defined the probability that each individual would be allocated to the treatment condition. We set the intercept such that, on average, the population’s probability of receiving the treatment was 0.5.

Survival times were drawn based on a Weibull model with a time-varying treatment effect.

In our base case, the hazard ratio varied by log time $f(Z,t)=\beta_{2}Z\mathrm{log(t)}$, we set $\beta_{1}$ = -0.69, $\beta_{2}$ = 0.25, with a sample size of 5000 total cases. For each individual, survival times under $Z=1$ and $Z=0$ were both drawn using the R package simsurv. [37] These true potential outcomes were used to calculate the true ATE for the estimands of interest. Observed treatment status was drawn from the binomial based on each individuals probability of receiving the treatment, as defined by our propensity score model. Covariate effect sizes ranged from 0 to 0.8 (see Appendix Table 1).

For each set of simulated outcomes, we first estimated the propensity scores using a logistic regression model, and then calculated respective IPTWs for each individual. Using these estimated weights, we fit each of the models described above in Section 3 using standard R packages, including flexsurv for parametric survival models and pseudo for pseudo-observations methods.

We repeated this simulation framework with various changes to Equation 2, for a total of five scenarios:

1. 1.

   Base case: $\beta_{1}$ = -0.69 ; $f(Z,t)=0.25Z\mathrm{log}(t)$

2. 2.

   Piecewise constant (PWC) hazard ratio: $\beta_{1}$ = 0; $f(Z,t)=-0.25Z\mathrm{I}(t\geq 2)$, where $\mathrm{I}(\cdot)$ is the indicator function

3. 3.

   Modest non-proportionality $\beta_{1}$ = -0.69 ; $f(Z,t)=0.125Z\mathrm{log}(t)$

4. 4.

   Modest treatment effect $\beta_{1}$ = -0.41 ; $f(Z,t)=0.25Z\mathrm{log}(t)$

5. 5.

   Base case, with a smaller sample size (N=500 per arm)

Senario 1 (base case) had a large treatment effect with substantial non-proportionality. Scenarios 2-4 were chosen to approximate the parameters and non-proportionality found in the clinical examples, and scenario 5 was used to assess how the methods perform with a more modest sample size.

For each method, we assessed performance based on the mean bias, standard error, and coverage for the ATE for treatment versus the control. For the base case, the bias was often largest for the standard Cox model, as expected given the substantial degree of non-proportionality in the simulated data. For example, considering median survival, bias for the Cox model was -0.347 years, while all other methods had biases less than -0.069. The generalized gamma AFT model also had large bias relative to other methods; at 10 years, it’s bias was actually greater than that of the Cox model (2.5% vs 2.0%, respectively). Even though this three-parameter model is flexible and allows for non-proportionality, fitting this mis-specified parametric model still introduced substantial biases into the results, often comparable in size to those found using the standard Cox model. In the base case, the Cox model with a treatment effect varying by log-time was correctly specified, as was the fully parametric AFT model with variable location and scale parameters. These models had lower biases than the standard Cox model (e.g. for median survival, biases were -0.042 and -0.056, respectively). However, the non-parametric Kaplan-Meier and pseudo-observation methods generally had the lowest biases (e.g. biases for median survival were -0.021 and -0.020, respectively). Finally, we note that the RMS was less sensitive to the choice of method than the other outcomes of interest; the relative differences in the sizes of the biases were smallest for this survival outcome. (see Figure 1)

As shown in Table 1, coverage was worst for the standard Cox model, followed by the generalized gamma AFT model. The time-varying Cox models did better, although coverage was surprisingly low for the log-time model at the 2-year outcome (0.82). Coverage for the non-parametric methods was close to the nominal level for all outcomes. Standard Errors (SEs) of the differences in treatment effects were usually slightly larger for the non-parametric methods than for the parametric and semi-parametric methods. (See Appendix Figure 1) Comparing the SEs for the Kaplan-Meier versus the Weibull model, the increases ranged from 0.3-24.4%, and were lowest for RMS and highest for 2-year survival.

For piecewise constant hazards (scenario 2), biases were typically lowest for the non-parametric methods and when using a PWC time-varying Cox model. An exception to this was for RMS, where the naïve Cox model and Weibull AFT model had the lowest bias. This likely occurred because bias from different timepoints were in different directions (positive at 2 and 5 years, negative at 10 years), and may have cancelled each other out when calculating the area under the survival curve. (Note that Figure 1 shows the absolute value of the average biases). So, although biases may have been lower in this simulation for these mis-specified models, such behavior should not be relied upon in other scenarios. Coverages (See Appendix Table 2) were close to the nominal when biases were low. Standard errors were generally similar to be base case, except for median survival, where they were somewhat larger.

When the effect of non-proportionality was modest (scenario 3), biases tended to be smaller than those of the base case, especially for estimates of survival probabilities at a given timepoint. When the treatment effect was modest (scenario 4), results were generally similar to those of the base case. In the small sample size case (scenario 5), the bias was similar to that of the base case, with the exception of RMS, where across methods, biases were larger than those of the base case. As expected, SEs were higher when the sample size was smaller.

The computational time and resources required to fit each model varied substantially. Depending on the size of the dataset, computational considerations may become important, especially given our use of the bootstrap to compute confidence intervals. The fastest method was weighted KM. To fit a Cox model with a time varying-effect of log-time, we found it infeasibly slow to split the dataset at each observed failure time given our large sample size. Instead, we limited the splits to occur at one-month intervals; we tested this approach and found that the results were nearly identical to those obtained by splitting at each failure time. The pseudo observations approach was also relatively slow, and highly memory intensive, even when we estimated $S(t)$ at one-month intervals. Table 2 shows the time required to fit each model once, and the total time required to calculate 500 bootstrap samples.

Soft-tissue sarcoma is a rare cancer of the connective tissue (e.g. fat, muscle, blood vessel), usually arising in the extremities. Stage III disease (high grade, large, or deep tumor) is typically treated surgically or with radiation therapy. It remains uncertain whether chemotherapy provides an overall survival benefit for these patients, and its use is optional according to national guidelines. [38] We studied this question using the NCDB, creating a cohort of Stage III sarcoma patients (with various histologies) treated with definitive surgical resection of the primary tumor, with or without chemotherapy pre/post operatively. We identified 5,337 cases with recorded overall survival data, of whom 28% were treated with chemotherapy. [8]

Due to the non-randomized nature of the data, there were substantial differences between treatment groups. For example, patients treated with chemotherapy were younger and had fewer recorded co-morbid conditions. Furthermore, the hazards for the treatment groups exhibited substantial non-proportionality, as is clearly evident upon examining the complementary log-log survival curves. (See Figure 2) Testing the Schoenfeld residuals provides further evidence of non-proportionality (P=0.047). We used IPTW to account for covariate imbalance between the treatment groups. The propensity score was estimated using all available covariates (including age, gender, insurance status, income, comorbidity score, tumor histology, grade, size, anatomic site, treating facility type, and travel distance.) The IPTW Kaplan-Meier curves are shown in Figure 2. In our original analysis, we also fit a time-varying Cox proportional hazards model, assuming a piecewise constant hazard ratio, which could change at two years. We found that there was little effect of chemotherapy during the first two years after diagnosis, but that patients treated with chemotherapy did better in the long term. These results motivated simulation scenario 2 (PWC).

Here, we apply each of the methods discussed in the sections above to the sarcoma data. Results for RMS, 2 year survival, and 5 year survival are shown in Table 3. Median and 10-year survival were not estimable from the observed data, so we excluded these outcomes. We found that the different methods yielded different effect size estimates; however, these differences were often modest and in some cases would not change statistical inferences. For example, $\Delta$RMS varied from a minimum of 0.409 (AFT with a generalized gamma distribution) to a maximum of 0.468, but all models showed an improvement in RMS for patients treated with chemotherapy. In other cases the inferences would differ based on the model chosen. At 2 years, the standard Cox model showed a difference of 3.6% between the treated and control arms, which was statistically significant. The results of both AFT models also found significantly higher survival in the chemotherapy group. However, the weighted Kaplan-Meier method showed a non-significant difference of -0.7%. Both time-varying Cox models and the pseudo-observations method also produced non-significant differences.

We would like to note that in this data analysis example, we found that the pseudo-observations method was sensitive to deviations in the mean of the estimated weights from 1. While the expectation of the mean of the weights is equal to 1, in practice there may be small deviations from the expected value. The pseudo-observations method seemed to be more affected by these deviations than the other methods we evaluated. This sensitivity can be explained by how the weights are applied and how the survival function is estimated. Each patient’s (binary) contribution to the survival function is weighted, and the weighted indicators are then averaged. Therefore, if the mean of the weights is not equal to 1, at baseline the survival function will not start at 1. To address this issue in the estimation of the survival function when using pseudo-observations, we centered the weights so their mean was exactly 1. However, we note that this was done heuristically, and not theoretically driven. Further exploration of this issue is warranted.

The second motivating clinical example for this work was a study of surgical options for early stage renal cancer. These small tumors can be treated with Radical Nephrectomy (RN) or Partial Nephrectomy (PN). PN preserves more kidney tissue, but surgery/recovery is more difficult and it is more likely to be given to healthier patients. [39, 9] Here, we focus on a subgroup from the original study: patients aged 51-60 with T1a tumors, which gives a sample size of 28,973 patients (61.1% treated with PN). In unadjusted analysis, patients given PN had improved survival over those given RN. However, these results are subject to confounding by indication. Examination of the complementary log-log plots shows a notable deviation from non-proportionality, which is confirmed by testing the Schoenfeld residuals (P$<$0.0001). [26]

When we fit an IPTW Cox PH model to these data, allowing the hazards to vary as a function of log-time, such that $h(t)=\lambda_{0}(t)\mathrm{exp}(\beta_{1}Z+\beta_{2}Z\mathrm{log}(t))$ the maximum likelihood estimates are $\hat{\beta_{1}}=-0.468$ and $\hat{\beta_{2}}=0.136$. Therefore, the non-proportionality observed here is close to our modest time-varying hazard simulation scenario (number 3), and the treatment effect is modest as well, similar to what we used in scenario 4.

Our results were largely simiar, regardless of what model was used. (See Table 4) This was particularly the case for the RMS and five-year survival. Two-year survival showed the largest range in results, with estimated differences ranging from 0.9% (Generalized Gamma AFT model) to 1.6% (Cox model with piece-wise time-varying effects). The standard Cox model gave an estimated difference of 1%, compared to 1.5% from the Kaplan-Meier method. Even with these small effects, we found that none of the 95% confidence limits crossed zero, so all models would lead to the same inference: a small benefit for PN over RN. The congruous inferences were partially driven by the very large sample size, and the resulting small confidence intervals; however, given the clear evidence of non-proportionality, it is notable that the point estimates for the standard Cox model were close to those of time-varying or non-parametric methods.