¿Por qué la regresión de crestas de glmnet me da una respuesta diferente a la del cálculo manual?

Question

Estoy utilizando glmnet para calcular las estimaciones de la regresión ridge. Obtuve algunos resultados que me hicieron sospechar que glmnet realmente hace lo que yo creo que hace. Para comprobarlo escribí un sencillo R script en el que comparo el resultado de la regresión ridge realizada por solve y el de glmnet, la diferencia es significativa:

n    <- 1000
p.   <-  100
X.   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

beta1 <- solve(t(X)%*%X+5*diag(p),t(X)%*%Y)
beta2 <- glmnet(X,Y, alpha=0, lambda=10, intercept=FALSE, standardize=FALSE, 
                family="gaussian")$beta@x
beta1-beta2

La norma de la diferencia suele estar en torno a 20, lo que no puede deberse a algoritmos numéricamente diferentes, debo estar haciendo algo mal. ¿Cuáles son los ajustes que tengo que establecer en glmnet para obtener el mismo resultado que con la cresta?

Answer 1

La diferencia que usted observa se debe a la división adicional por el número de observaciones, N, que GLMNET utiliza en su función objetivo y a la estandarización implícita de Y por su desviación estándar muestral, como se muestra a continuación.

$$ \frac{1}{2N}\left\|\frac{y}{s_y}-X\beta\right\|^2_{2}+\lambda\|\beta\|^2_{2}/2 $$

donde utilizamos $1/n$ en lugar de $1/(n-1)$ para $s_y$ , $$ s_y=\frac{\sum_i(y_i-\bar{y})^2}{n} $$

Diferenciando con respecto a beta, poniendo la ecuación a cero,

$$ X^TX\beta-\frac{X^Ty}{s_y}+N\lambda\beta =0 $$

Y resolviendo para beta, obtenemos la estimación,

$$ \tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}\frac{X^Ty}{s_y} $$

Para recuperar las estimaciones (y sus correspondientes penalizaciones) sobre la métrica original de Y, GLMNET multiplica tanto las estimaciones como las lambdas por $s_y$ y devuelve estos resultados al usuario,

$$ \hat{\beta}_{GLMNET}=s_y\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}X^Ty $$ $$ \lambda_{unstd.}=s_y\lambda $$

Compare esta solución con la derivación estándar de la regresión de cresta.

$$ \hat{\beta}= (X^TX+\lambda I_p)^{-1}X^Ty $$

Observe que $\lambda$ se escala por un factor adicional de N. Además, cuando utilizamos el predict() o coef() la penalización se va a escalar implícitamente por $1/s_y$ . Es decir, cuando utilizamos estas funciones para obtener las estimaciones de los coeficientes de algunos $\lambda^*$ estamos obteniendo efectivamente estimaciones para $\lambda=\lambda^*/s_y$ .

Basándose en estas observaciones, la penalización utilizada en GLMNET necesita ser escalada por un factor de $s_y/N$ .

set.seed(123)

n    <- 1000
p   <-  100
X   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

beta1 <- solve(t(X)%*%X+10*diag(p),t(X)%*%(Y))[,1]

fit_glmnet <- glmnet(X,Y, alpha=0, standardize = F, intercept = FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

           [,1]        [,2]
[1,]  0.23793862  0.23793862
[2,]  1.81859695  1.81859695
[3,] -0.06000195 -0.06000195
[4,] -0.04958695 -0.04958695
[5,]  0.41870613  0.41870613
[6,]  1.30244151  1.30244151
[7,]  0.06566168  0.06566168
[8,]  0.44634038  0.44634038
[9,]  0.86477108  0.86477108
[10,] -2.47535340 -2.47535340

Los resultados se generalizan con la inclusión de un intercepto y variables X estandarizadas. Modificamos una matriz X estandarizada para incluir una columna de unos y la matriz diagonal para tener una entrada cero adicional en la posición [1,1] (es decir, no penalizamos el intercepto). A continuación, puede desestandarizar las estimaciones por sus respectivas desviaciones estándar de la muestra (de nuevo, asegúrese de que está utilizando 1/n al calcular la desviación estándar).

$$ \hat\beta_{j}=\frac{\tilde{\beta_j}}{s_{x_j}} $$

$$ \hat\beta_{0}=\tilde{\beta_0}-\bar{x}^T\hat{\beta} $$

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)
X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}
X_scaled_ones <- cbind(rep(1,n), X_scaled)

