How to meet Cox proportional hazard model assumptions for seedling survival data in R

user_11780094 0
Viewed 114

I am examining seedling survival using Cox PH Regression in the Survival package in R. The survival of individual plants from 9 species have been tracked biannually for several years. Most simply, I want to 1) compare survival rates and hazard ratios between species, and later, 2) perform multivariate, mixed effect analysis based on various treatments, traits, and plots. When I run the cox model I find that my basic models do not meeting model assumptions (proportional hazards assumption). I am curious if my problem is because I need a time-dependent interaction, or if it is something simple such as data or model formatting .

My data looks like this and a subset my data is pasted further down in my post. (SpCode: species, Time: months, Event: 0 = live, 1 = dead, Assessment: assessment period, 1 = first, 6 = last):

> head(myData)
SpCode Time Event Assessment
1 A 28 0 6
2 B 28 0 6
3 B 28 0 6
4 A 28 0 6
5 C 28 0 6
6 D 28 0 6

Here is my code to run these models

myCPH <- coxph(Surv(Time, Event) ~ SpCode, data = myData) #runs model
summary(myCPH) #examine output - hazard ratios (expressed as exp(coef)) actually makes sense relative to reference species "A"
cox.zph(myCPH) # test proportionality assumption, should not be significant

I have experimented with including a time-dependent interaction but am wondering if I am applying it correctly?: myTDCPH <- coxph(Surv(Time, Event) ~ SpCode*Assessment, data = myData). Here I have used SpCode*Assessment as my interaction since Time cannot be on both the left and right sides of the formula. This approach doesn't meet model assumptions and the summary output makes even less sense. Moreover I get the following warning

