2
$\begingroup$

I am trying to compare traditional Cox model and LASSO Cox in data with a counting process structure (see below for the data). I fitted a LASSO Cox model with lambda = 0, which in theory should lead to the same coefficients as the traditional Cox but doesn’t in practice.

Q1: why coxph() and glmnet() produce different coeffieicnets?

I also notice that glmnet() reports a warning message saying that cox.fit algorithm did not converge.

Q2: Why coxph() fits the model without any convergence issue but glmnet() has?

Really appreciate the help.

# load package
library(tidyverse)
library(survival)
library(glmnet)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
#> Loaded glmnet 4.1-3


# import data
data_death <-
  structure(
    list(
      person_id = c(1L, 2L, 2L, 3L, 3L, 4L, 4L, 4L,
                    5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L),
      age = c(20, 21, 21, 19, 19,
              22, 22, 22, 20, 20, 20, 20, 24, 24, 24, 24),
      female = c(0L, 1L,
                 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L),
      time0 = c(0L,
                0L, 1L, 0L, 1L, 0L, 1L, 7L, 0L, 1L, 7L, 10L, 0L, 1L, 7L, 10L),
      time1 = c(1L, 1L, 4L, 1L, 7L, 1L, 7L, 10L, 1L, 7L, 10L, 12L,
                1L, 7L, 10L, 13L),
      death = c(1L, 0L, 0L, 0L, 1L, 0L, 0L,
                1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L)
    ),
    row.names = c(NA,-16L),
    class = c("tbl_df", "tbl", "data.frame")
  )


data_death
#> # A tibble: 16 x 6
#>    person_id   age female time0 time1 death
#>        <int> <dbl>  <int> <int> <int> <int>
#>  1         1    20      0     0     1     1
#>  2         2    21      1     0     1     0
#>  3         2    21      1     1     4     0
#>  4         3    19      0     0     1     0
#>  5         3    19      0     1     7     1
#>  6         4    22      1     0     1     0
#>  7         4    22      1     1     7     0
#>  8         4    22      1     7    10     1
#>  9         5    20      0     0     1     0
#> 10         5    20      0     1     7     0
#> 11         5    20      0     7    10     0
#> 12         5    20      0    10    12     0
#> 13         6    24      1     0     1     0
#> 14         6    24      1     1     7     0
#> 15         6    24      1     7    10     0
#> 16         6    24      1    10    13     1



# fit traditional cox model
model_cox <- 
  coxph(Surv(time = time0, 
             time2 = time1, 
             event = death,
             type = "counting") ~ female + age, 
        data = data_death)



# fit lasso cox model (with a penalty of 0)
model_lasso <- 
  glmnet(x = data_death %>% select(age, female) %>% as.matrix(),
         y = Surv(time = data_death$time0,
                  time2 = data_death$time1,
                  event = data_death$death,
                  type = "counting"),
         family = "cox",
         lambda = 0) 
#> Warning: cox.fit: algorithm did not converge


# compare model coefficient
model_cox
#> Call:
#> coxph(formula = Surv(time = time0, time2 = time1, event = death, 
#>     type = "counting") ~ female + age, data = data_death)
#> 
#>           coef exp(coef) se(coef)      z     p
#> female  1.5446    4.6860   2.7717  0.557 0.577
#> age    -0.9453    0.3886   1.0637 -0.889 0.374
#> 
#> Likelihood ratio test=1.65  on 2 df, p=0.4378
#> n= 16, number of events= 4
model_lasso$beta
#> 2 x 1 sparse Matrix of class "dgCMatrix"
#>                  s0
#> age    -0.009869489
#> female -0.872824209

Created on 2021-11-30 by the reprex package (v2.0.0)

$\endgroup$
2
  • $\begingroup$ Does the reproducible example in the section of the Cox glmnet vignette on "Cox models for start-stop data" work for you? $\endgroup$
    – EdM
    Dec 6, 2021 at 17:10
  • $\begingroup$ Yes. The example there works for me perfectly. Fitting glmnet on the start-stop data there doesn't produce any warning about model convergence. Also, the coefficients of coxph() and glmnet() are almost the same. $\endgroup$
    – zeming
    Dec 8, 2021 at 16:52

1 Answer 1

2
$\begingroup$

This seems to have something to do with the combination of the peculiarities of this data set and the different ways that coxph() and glmnet() fit models. You have substantial collinearity. All males have ages of 19 or 20, while all females are 21 or older. I can't say why that might pose a problem for glmnet() but not for coxph(), but presumably the open-source code could provide an explanation.

If you work with standardized values for both age and female you can get glmnet() to fit the model without warning:

data_death[,"stdAge"] <- (data_death$age - mean(data_death$age))/sd(data_death$age)
data_death[,"stdFemale"] <- (data_death$female - mean(data_death$female))/sd(data_death$female)

> model_lasso_stdAgeFemale <- 
  glmnet(x = as.matrix(data_death[,c("stdAge","stdFemale")]),
         y = Surv(time = data_death$time0,
                  time2 = data_death$time1,
                  event = data_death$death,
                  type = "counting"),
         family = "cox",
         lambda = 0)
 
model_lasso_stdAgeFemale$beta
# 2 x 1 sparse Matrix of class "dgCMatrix"
#                   s0
# stdAge    -1.7204162
# stdFemale  0.7906132

and back-correcting for the standardization:

0.7906132/sd(data_death$female)
# [1] 1.543119
-1.7204162/sd(data_death$age)
# [1] -0.9446754

you get close to your coxph coefficients. The glmnet() function is supposed to standardize predictor values by default; can't say what's going on here. The coxph() code does standardize internally by default. From the manual page:

The routine internally scales and centers data to avoid overflow in the argument to the exponential function. These actions do not change the result, but lead to more numerical stability.

$\endgroup$
3
  • $\begingroup$ Apologies for the late reply. Thank you very much for your answer. In the data for my own research (which I can't share here), I did standardised all predictors there but still have convergence issue. I will have a look at that data to make sure I did complete the standardisation. $\endgroup$
    – zeming
    Dec 15, 2021 at 10:35
  • $\begingroup$ @zeming I got the glmnet() error with your data if I only standardized age and not the categorical female; glmnet() and coxph() clearly differ in some fitting details that I don't understand. You'd have to look at the code. Having categorical predictor levels (or combinations in interactions) without events can also lead to lack of convergence. That wasn't the case in the example data you showed here, but it happens in practice and it's a bigger problem when you use cross validation to choose the penalty. Look at the data for that issue and for severe multicollinearity. $\endgroup$
    – EdM
    Dec 15, 2021 at 13:58
  • $\begingroup$ Thank you for your answer. I will have a check on my data. Have a good day. $\endgroup$
    – zeming
    Dec 16, 2021 at 9:20

Your Answer

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

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