Convergence issue when fitting LASSO Cox using glmnet() in R

Question

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)}

Answer 1

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.