Warning message: In fitter(X, Y, istrat, offset, init, control, weights = weights, : Ran out of iterations and did not converge

Here is a Dropbox link to some example data. I have included a large, 1000 individual subset of dummy data (~10% of total dataset) since I don't seem to get this error when I run this on a small dataset (<200 individuals)

In summary, my questions for the forum are

  1. Most importantly, how can I meet model assumptions given my data? Is my time-dependent interaction formatted correctly?
  2. Any chance my SpCode covariate (which contains 9 factors) could be causing a problem?
  3. Should my data be formatted differently? I have seen the survival object formatted as (TimeStart, TimeEnd, Event), but I haven't had success when I try this.

Thank you

cox-regression survival-analysis r
1 Answers

If you plot the survival curves you see that there is a big event rate in the first irval that varies across the species but then the event rates look fairly steady and roughly equal. So that's what the cox.zph is picking up: the inhomogeneity of relative rates across time:

library(ggfortify)
library(survival)
fit <- survfit(formula = Surv(Time, Event) ~ SpCode, data = myData)
autoplot(fit)

You could also have tracked this down by doing separate summary calls by species: (the distribution of Time values are clearly quite different across species categories)

by(myData, myData$SpCode, summary)
myData$SpCode: A
     SpCode         Time           Event          Assessment   
 A      :117   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 B      :  0   1st Qu.: 0.50   1st Qu.:0.0000   1st Qu.:0.000  
 C      :  0   Median : 4.00   Median :1.0000   Median :1.000  
 D      :  0   Mean   :12.52   Mean   :0.6752   Mean   :2.581  
 E      :  0   3rd Qu.:28.00   3rd Qu.:1.0000   3rd Qu.:6.000  
 F      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0                                                   
------------------------------------------------------------------------------------------------- 
myData$SpCode: B
     SpCode         Time           Event          Assessment   
 B      :114   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 A      :  0   1st Qu.:28.00   1st Qu.:0.0000   1st Qu.:6.000  
 C      :  0   Median :28.00   Median :0.0000   Median :6.000  
 D      :  0   Mean   :24.09   Mean   :0.2193   Mean   :5.088  
 E      :  0   3rd Qu.:28.00   3rd Qu.:0.0000   3rd Qu.:6.000  
 F      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0                                                   
------------------------------------------------------------------------------------------------- 
myData$SpCode: C
     SpCode         Time           Event          Assessment   
 C      :120   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 A      :  0   1st Qu.:28.00   1st Qu.:0.0000   1st Qu.:6.000  
 B      :  0   Median :28.00   Median :0.0000   Median :6.000  
 D      :  0   Mean   :26.24   Mean   :0.1083   Mean   :5.583  
 E      :  0   3rd Qu.:28.00   3rd Qu.:0.0000   3rd Qu.:6.000  
 F      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0                                                   
------------------------------------------------------------------------------------------------- 
myData$SpCode: D
     SpCode         Time           Event          Assessment   
 D      :120   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 A      :  0   1st Qu.: 4.00   1st Qu.:0.0000   1st Qu.:1.000  
 B      :  0   Median :28.00   Median :0.0000   Median :6.000  
 C      :  0   Mean   :18.78   Mean   :0.4917   Mean   :3.825  
 E      :  0   3rd Qu.:28.00   3rd Qu.:1.0000   3rd Qu.:6.000  
 F      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0                                                   
------------------------------------------------------------------------------------------------- 
myData$SpCode: E
     SpCode         Time           Event          Assessment   
 E      :240   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 A      :  0   1st Qu.: 0.50   1st Qu.:0.0000   1st Qu.:0.000  
 B      :  0   Median : 0.50   Median :1.0000   Median :0.000  
 C      :  0   Mean   :11.34   Mean   :0.6625   Mean   :2.321  
 D      :  0   3rd Qu.:28.00   3rd Qu.:1.0000   3rd Qu.:6.000  
 F      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0                                                   
------------------------------------------------------------------------------------------------- 
myData$SpCode: F
     SpCode         Time           Event          Assessment   
 F      :115   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 A      :  0   1st Qu.:12.00   1st Qu.:0.0000   1st Qu.:2.000  
 B      :  0   Median :24.00   Median :1.0000   Median :4.000  
 C      :  0   Mean   :19.27   Mean   :0.6087   Mean   :3.809  
 D      :  0   3rd Qu.:28.00   3rd Qu.:1.0000   3rd Qu.:6.000  
 E      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0                                                   
------------------------------------------------------------------------------------------------- 
myData$SpCode: G
     SpCode         Time           Event         Assessment   
 G      :122   Min.   : 0.50   Min.   :0.000   Min.   :0.000  
 A      :  0   1st Qu.: 4.00   1st Qu.:0.000   1st Qu.:1.000  
 B      :  0   Median :16.00   Median :1.000   Median :3.000  
 C      :  0   Mean   :17.19   Mean   :0.623   Mean   :3.484  
 D      :  0   3rd Qu.:28.00   3rd Qu.:1.000   3rd Qu.:6.000  
 E      :  0   Max.   :28.00   Max.   :1.000   Max.   :6.000  
 (Other):  0                                                  
------------------------------------------------------------------------------------------------- 
myData$SpCode: H
     SpCode         Time           Event          Assessment   
 H      :116   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 A      :  0   1st Qu.:28.00   1st Qu.:0.0000   1st Qu.:6.000  
 B      :  0   Median :28.00   Median :0.0000   Median :6.000  
 C      :  0   Mean   :23.35   Mean   :0.2414   Mean   :4.966  
 D      :  0   3rd Qu.:28.00   3rd Qu.:0.0000   3rd Qu.:6.000  
 E      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0                                                   
------------------------------------------------------------------------------------------------- 
myData$SpCode: I
     SpCode         Time           Event          Assessment   
 I      :118   Min.   : 0.50   Min.   :0.0000   Min.   :0.000  
 A      :  0   1st Qu.: 4.00   1st Qu.:0.0000   1st Qu.:1.000  
 B      :  0   Median :28.00   Median :0.0000   Median :6.000  
 C      :  0   Mean   :18.77   Mean   :0.4915   Mean   :3.881  
 D      :  0   3rd Qu.:28.00   3rd Qu.:1.0000   3rd Qu.:6.000  
 E      :  0   Max.   :28.00   Max.   :1.0000   Max.   :6.000  
 (Other):  0

enter image description here

Multipart questions are not welcome on SO especially when they are really methods questions. Furthermore I have no idea what you are asking whether you are applying a time dependnet interaction correctly. I see nothing that resembles such a modeling effort. Do remember that we are not looking over your shoulder. Only the Russian bots are doing that.

user_18556770
Related