Por que a regressão glmnet ridge está me dando uma resposta diferente do cálculo manual?

28

Estou usando o glmnet para calcular estimativas de regressão de crista. Eu obtive alguns resultados que me fizeram suspeitar que o glmnet está realmente fazendo o que eu acho que ele faz. Para verificar isso, escrevi um script R simples, onde comparo o resultado da regressão de cume feita por resolve e a do glmnet, a diferença é 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

A norma da diferença é geralmente em torno de 20, o que não pode ser devido a algoritmos numericamente diferentes, devo estar fazendo algo errado. Quais são as configurações que tenho que definir glmnetpara obter o mesmo resultado que o cume?

John
fonte
1
Você já viu essa pergunta ?
Cdeterman
1
Sim, mas ainda não recebo o mesmo resultado usando a normalização.
John John
Você poderia postar seu código então?
shadowtalker
Acabei de ter o mesmo problema! a = data.frame (a = jitter (1:10), b = jitter (1:10), c = jitter (1:10), d = jitter (1:10), e = jitter (1:10) , f = instabilidade (1:10), g = amostra (instabilidade (1:10)), y = seq (10.100,10)); coef (lm.ridge (y ~ a + b + c + d + e + f + g, a, lambda = 2,57)); coef (glmnet (as.matrix (a [, 1: 7]), a $ y, família = "gaussiana", alfa = 0, lambda = 2,57 / 10)) Os resultados diferem bastante e se tornam muito mais semelhantes quando Eu uso lambdas muito mais altas para o glmnet.
A11msp
Intrigante. Os coeficientes parecem diferir aproximadamente pelo fator de 10.
Tomka

Respostas:

27

A diferença que você está observando se deve à divisão adicional pelo número de observações, N, que o GLMNET usa em sua função objetivo e padronização implícita de Y por seu desvio padrão da amostra, como mostrado abaixo.

12NysyXβ22+λβ22/2

onde usamos no lugar de 1 / ( n - 1 ) para s y , s y = i ( y i - ˉ y ) 21/n1/(n1)sy

sy=i(yiy¯)2n

Ao diferenciar em relação ao beta, definir a equação para zero,

XTXβXTysy+Nλβ=0

E resolvendo para beta, obtemos a estimativa,

β~GLMNET=(XTX+NλIp)1XTysy

Para recuperar as estimativas (e as respectivas penalidades correspondentes) na métrica original de Y, o GLMNET multiplica as estimativas e as lambdas por e retorna esses resultados ao usuário,sy

Xunstd. =s

β^GLMNET=syβ~GLMNET=(XTX+NλIp)1XTy
λunstd.=syλ

Compare esta solução com a derivação padrão da regressão de crista.

β^=(XTX+λIp)1XTy

Observe que é dimensionado por um fator extra de N. Além disso, quando usamos a função ou , a penalidade será implicitamente dimensionada em 1 / s y . Ou seja, quando usamos essas funções para obter as estimativas do coeficiente para alguns λ , estamos efetivamente obtendo estimativas para λ = λ / s y .λpredict()coef()1/syλλ=λ/sy

Com base nestas observações, a pena utilizados em GLMNET necessita de ser dimensionado por um factor de .sy/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

Os resultados generalizam para a inclusão de uma interceptação e variáveis ​​X padronizadas. Modificamos uma matriz X padronizada para incluir uma coluna de unidades e a matriz diagonal para ter uma entrada zero adicional na posição [1,1] (isto é, não penaliza a interceptação). Em seguida, você pode padronizar as estimativas pelos respectivos desvios padrão da amostra (verifique novamente se você está usando 1 / n ao calcular o desvio padrão).

β^j=βj~sxj

β^0=β0~x¯Tβ^
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

Adicionado código para mostrar X padronizado sem interceptação:

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
skijunkie
fonte
3
+6. Bem-vindo ao CV e obrigado por responder a essa pergunta antiga de maneira tão clara.
Ameba diz Reinstate Monica
1
ββ~
Percebo também que, na segunda parte, em que você disse "Os resultados generalizam para a inclusão de uma interceptação e variáveis ​​X padronizadas"; Para esta parte, se você excluir a interceptação, seguindo os mesmos cálculos, os resultados do glmnet se tornarão diferentes do cálculo manual.
user1769197
Correto, atualizei a solução com a matriz de identidade no lugar de βcomo necessário. Eu verifiquei a solução para o X padronizado sem interceptação e ainda obtenho resultados idênticos (consulte o código adicional acima).
skijunkie
3

De acordo com https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , quando a família estiver gaussian, glmnet()deve minimizar

(1)12nEu=1n(yEu-β0 0-xEuTβ)2+λj=1p(α|βj|+(1-α)βj2/2).

Ao usar glmnet(x, y, alpha=1)para ajustar o laço com as colunas emx padronizada, a solução para a penalidade relatada λ é a solução para minimizar

12nEu=1n(yEu-β0 0-xEuTβ)2+λj=1p|βj|.
No entanto, pelo menos em glmnet_2.0-13, ao usar glmnet(x, y, alpha=0)para ajustar a regressão de crista, a solução para uma penalidade relatadaλ é a solução para minimizar
12nEu=1n(yEu-β0 0-xEuTβ)2+λ12syj=1pβj2.
Onde sy é o desvio padrão de y. Aqui, a penalidade deveria ter sido relatada comoλ/sy.

O que pode acontecer é que a função primeiro padronize y para y0 0 e depois minimiza

2)12nEu=1n(y0 0Eu-xEuTγ)2+ηj=1p(α|γj|+(1-α)γj2/2),
que efetivamente é minimizar
12nsy2Eu=1n(yEu-β0 0-xEuTβ)2+ηαsyj=1p|βj|+η1-α2sy2j=1pβj2,
ou equivalente, para minimizar
12nEu=1n(yEu-β0 0-xEuTβ)2+ηsyαj=1p|βj|+η(1-α)j=1pβj2/2)

Para o laço (α=1), dimensionamento η de volta para relatar a penalidade como ηsyfaz sentido. Então para todosα, ηsy deve ser relatado como uma penalidade para manter a continuidade dos resultados α. Provavelmente, essa é a causa do problema acima. Isso se deve em parte ao uso de (2) para resolver (1). Apenas quandoα=0 0 ou α=1 existe alguma equivalência entre os problemas (1) e (2) (ou seja, uma correspondência entre os λ em (1) e o ηem 2)). Para qualquer outroα(0 0,1), os problemas (1) e (2) são dois problemas de otimização diferentes e não há correspondência individual entre os λ em (1) e o η em 2).

Chun Li
fonte
1
Não vejo onde sua resposta difere da anterior. Você poderia explicar, por favor?
Firebug 23/02
1
@Firebug Eu queria esclarecer por que a função relata o lambda dessa maneira, que parece não natural quando vista apenas da perspectiva da regressão de crista, mas faz sentido (ou tem que ser dessa maneira) quando vista da perspectiva de todo o espectro incluindo o cume e o laço.
Chun Li