Cox models assume proportionality. These assumptions can often be formally tested (an example in R is cox.zph from the survival package).
Is there any statistical method for seeking out violations of the proportional hazards assumptions for a Cox models that are fitted with an ElasticNet-style procedure?
You need to extract scaled Schoenfeld residuals from a penalized (via ridge, LASSO, or elastic net) Cox model returned, say, by the glmnet() function. The problem is that the object returned by the glmnet() function isn't itself a Cox model; it just contains the set of penalized coefficients for such a model. This also poses problems for predictions from such models.
Thus you have to specify a Cox model that uses the set of penalized coefficients returned by the elastic net. You can do that (without undoing the elastic-net penalization) by forcing the Cox model software to start with those coefficient values while preventing it from finding the maximum partial likelihood solution. See this answer, for example. Tests on coefficients won't work, as you have prevented the model from finding the optimal solution. You can, however, extract scaled Schoenfeld residuals.
Here's an example with LASSO, using a case from the glmnet() help page. The same approach applies to any situation in which you need to specify a fixed set of Cox regression coefficients.
library(glmnet)
library(survival)
set.seed(10101)
N = 1000
p = 30
nzc = p/3
x = matrix(rnorm(N * p), N, p)
beta = rnorm(nzc)
fx = x[, seq(nzc)] %*% beta/3
hx = exp(fx)
ty = rexp(N, hx)
tcens = rbinom(n = N, prob = 0.3, size = 1) # censoring indicator
y = cbind(time = ty, status = 1 - tcens) # y=Surv(ty,1-tcens) with library(survival)
## above is from glmnet() help page
## find lambda min for LASSO by cross validation, get fit there
cv.fit <- cv.glmnet(x, y, family = "cox")
bestlam <- cv.fit$lambda.min
fitAtLmin <- glmnet(x,y,family="cox",lambda=bestlam)
## fix those coefficient values into a Cox model
coxL <- coxph(Surv(ty,1-tcens)~x,
init=as.numeric(coef(fitAtLmin)), ## specify coefficient values
iter.max=0) ## force the software to keep those values
cox.zph(coxL)
# chisq df p
# x 41.6 30 0.077
# GLOBAL 41.6 30 0.077
With a matrix provided as the predictor to coxph() as it was to glmnet(), you only get a single "predictor" evaluated by cox.zph() even though you have regression coefficients (the values you specified) for all the individual predictors. If you want to evaluate the predictors individually, provide them individually to coxph(). The GLOBAL test is the same either way (not shown).
A couple of warnings. First, I'm not sure how reliable this evaluation of proportional hazards (PH) is when coefficients are known to be biased, as they are in penalized regressions. I think that trends over time should be visualized OK, as the plot of residuals for each predictor is based on the corresponding fixed-in-time model coefficient. Second, as elastic net models typically involve large numbers of predictors and large numbers of cases, you might have "significant" violations of PH that don't matter in practice. You thus have to apply your knowledge of the subject matter to evaluate the results.
