Resultados diferentes para glmnet cuando se estandariza = FALSE

Question

Hice esta pregunta en stackoverflow, pero en este momento no sé exactamente a qué lugar pertenece porque es una pregunta relacionada con el proceso de estandarización de glmnet y el lazo.

Estoy corriendo glmnet y tratando de ver la diferencia al usar standardize = FALSE con variables pre-estandarizadas.

library(glmnet)
data(QuickStartExample)
### Standardized Way
x_standardized <- scale(x, center = TRUE, scale = TRUE)
y_standardized <- scale(y, center = TRUE, scale = TRUE)
cv_standardized <- cv.glmnet(x_standardized, y_standardized,
                             intercept = FALSE,
                             standardize = FALSE, standardize.response = FALSE)
destandardized_coef <- coef(cv_standardized)[-1] * sd(y) / apply(x, 2, sd)
destandardized_coef
mean(y) - sum(destandardized_coef * colMeans(x))

### Let glmnet Stanardize
cv_normal <- cv.glmnet(x, y)
coef(cv_normal, cv_normal$lambda.min) %>% as.numeric()

Inicialmente, estoy estandarizando los datos yo mismo y transformando los coeficientes hacia atrás. Esperaría ver los mismos resultados, pero por alguna razón estoy obteniendo coeficientes ligeramente diferentes.

Mi pregunta es, ¿cómo puedo extraer los mismos resultados, y por qué los coeficientes son actualmente diferentes de esta manera?

EDIT: Aquí está mi código actualizado basado en la respuesta de @user2974951, parece que glmnet estandariza de forma diferente a la scale función en R . Quiero hacer notar que estoy ajustando el modelo inicial sin el intercepto ya que luego lo estoy calculando a mano. En última instancia, esto no debería importar ya que proporciono un modelo estandarizado x y y el intercepto debe ser cero.

set.seed(123)
library(glmnet)
library(tidyverse)
data(QuickStartExample)
# standardize data
n <- length(y)
y_mean <- mean(y)
y_centered <- y - mean(y)
y_sd <- sqrt(sum((y - y_mean) ^ 2) / n)
ys <- y_centered / y_sd
X_centered <- apply(x, 2, function(x) x - mean(x))
Xs <- apply(X_centered, 2, function(x) x / sqrt(sum(x ^ 2) / n))
ys <- y_centered / sqrt(sum(y_centered ^ 2) / n)

set.seed(123)
cv_standardized <- cv.glmnet(Xs, ys,
                             intercept = FALSE,
                             standardize = FALSE,
                             standardize.response = FALSE)
destandardized_coef <- coef(cv_standardized)[-1] * sd(y) / apply(x, 2, sd)
personal_standardize <- c(mean(y) - sum(destandardized_coef * colMeans(x)),
                          destandardized_coef)

set.seed(123)
cv_normal <- cv.glmnet(x, y, intercept = TRUE)
glmnet_standardize <- coef(cv_normal, cv_normal$lambda.min)

Resultados:

cbind(personal_standardize, glmnet_standardize)
21 x 2 sparse Matrix of class "dgCMatrix"
            personal_standardize  glmnet_standardize
(Intercept)           0.15473991  0.145832036
V1                    1.29441444  1.340981414
V2                    .           .          
V3                    0.62890756  0.708347139
V4                    .           .          
V5                   -0.77785401 -0.848087765
V6                    0.48387954  0.554823781
V7                    .           0.038519738
V8                    0.28419264  0.347221319
V9                    .           .          
V10                   .           0.010034050
V11                   0.09130386  0.186365264
V12                   .           .          
V13                   .           .          
V14                  -1.03241106 -1.086006902
V15                   .          -0.004444318
V16                   .           .          
V17                   .           .          
V18                   .           .          
V19                   .           .          
V20                  -0.96190123 -1.069942845

Gracias de antemano.

Answer 1

Tuve un pequeño problema con una errata pero pude solucionarlo. Basado en la respuesta de @user2974951 y la página de Github para glmnet Parece que están usando 1/n estandarización en lugar de 1/n-1 . He editado el código y he conseguido los mismos resultados.

set.seed(123)
library(glmnet)
library(tidyverse)
data(QuickStartExample)
# standardize data
n <- length(y)
y_centered <- scale(y, center = TRUE, scale = FALSE)
y_sd <- sqrt(sum((y - mean(y)) ^ 2) / n)
ys <- y_centered / y_sd
x_centered <- scale(x, center = TRUE, scale = FALSE)
x_sd <- apply(x, 2, function(x) sqrt(sum((x - mean(x)) ^ 2) / n))
xs <- matrix(ncol = ncol(x), nrow = nrow(x))
for (i in seq_len(ncol(x))) {
  xs[, i] <- x_centered[, i] / x_sd[i]
}
### glmnet using standardized data
set.seed(123)
cv_standardized <- cv.glmnet(xs, ys,
                             intercept = FALSE,
                             standardize = FALSE,
                             standardize.response = FALSE)
destandardized_coef <- coef(cv_standardized, cv_standardized$lambda.min)[-1] * y_sd /
  apply(x, 2, function(x) sqrt(sum((x - mean(x)) ^ 2) / n))
personal_standardize <- c(mean(y) - sum(destandardized_coef * colMeans(x)),
                          destandardized_coef)

### Letting glmnet standardize data
set.seed(123)
cv_normal <- cv.glmnet(x, y, intercept = TRUE)
coef(cv_normal, cv_normal$lambda.min) %>% as.numeric()
glmnet_standardize <- coef(cv_normal, cv_normal$lambda.min)

Resultados:

cbind(personal_standardize, glmnet_standardize)
21 x 2 sparse Matrix of class "dgCMatrix"
            personal_standardize            1
(Intercept)          0.145832036  0.145832036
V1                   1.340981414  1.340981414
V2                   .            .          
V3                   0.708347139  0.708347139
V4                   .            .          
V5                  -0.848087765 -0.848087765
V6                   0.554823781  0.554823781
V7                   0.038519738  0.038519738
V8                   0.347221319  0.347221319
V9                   .            .          
V10                  0.010034050  0.010034050
V11                  0.186365264  0.186365264
V12                  .            .          
V13                  .            .          
V14                 -1.086006902 -1.086006902
V15                 -0.004444318 -0.004444318
V16                  .            .          
V17                  .            .          
V18                  .            .          
V19                  .            .          
V20                 -1.069942845 -1.069942845

Incluso puede obtener el desestandarizado lambda del método de estandarización personal por:

> cv_standardized$lambda.min * y_sd
[1] 0.06284188
> cv_normal$lambda.min
[1] 0.06284188

También, set.seed(123) me ayudó a encontrar la solución a este problema. Gracias por la ayuda.