lasso with *extended* cox regression (time-varying covariates using counting process notation)

Question

I'm trying to find a way to build a predictive model the development of a disease. However, some of our predictors are time-varying (aka time-dependent) -- for example, the appearance of other, age-related comorbidities.

We have a dataset with on the order of 1 million records, with about 15 million years of pt-time in follow-up

However, we have a fairly large number of potential predictors we could consider (on the order of 100 depending on how specifically we code comorbidities and concomitant medications) and the disease is quite rate (on the order of a couple of thousand events observed). Though this might suggest reasonable "events per variable" ratios, we still worry about over-fitting and would like to produce more parsimonious models using "the most important" predictors.

So, our current attempts involve penalized regression -- as in the lasso using the glmnet package in R -- using the counting process notation Surv(tstart, tstop, status) as a way to handle the time-varying covariates.

But, glmnet does not appear to be set up to handle interval censoring.

The glmnet package in R is able to address conventional cox regression, where all intervals start with an implicit time 0. However, when using the Surv(tstart, tstop, status) counting process notation, the software chokes in a way that suggests it was not designed to handle this case.

> l.fit <- glmnet(x, y, family="cox", alpha=1)
Error: Cox model requires a matrix with columns 'time' (>0) and 'status'  (binary) as a response; a 'Surv' object suffices

The documentation appears mute on this example, unless I'm missing something.

We would be grateful for suggestions on other ways to proceed, or perhaps alternative ways to specify the model such that glmnet could be used?

Thanks for any pointers or suggestions!

Rob

--

Here's the full example:

# Create Data
jasa$subject <- 1:nrow(jasa)
tdata <- with(jasa, data.frame(subject = subject, 
          futime= pmax(.5, fu.date - accept.dt), txtime= ifelse(tx.date== fu.date, 
          (tx.date -accept.dt) -.5, 
          (tx.date - accept.dt)), fustat = fustat )) sdata <- tmerge(jasa, tdata, id=subject, 
          death = event(futime, fustat), 
          trt = tdc(txtime), options= list(idname="subject")) sdata$age <- sdata$age -48 sdata$year <- as.numeric(sdata$accept.dt - as.Date("1967-10-01"))/365.25

# Extended Cox formula;
regular.fit <- coxph(Surv(tstart, tstop, death) ~ 
        age*trt + surgery + year, data= sdata)
summary(regular.fit)

# create matrix
x <- model.matrix(~age*trt + surgery + year, sdata)

# create response
y <- Surv(sdata$tstart, sdata$tstop, sdata$death)

# fit lasso
l.fit <- glmnet(x, y, family="cox", alpha=1) 
# Another attempt 
l.fit <- glmnet(x, Surv(sdata$tstart, sdata$tstop, sdata$death), family="cox", alpha=1) 

Answer 1

With a reasonable ratio of events to potential predictors, as you have here, a better approach than LASSO might be standard Cox regression. In outline: use knowledge of the subject matter to choose the predictors to consider, develop the model (with appropriate attention to linearity, proportional hazards assumptions, and so forth), use backward stepwise regression if you want a more parsimonious model, use bootstrapping to validate the model (examining the degree of optimism/overfitting), and simplify if necessary to make the model practically useful.

Chapter 4 of Frank Harrell's Regression Modeling Strategies provides a thorough description of how to proceed with regression modeling in such cases, and Chapter 21 illustrates the approach in practice with Cox modeling. Although stepwise selection, particularly unpenalized forward selection, has dangers, Harrell recognizes that model simplification is sometimes useful. Step 13 out of 20 in his approach (Section 4.12.1) is:

Do limited backwards step-down variable selection if parsimony is more important than accuracy. The cost of doing any aggressive variable selection is that the variable selection algorithm must also be included in a resampling procedure to properly validate the model or to compute confidence limits and the like.

Further steps deal with and correct for the potential overfitting. Finally, step 20 is: "Develop simplifications to the full model by approximating it to any desired degrees of accuracy (Section 5.5)." This process allows development of a practically useful model that minimizes the danger of overfitting. Harrell's rms package provides the necessary software tools.

Note, however, that time-dependent covariates will pose a problem with getting calibration curves for the model. (The validate function, which provides useful information on optimism, does seem to work with time-dependent covariates.) As Harrell said in a comment on a question* several years ago:

I should make it clear in the documentation but calibrate doesn't understand time-dependent covariates. As Therneau has stated frequently, estimation of survival probabilities in the presence of time-dependent covariates is not a simple thing to conceptualize.

That might be part of the reason why many standard tools that involve survival predictions don't handle time-dependent covariates. For example, if you have information on a covariate for a patient at 10 years, you already know that she survived for 10 years. What exactly are you predicting about that patient? For predictions, the baseline hazard needs to be adjusted individually for each patient to take into account the time-varying covariate values. The survfit function in the R survival package can do this, but many other survival-analysis functions and packages don't even try to handle the counting-process data structure. This raises the question whether a different approach like multi-state modeling might be more appropriate to your study.

If you are still worried about overfitting and you wish to use penalized methods (LASSO or ridge regression), the standard R coxph function allows for penalized maximum likelihood estimation. A ridge option (penalizing the sum of squares of coefficients) is provided, and it is possible to provide user-defined penalty functions so that you could write your own LASSO penalty (on the sum of absolute values of the coefficients) if you can't find a pre-written function that can do LASSO on data with time-dependent covariates.

Instead of cross-validation to choose penalty values, as is typically used with cv.glmnet, Harrell notes in Section 9.10 that optimizing a modified AIC performs well. In particular "for the lasso, with a fixed penalty parameter $\lambda$, the number of nonzero coefficients $k_{\lambda}$ is an unbiased estimate of the degrees of freedom" (Statistical Learning with Sparsity, page 18), simplifying calculation of the modified AIC.

Penalized methods, however, require some care in standardizing predictors so that penalization is applied fairly and measurement scales don't matter. There's a particular problem with trying to standardize categorical predictors, as the relative prevalence of categories and (for predictors with 3 or more levels) the choice of reference category will affect the standardization. This is discussed extensively in Section 9.10 of Harrell's book. Furthermore, LASSO will tend to choose one from among several correlated predictors, which isn't typically a problem with prediction but can seem somewhat arbitrary as choices can vary among different data samples.

If you still wish to use LASSO with cross-validation to choose penalties, there might be a way to proceed even with time-dependent covariates. The way that cv.glmnet judges Cox model fit is with the deviance (minus twice the log partial likelihood). The program uses a trick to handle the possibility that the number of events per fold is small, explained in section 3.5.1 of Statistical Learning with Sparsity. With the large scale of your data set, however, you could be able to get away with simply extracting the coefficients developed on the subset of cases that omits each fold, using those coefficients to calculate the deviance of the corresponding omitted cases, and summing over all folds.

---

*since deleted and potentially invisible, depending on reputation on this site