Testing the proportional hazards assumption with a time varying covariate

Question

One of the cohort studies I have carried out recently involved examining the association between post-diagnostic beta-blocker use and breast cancer outcomes in a cohort of breast cancer patients. Briefly, I followed patients up from their breast cancer diagnosis until death or end of 2017, and treated beta-blocker use as a time-varying covariate. That is, people were considered nonusers of beta blockers until their first post-diagnostic beta blocker dispensing, and once they were a user they were a user until the end of follow-up. For people who used beta blockers, they have two rows of data in my STATA dataset, with the rows split based on when they first used the beta blocker. An example of this data is shown below:

In this dataset, 'time5' represents each patient's entire follow-up time, and 'timetofbetablocker' is the time between breast cancer diagnosis and the patient's first beta blocker dispensing. As you can see, I've split the data for people who used BBs (patient 476 is the only one to have used a BB in this sample data), so in row 9 patient 476 is counted as using a BB between 8.85 years and 10.82 years. I've then created a variable (the second to last column in this sample data) to run a Cox model on, representing the time-varying nature of BB use (code if for STATA, sorry R users):

gen timetofbetablockerana=1 if _t0>=timetofbetablocker & timetofbetablocker!=.

replace timetofbetablockerana=0 if timetofbetablockerana==.

I can then run a simple (unadjusted) Cox model for BBs using stcox timetofbetablockerana

A question I have been grappling with for a while now is how I adequately test the PH assumption with this time-varying covariate. I don't think a simple estat phtest (testing the proportional-hazards assumption on the basis of Schoenfeld residuals) in STATA suffices for this, because the variable is time-varying. In slide 32 of this presentation (https://ms.uky.edu/~mai/sta635/Cox%20model.pdf), they suggest interacting the time-varying covariate with time itself and then testing the assumption on this variable. So, I was thinking of interacting 'timetofbetablockerana' with 'timetofbetablocker' to allow the value of the analysis variable to vary based on when people started using BBs, and then including this interaction term in the model. e.g, gen interaction= timetofbetablocker* timetofbetablockerana, replace interaction=0 if interaction==., and then stcox timetofbetablockerana interaction and estat phtest, detail to test the assumption on the interaction term.

Is anyone able to let me know if I'm on the right track here? If this isn't correct, does anyone have any other ideas/suggestions on how I could go about this? Further, even with this time-varying covariate in the model, if I were to test the PH assumption on time-independent/fixed covariates (which are additionally included in the fully adjusted model), would I simply go about it the usual way, even in the presence of time-varying covariates in the same model?

Answer 1

I don't know the origin of the idea that Schoenfeld residuals can't be used for testing proportional hazards (PH) with time-varying covariate values. As Therneau and Grambsch say in Section 4.6:

Since the [Schoenfeld] residuals are defined at each unique death, their definition and computation is unchanged by a counting process formulation of the data set.

The "counting process formulation of the data set" is precisely what's used for time-varying covariates. The Cox model simply uses whatever covariate values are currently in place for all individuals at risk at each event time. Whether a value for some individual might be different from some prior value doesn't matter (provided that there's no more than 1 event possible per individual).

That fact, however, also points out the difficulty in formulating a useful time-varying covariate. Since only the current value of a covariate at an event time is used, the covariate value used for the model must incorporate whatever aspect of covariate history is associated with outcome at any particular time.

I don't know whether the STATA PH test based on Schoenfeld residuals might have some issues with time-varying covariates, but there's not a theoretical problem.

In response to/ incorporating comments

The time dependence vignette for the base R survival package deals (explicitly or implicitly) with most of the issues raised in comments. I think that all contained in this answer and my comments comes pretty directly from there or from the classic Therneau and Grambsch text linked from the first paragraph.

There is no problem with evaluating scaled Schoenfeld residuals for binary or categorical predictors, or for a mix of time-constant and time-varying covariates.

A single yes/no time-varying exposure variable as you use for the beta-blockers assumes that the association of outcome with drug exposure is independent of the duration of exposure: that the very first day on the drug provides the same effect as being on it for 10 years. A rolling-average estimate of drug exposure over a biologically reasonable time span might be better. In practice, a 0/1 no-drug/drug covariate might work well enough.

Constructing a useful time-varying covariate and testing it for PH are two different things. Whether you use a 0/1 drug exposure indicator or some type of running average, the PH test will determine whether the association of that marker with outcome is independent of time since diagnosis. PH might be violated, for example, if the drug made a lot of difference soon after diagnosis but had less association with outcome later on.

I don't know what's available in STATA for testing the PH assumption via scaled Schoenfeld residuals. The vignette linked above describes what's available in the basic R survival package. Numeric tests for PH depend on the choice of time scale, something that isn't always obvious at first.

There are several options for time-scale transformation available in R; transformation based on Kaplan-Meier estimates is a common choice. What's shown in your graph suggests that STATA is using a logarithmic transformation of time instead. In that time scale, at least, it looks like PH is a reasonable assumption, given the p-value that you report and the corresponding flat estimate of what I take to be a smoothed representation of the scaled residuals (see the "bandwidth" notation near the bottom of the graph).

The R survival package does have an overall display and test of the PH assumption for multi-level categorical predictors. I don't know whether that's available in STATA. If a categorical predictor sufficiently violates PH and isn't of interest except as something for which you wish to control, stratification on that predictor can be a useful choice. That said, the PH assumption is less critical for such control predictors; see the discussion on this page.

For a predictor of primary interest that fails PH, the vignette notes that partitioning of survival time into a few blocks can sometimes allow for PH within each time block, with a separate HR for each time block. Like any categorization of a continuous variable, however, that is somewhat arbitrary at best.