**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](http://gforge.se/2016/03/dealing-with-non-proportional-hazards-in-r/) and in [my answer](https://stats.stackexchange.com/questions/18576/time-dependent-coefficients-in-r-how-to-do-it) to my own question.