Since the seasonal dummy variables are static by nature, and their coefficients clearly vary with the time variable, how much does it matter?
The value that you get is a form of average over time. Unfortunately as you naturally have more cases in time-to-event analyses early on, you can't simply say that the effect is balanced throughout time. From my experience the initial period has a much heavier impact on the estimate than the later and it
I get that statistically, it means something, but does the violation of PH assumption invalidate the (intuitively appealing) result that non-response is more likely to happen in the summer and winter?
Again, it probably doesn't but you can't be sure until you've checked. Something is definitely happening over time and you should at least have a look at the residual plot (plot(cox.zph(...))). It isn't entirely surprising that you have a problem with the PH since seasons are part of the time variable and there will be situations where early summer and late spring are similar.
If so, is there a way to handle this so that the PH assumption is not violated? I know about using the tt transform, but I can't seem to figure out the exact form for the function.
The tt transform is tricky to use with big data. It explodes the matrix and can get a little messy, e.g. if you modify the lung example in the survival package:
library(survival)
coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,
tt=function(x,t,...) {
print(length(x))
pspline(x + t/365.25)
})
It prints 15809 while there are only 228 rows in the original dataset. The principle of the tt() is that it feeds the variables into the transformation function where you are free to use time any way you wish. Note that you can also have different transformation functions depending for each variable:
library(survival)
coxph(Surv(time, status) ~ tt(ph.ecog) + tt(age), data=lung,
tt=list(
function(x,t,...) {
cbind(x, x + t/365.25, (x + t/365.25)^2)
},
function(x,t,...) {
pspline(x + t/365.25)
}),
x=T,
y=T) -> fit
head(fit$x)
Gives:
tt(ph.ecog)x tt(ph.ecog) tt(ph.ecog) tt(age)1 tt(age)2 tt(age)3 tt(age)4
6 1 3.4 11.7 0 0 0 0.000
3 0 2.4 5.8 0 0 0 0.020
38 1 3.4 11.7 0 0 0 0.000
5 0 2.4 5.8 0 0 0 0.000
6.1 1 3.2 10.4 0 0 0 0.000
3.1 0 2.2 5.0 0 0 0 0.026
tt(age)5 tt(age)6 tt(age)7 tt(age)8 tt(age)9 tt(age)10 tt(age)11 tt(age)12
6 0.00 0.00000 0.000 0.0052 0.359 0.58 0.053 0
3 0.48 0.48232 0.021 0.0000 0.000 0.00 0.000 0
38 0.00 0.00087 0.266 0.6393 0.094 0.00 0.000 0
5 0.03 0.51933 0.437 0.0136 0.000 0.00 0.000 0
6.1 0.00 0.00000 0.000 0.0078 0.388 0.56 0.044 0
3.1 0.50 0.45457 0.016 0.0000 0.000 0.00 0.000 0
I therefore try to avoid this solution and use the time-split approach that I wrote about here and in my answer to my own question.
One issue here is your choice of reference level for the Group variable, G1. The regression coefficients for other Groups are with respect to that reference level, and as I understand it the same is true for the "significant" non-proportionalities seen for the other Groups. Note that this type of summary does not provide a test for the significance of the Group variable as a whole. Had you chosen a different reference level, much of the difficulty with non-proportionality might have been isolated to just one or two Groups. It's important to think about the subject-matter content of your data; there might be good reasons why some groups have different hazard time courses than others.
Also, be careful about how you interpret the p-values for the cox.zph tests. A low p-value for a coefficient is evidence that the proportional hazards assumption doesn't hold for its associated predictor, but a "non-significant" p-value is not proof that the proportional hazards assumption is met. As with any statistical test, a non-significant p-value might simply mean two few cases or too much variability to argue against the null hypothesis of a proportional hazard. It's hard to tell from your graph, but that might explain why crossing plots have p-values that do not rule out the PH assumption.
Best Answer
If the covariate is very influential you sometimes benefit by adjusting for it as a covariate and also as a stratification factor. Stratification will not hurt the effective sample size very much for estimating relative hazards, but will lower the precision of survival curve estimates. But you need to state how you avoided linearity assumptions for continuous covariates. Frequently, nonlinearities can be more impactful than non-proportional hazards.
Sometimes when you change models because of a failure to satisfy proportional hazards for one variable, you make the fit worse for other variables. It can be helpful to get a sense of the overall "average" structure. In RMS I have a detailed case study in a latter chapter where I examine the fit of a parametric survival model, and used an accelerated failure time model instead of the Cox model.