Test Cox proportional hazard assumption (Bad Schoenfeld residuals)

Question

Using R I generated a Cox model looking like this

> summary(cox5 <- coxph(surv_ob_2 ~ age + AMP + nation + GESL + year, data=sd)) # now gender makes sense
Call:
coxph(formula = surv_ob_2 ~ age + AMP + nation + GESL + year, 
    data = sd)

  n= 35800, number of events= 31873 

                             coef  exp(coef)   se(coef)       z             Pr(>|z|)    
age                    -0.0136536  0.9864392  0.0005184 -26.336 < 0.0000000000000002 ***
AMPofficial            -0.6058868  0.5455904  0.0338201 -17.915 < 0.0000000000000002 ***
AMPmarginal employment  0.5032381  1.6540687  0.0155578  32.346 < 0.0000000000000002 ***
AMPfreelancer          -0.7277312  0.4830036  0.0342688 -21.236 < 0.0000000000000002 ***
nationFormerYugoslavia  0.2758670  1.3176726  0.0304818   9.050 < 0.0000000000000002 ***
nationGermany           0.2354548  1.2654842  0.0341056   6.904     0.00000000000507 ***
nationother             0.1421815  1.1527859  0.0144441   9.844 < 0.0000000000000002 ***
GESLM                   0.1089281  1.1150822  0.0114915   9.479 < 0.0000000000000002 ***
year                    0.0526784  1.0540906  0.0009662  54.521 < 0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                       exp(coef) exp(-coef) lower .95 upper .95
age                       0.9864     1.0137    0.9854    0.9874
AMPofficial               0.5456     1.8329    0.5106    0.5830
AMPmarginal employment    1.6541     0.6046    1.6044    1.7053
AMPfreelancer             0.4830     2.0704    0.4516    0.5166
nationFormerYugoslavia    1.3177     0.7589    1.2413    1.3988
nationGermany             1.2655     0.7902    1.1837    1.3530
nationother               1.1528     0.8675    1.1206    1.1859
GESLM                     1.1151     0.8968    1.0902    1.1405
year                      1.0541     0.9487    1.0521    1.0561

Concordance= 0.623  (se = 0.002 )
Rsquare= 0.17   (max possible= 1 )
Likelihood ratio test= 6690  on 9 df,   p=0
Wald test            = 6141  on 9 df,   p=0
Score (logrank) test = 6285  on 9 df,   p=0

I now want to test if the proportional hazard assumption is met. Therefore I run

> print(test.cox5 <- cox.zph(cox5, transform=rank)) ## Gender does not meet ph-assumption
                           rho   chisq             p
age                     0.1099  435.85 0.00000000000
AMPofficial             0.1001  312.68 0.00000000000
AMPmarginal employment  0.0169    9.26 0.00233833597
AMPfreelancer           0.0332   34.78 0.00000000369
nationFormerYugoslavia  0.0171    9.31 0.00227422735
nationGermany           0.0105    3.50 0.06133812325
nationother            -0.0226   16.14 0.00005894624
GESLM                   0.0157    7.88 0.00500199726
year                    0.1047  226.84 0.00000000000
GLOBAL                      NA 1193.78 0.00000000000

These results look very bad to me. Do I have to assume that only nationGermany (hardly) meets the ph-assumption, or is there another explanation as well?

My teacher supposed that a possible reason could be the large sample size (n=35800). Is he right?

What else can I do to check if the ph-assumption is met and to improve the model?

If you know any literature covering this issue, I would be interested as well.

[edit]
The data I analyse comes from a register based labor market database. I try to model the length of employments. The database covers several decades, therefore also macroeconomic circumstances could have changed with time. This is why I added age and year in the model. I am however unsure if this is a good idea...

Answer 1

It is likely that the large sample size is responsible for the seemingly strong evidence against the PH assumption. P-values are a function of sample size, and their usefulness declines when sample size grows very large as the null hypothesis is never exactly true. They don't help too much with your question here, which is not "is the PH assumption satisfied" but "is the deviation from the PH assumption so large that inference is impaired".

One way to assess this for categorical variables, which your model seems to mostly contain, is by a log-minus-log plot. This is explained e.g. in this book and easily implemented in R using the rms library

library(rms)

myfit = survfit(Surv(time, status) ~ catvar)
survplot(myfit, loglog=T)

When the PH assumption holds, the lines are parallel, and their vertical distance is the log hazard ratio.

Convergent curves are seen when the difference between the groups decreases with time, and divergent curves when it increases, which indicates some deviation from the PH assumption. Crossing of the curves indicates a more severe deviation, with the effect of group membership changing signs.