I fitted a Cox PH Model and upon examination of the Schoenfeld residuals it became apparent that the proportionality assumption was violated by several variables. Following Box-Steffensmeier & Jones (2004) I included interactions with time for these covariates. However I'm not sure how to interpret this. So its obvious that this indicates time-dependency of the covariate in the sense that the effect inflates/deflates through over some function of time, but I work with sociological data and my theory indicates no time-dependency of the effects in whatever direction. So if I get that right I should therefore consider the time-dependency to come from some kind of unobserved heterogeneity? Due to the nature of the data I can also not implement frailty or fixed effects to account for this. So how do I interpret a coefficient that increases/decreases as time progresses given that the theory does not indicate this?
There are quite a few reasons why proportionality does not hold. One of them is omission of an important variable and thus you get a wrong model that may appear as time-varying coefficients. Let's have a look at a simple situation when there are two covariates $x_1$ and $x_2$. The true model is $$\lambda(t) = \lambda_0(t)\exp(x_1\beta_1+x_2\beta_2).$$
You might omit e.g. $x_2$ when modeling the data, in which case you get $$\lambda(t) = \lambda_0(t)\exp x_1\beta.$$
The misspecified model is wrong and suffers from two drawbacks.
- When $x_2$ is omitted, proportional hazards does not hold.
- It can be shown that the estimate of $\beta$ based on the misspecified model is biased as an estimate of $\beta_1$.
The other problem can be wrong functional form of one of the covariates. Diagnostics applied to the residuals from the incorrect model in this case will also suggest presence of nonproportionality. Unless your covariates change with time, you can use martingale residuals to check the functional form.
As for the interpretation, it can be challenging. Let's assume that to model nonproportionality by time-dependent covariates, you create a time-dependent covariate $X^*(t)$ so that $$\beta(t)X = \beta X^*(t).$$
Now, if $X^*(t) = tX$, then the hazard ratio is time-dependent and (assuming other variables are constant) effect of a unit increase in the variable for the hazard rate is $$HR = \exp(\beta t),$$
which you can plot on time scale for the variable $X$.
You can try different methods to deal with nonproportionality such as stratification of covariates with nonproportional effects (which might be a wrong idea for quantitative variables because of possible bias and also stratified analyses are less efficient).
Another way is to partition the time axis, in which case you can choose a time point suggested by the residual plots and fit two models for different time intervals. It also has some drawbacks because you truncate your data (subjects are sampled from conditional distribution) and the models may answer different questions that those you are interested in.
The last option is to use a different model. An accelerated failure time model might be more appropriate for your data.