beta3 <- solve(t(X_scaled_ones)%*%X_scaled_ones+1000*diag(x = c(0, rep(1,p))),t(X_scaled_ones)%*%(Y))[,1]
beta3 <- c(beta3[1] - crossprod(mean_x,beta3[-1]/sd_x), beta3[-1]/sd_x)

fit_glmnet2 <- glmnet(X,Y, alpha=0, thresh = 1e-20)
beta4 <- as.vector(coef(fit_glmnet2, s = sd_y*1000/n, exact = TRUE))

cbind(beta3[1:10], beta4[1:10])
             [,1]        [,2]
 [1,]  0.24534485  0.24534485
 [2,]  0.17661130  0.17661130
 [3,]  0.86993230  0.86993230
 [4,] -0.12449217 -0.12449217
 [5,] -0.06410361 -0.06410361
 [6,]  0.17568987  0.17568987
 [7,]  0.59773230  0.59773230
 [8,]  0.06594704  0.06594704
 [9,]  0.22860655  0.22860655
[10,]  0.33254206  0.33254206

Se ha añadido código para mostrar la X estandarizada sin intercepción:

set.seed(123)

n <- 1000
p <-  100
X <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)

X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}

beta1 <- solve(t(X_scaled)%*%X_scaled+10*diag(p),t(X_scaled)%*%(Y))[,1]

fit_glmnet <- glmnet(X_scaled,Y, alpha=0, standardize = F, intercept = 
FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

             [,1]        [,2]
 [1,]  0.23560948  0.23560948
 [2,]  1.83469846  1.83469846
 [3,] -0.05827086 -0.05827086
 [4,] -0.04927314 -0.04927314
 [5,]  0.41871870  0.41871870
 [6,]  1.28969361  1.28969361
 [7,]  0.06552927  0.06552927
 [8,]  0.44576008  0.44576008
 [9,]  0.90156795  0.90156795
[10,] -2.43163420 -2.43163420

Answer 2

Según https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html cuando la familia es gaussian , glmnet() debe minimizar $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda\sum_{j=1}^p(\alpha|\beta_j| +(1-\alpha)\beta_j^2/2). \tag{1}$$

Al utilizar glmnet(x, y, alpha=1) para ajustar el lazo con las columnas en $x$ estandarizada, la solución para la penalización reportada $\lambda$ es la solución para minimizar $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \sum_{j=1}^p |\beta_j|.$$ Sin embargo, al menos en glmnet_2.0-13 cuando se utiliza glmnet(x, y, alpha=0) para ajustar la regresión de cresta, la solución para una penalización informada $\lambda$ es la solución para minimizar $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \frac{1}{2s_y} \sum_{j=1}^p \beta_j^2.$$ donde $s_y$ es la desviación estándar de $y$ . En este caso, la sanción debería haberse comunicado como $\lambda/s_y$ .

Lo que puede ocurrir es que la función primero normalice $y$ a $y_0$ y luego minimiza $$\frac{1}{2n} \sum_{i=1}^n (y_{0i}-x_i^T\gamma)^2 +\eta \sum_{j=1}^p(\alpha|\gamma_j| +(1-\alpha)\gamma_j^2/2), \tag{2}$$ que efectivamente es minimizar $$\frac{1}{2n s_y^2} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta \frac{\alpha}{s_y} \sum_{j=1}^p |\beta_j| +\eta \frac{1-\alpha}{2s_y^2} \sum_{j=1}^p \beta_j^2,$$ o, de forma equivalente, para minimizar $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta s_y \alpha \sum_{j=1}^p |\beta_j| +\eta (1-\alpha) \sum_{j=1}^p \beta_j^2/2.$$

Para el lazo ( $\alpha=1$ ), escalando $\eta$ para informar de la sanción como $\eta s_y$ tiene sentido. Entonces para todos $\alpha$ , $\eta s_y$ tiene que ser reportado como la pena para mantener la continuidad de los resultados a través de $\alpha$ . Esta es probablemente la causa del problema anterior. Esto se debe, en parte, a la utilización de (2) para resolver (1). Sólo cuando $\alpha=0$ o $\alpha=1$ existe cierta equivalencia entre los problemas (1) y (2) (es decir, una correspondencia entre los $\lambda$ en (1) y el $\eta$ en (2)). Para cualquier otro $\alpha\in(0,1)$ Los problemas (1) y (2) son dos problemas de optimización diferentes, y no existe una correspondencia unívoca entre los $\lambda$ en (1) y el $\eta$ en (2).