Survival analysis with time dependent covariates and non-proportional hazards in R

Question

I am attempting to do a survival analysis which will examine the effects of both rainfall (a time-dependant variable) and altitude on nest survival in a species of wasp found in NW Ecuador.

I have relatively large dataset where both the presence/absence of nests across an altudinal transect were monitored periodically every 5 days. Rainfall data were also collected every 5 days from 9 different raingauges along the same gradient.

In order to obtain rainfall estimations for specific nests, I have modelled the rainfall data using a generalized additive model (GAM) with both altitude, date checked (the Date the rainfall measurement was taken) and their interaction as explanatory variables. Using the predict function I've obtained rainfall estimations for each individual nest at each time checked.

Here is a sample of the resulting data in long format:

Nest.No   Site    Altitude  Date.Found   End.Date   Days.Old Censored Check.No Time.Start Time.Stop Est.Daily.RF
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          20      23      28           50
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          16       3      8            43
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          25      48      53           32
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          26      53      58           30
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          21      28      33           28
  625   Mashpi      1062    18/01/2017  20/03/2017      63       1          27      58      63           25
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          18      13      18           24
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          23      38      43           18
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          24      43      48           17
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          19      18      23            8
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          17       8      13            8
  625   Mashpi      1062    18/01/2017  20/03/2017      63       0          22      33      38            7
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          20      17      22           49
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          21      22      27           30
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          25      42      47           30
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          26      47      52           26
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          18       7      12           25
  674   Mashpi      1129    24/01/2017  20/03/2017      57       1          27      52      57           24
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          24      37      42           19
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          23      32      37           18
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          16       0       2           43
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          19      12      17            8
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          17       2       7            7
  674   Mashpi      1129    24/01/2017  20/03/2017      57       0          22      27      32            6
  800   Mashpi      1180    18/02/2017  20/03/2017      32       0          25      17      22           28
  800   Mashpi      1180    18/02/2017  20/03/2017      32       0          26      22      27           24
  800   Mashpi      1180    18/02/2017  20/03/2017      32       1          27      27      32           22
  800   Mashpi      1180    18/02/2017  20/03/2017      32       0          24      12      17           20
  800   Mashpi      1180    18/02/2017  20/03/2017      32       0          23       7      12           18
  800   Mashpi      1180    18/02/2017  20/03/2017      32       0          21       0       2           30
  800   Mashpi      1180    18/02/2017  20/03/2017      32       0          22       2       7            7
  801   Mashpi      1173    18/02/2017  20/03/2017      32       0          25      17      22           28
  801   Mashpi      1173    18/02/2017  20/03/2017      32       0          26      22      27           25
  801   Mashpi      1173    18/02/2017  20/03/2017      32       1          27      27      32           22
  801   Mashpi      1173    18/02/2017  20/03/2017      32       0          24      12      17           20
  801   Mashpi      1173    18/02/2017  20/03/2017      32       0          23       7      12           18
  801   Mashpi      1173    18/02/2017  20/03/2017      32       0          21       0       2           30
  801   Mashpi      1173    18/02/2017  20/03/2017      32       0          22       2       7            7
  827   Mashpi      1201    23/02/2017  20/03/2017      27       0          25      12      17           28
  827   Mashpi      1201    23/02/2017  20/03/2017      27       0          26      17      22           24
  827   Mashpi      1201    23/02/2017  20/03/2017      27       1          27      22      27           22
  827   Mashpi      1201    23/02/2017  20/03/2017      27       0          24       7      12           21
  827   Mashpi      1201    23/02/2017  20/03/2017      27       0          23       2       7           18
  827   Mashpi      1201    23/02/2017  20/03/2017      27       0          22       0       2            7
  828   Mashpi      1107    23/02/2017  20/03/2017      27       0          25      12      17           31
  828   Mashpi      1107    23/02/2017  20/03/2017      27       0          26      17      22           28
  828   Mashpi      1107    23/02/2017  20/03/2017      27       1          27      22      27           25
  828   Mashpi      1107    23/02/2017  20/03/2017      27       0          24       7      12           18
  828   Mashpi      1107    23/02/2017  20/03/2017      27       0          23       2       7           18
  828   Mashpi      1107    23/02/2017  20/03/2017      27       0          22       0       2            6

I have fitted a Cox's Proportional Hazards Regression Model to the data using the code:

cfit<-coxph(Surv(Time.Start,Time.Stop,Censored)~Est.Daily.RF+Altitude, data=Mashpi)

The resulting call

Call:
coxph(formula = Surv(Time.Start, Time.Stop, Censored) ~ Est.Daily.RF + 
Altitude, data = Mashpi)

               coef      exp(coef)  se(coef)   z     p
 Est.Daily.RF  0.012874  1.012957  0.005474  2.35 0.019
 Altitude     -0.000393  0.999607  0.000650 -0.60 0.545

 Likelihood ratio test=5.86  on 2 df, p=0.0533
 n= 2030, number of events= 245 

So I can see that I have a significant effect of estimated daily rainfall on nest survival. However I have several questions which I haven't been able to find clear answers to. I'm posting here with the hope that somebody more experienced in these types of analyses might be able to give me some pointers.

Q1. My data includes some nests which are both right and left-censored. By this I mean that for some of the individuals nests I do not know the original founding date. For these nests "Time" starts on the day the nest was found rather than the date founded. For other nests which were founded during my field season I know the exact founding date but for many of these I do not know the exact end date. The right-censored data seems to be accounted for in my model but the left-censored data does not, is it necessary or even possible to include this information in my analysis?

