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!



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)

# 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) 
  • $\begingroup$ Note that questions about software packages per se are off topic here. If you could edit the question to say more about how you wish to use the penalization provided by glmnet and the scale and goal of your actual problem (number of cases, number of events, predictors, goals of modeling, intended use of model, etc) then there might be a statistical issue that is on topic (e.g., the criteria that glmnet uses for cross-validation, the difficulties in making predictions with time-varying covariates) and that could point you toward a solution for your problem. $\endgroup$
    – EdM
    May 10, 2019 at 18:07
  • $\begingroup$ Very useful addition to the question. One further question before I try an answer: just how many and what type of predictors are you considering? With, say, 2000 events it might be possible to include on the order of 100-200 predictors in a standard Cox model. That should be enough to handle most clinical variables. Or are you trying to incorporate gene-expression data or some other high-dimensional data, such that the number of predictors is greater than the number of events? $\endgroup$
    – EdM
    May 13, 2019 at 14:29

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

  • $\begingroup$ Thank you for the detailed response. Your advice is very sensible and lines up with the basic approach we had adopted. In this particular phase of the project, we had wanted to try a penalized approach and examine the potential for overfitting. I had not thought to return to coxph and write my own penalty function - thank you for that idea. $\endgroup$ May 14, 2019 at 2:05
  • $\begingroup$ @RobertBalshaw if you do use penalization like LASSO or ridge, think very carefully about how to handle the categorical predictors. Whether or how to standardize them before penalization is problematic. Again, see Section 9.10 of Harrell's book. This page provides examples of how category prevalence can affect the relative weights applied to predictors in penalization. $\endgroup$
    – EdM
    May 14, 2019 at 14:35
  • $\begingroup$ Again - thank you. The notion of standardization is very much on the table right now. My own feeling is that the 0-1 for binary predictors is already appropriately "scaled" and we are trying to work with that as a starting point. Multiple category predictors are the real problem - and we may just avoid them or work out a way to express them using binary predictors. $\endgroup$ May 14, 2019 at 16:43
  • $\begingroup$ @RobertBalshaw see Verweij and van Houwelingen, Statistics in Medicine 13:2427, 1994. In addition to being a seminal paper on Cox regression cross-validation, it shows reasonable ways to handle multi-category and ordinal predictors. Harrell also discusses this in his Section 9.10. So don't give up on multiple-category predictors. Just be smart with them. $\endgroup$
    – EdM
    May 14, 2019 at 17:05
  • $\begingroup$ A simple upvote is insufficient to express what an amazing answer this is. $\endgroup$
    – matmat
    Jan 5, 2022 at 6:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.