Would one potential option to include a variable in my model which indicated whether or not the nest was newly founded when entering into the observations?

Q2. When I test the proportional hazards assumption for a Cox regression model fit using the cox.zph(cfit), I obtain the following results:

                rho   chisq  p
 Est.Daily.RF  0.0465 0.539 0.4627
 Altitude      0.1411 4.718 0.0299
 GLOBAL        NA     5.162 0.0757

So altitude appears to be violating this assumption but rainfall does not. As altitude does not seem to be having a significant result on survival do I need to adjust my model to account for the non-proportional risks associated with altitude?

I understand that to correct for this non-proportionality I can include an interaction term between Altitude and Start.Time. However, if I do this and re-test the propoprtional hazards assumption all the variables now appear to be violating this assumption

vfit<-coxph(Surv(Time.Start,Time.Stop,Censored)~Est.Daily.RF+Altitude+Altitude:Time.Start, data=Mashpi)

zp.M<- cox.zph(vfit) 

                     rho   chisq   p
Est.Daily.RF        0.0665  1.03   3.11e-01
Altitude            0.3470 33.68   6.51e-09
Altitude:Time.Start 0.1126  4.15   4.17e-02
GLOBAL                  NA 88.84 3.88e-19

Moreover, as some of my nests are left-censored, not all of the nests are the same age for corresponding time points. i.e a nest which is left-censored may be much older for time 0-4 that a new nest. Am I right in thinking that this means I'm unlikely to detect any non-proportional risks associated with rainfall?

Q4. As mentioned, my individual Est.Daily.Rainfall datapoints have been generated using a model based the highly significant interaction effect of Check.Number:Altitude to make its predictions. This variable is therefore highly correlated by Check.Number. Would it be correct to include cluster(Check.Number) in my coxph model? In the examples I have seen this function is used to cluster individuals into groups (eg clinic) rather than to cluster individual observations of a time-dependent variable.

Note: doing so gives me highly significant results for all varibles and also removes non-proportionality from the model.

Q3. I'm aware that this model is only appropriate for time-dependent variables of an exogenous nature. However I'm now starting to think that as rainfall is unpredictable, contains error in my estimations and, for many, the complete rainfall path is not fully observed, this varaible might actually be endogenous. Would a joint model be more appropriate?

Q4. If so, how I might go about using a joint modelling approach? Is it as simple as fitting a linear model with rainfall~time and including this with my original Cox's Proportional Hazards Regression Model in a joint model?

This is the first time I've done any type of survival analysis and I'm feeling very unsure of the correct way to proceed. Any help would be greatly appreciated.

Answer 1

I've obtained rainfall estimations for each individual nest at each time checked

Remember that a proportional hazards model with time-varying covariates only uses the instantaneous covariate values at event times. So you are implicitly modeling nest survival as a function of that particular day's rainfall, not as a function of some measure of recent cumulative rainfall (which to me as a non-expert would seem to make more sense). Make sure you have decided just what you want to be modeling.

My data includes some nests which are both right and left-censored.

It's important to be clear about what is censored. It's not nests, it's survival times. For nests that already existed at the start of the study, you still have a lower limit to the total nest survival time when the nest is gone. That's right censoring of their survival times. So instead of marking an event at the time such a nest is gone, keep that value marked as censored. That will, however, decrease the number of events and thus the power of your model.

altitude appears to be violating this [PH] assumption but rainfall does not

You have to apply your knowledge of the subject matter and perhaps modify your model. If you expected based on your knowledge of the subject matter that altitude would be an important predictor, it wouldn't be right just to drop it from the model. Sometimes an apparent lack of PH for a continuous predictor comes from improper specification of its association with outcome. A flexible regression spline fit for altitude might both show something interesting and fix the PH problem. Alternatively, you could evaluate the plots of scaled Schoenfeld residuals to consider whether the magnitude of the lack of PH is large enough to matter. That's a judgment call.

Would it be correct to include cluster(Check.Number) in my coxph model?

Maybe. If observations aren't independent then intra-cluster correlations need to be taken into account somehow. I haven't yet thought through how well this would work when the issue is correlated errors in predictor variables. A Cox model also assumes no error in the predictor variables. Your next question indicates what might be a better approach.

Would a joint model be more appropriate?

A joint model that combines survival analysis with longitudinal modeling of covariates is an accepted way to deal with correlations and with covariates that are missing or measured with some error. See this paper from a quarter-century ago and this recent tutorial. The CRAN survival view provides links to the joineR, JointModel, and JM packages. I don't have experience with such modeling, though.

Answer 2

I am not an expert on survival analysis but I am also working on a similar problem now and here's my take on your situation.

Q1. My data includes some nests which are both right and left-censored. By this I mean that for some of the individuals nests I do not know the original founding date. For these nests "Time" starts on the day the nest was found rather than the date founded. For other nests which were founded during my field season I know the exact founding date but for many of these I do not know the exact end date. The right-censored data seems to be accounted for in my model but the left-censored data does not, is it necessary or even possible to include this information in my analysis?

My understanding is that if you have both left and right censored data, you need to specify it in the type argument of the Surv function.

So altitude appears to be violating this assumption but rainfall does not. As altitude does not seem to be having a significant result on survival do I need to adjust my model to account for the non-proportional risks associated with altitude?

Since altitude does not have a significant impact on survival, you can discard it from your analysis, so it doesn't need to respect the proportional hazards assumption. See this example.