Mengapa regresi ridge glmnet memberi saya jawaban berbeda dari perhitungan manual?

28

Saya menggunakan glmnet untuk menghitung estimasi regresi ridge. Saya mendapatkan beberapa hasil yang membuat saya curiga bahwa glmnet benar-benar melakukan apa yang menurut saya benar. Untuk memeriksa ini saya menulis sebuah skrip R sederhana di mana saya membandingkan hasil regresi ridge yang dilakukan oleh resol dan yang di glmnet, perbedaannya signifikan:

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

Norma perbedaan biasanya sekitar 20 yang tidak dapat disebabkan oleh algoritma yang berbeda secara numerik, saya pasti melakukan sesuatu yang salah. Pengaturan apa yang harus saya atur glmnetuntuk mendapatkan hasil yang sama seperti dengan ridge?

John
sumber
1
Pernahkah Anda melihat pertanyaan ini ?
cdeterman
1
Ya, tapi saya masih belum mendapatkan hasil yang sama menggunakan normalisasi.
John
Bisakah Anda memposting kode Anda?
shadowtalker
Saya baru saja mengalami masalah yang sama! a = data.frame (a = jitter (1:10), b = jitter (1:10), c = jitter (1:10), d = jitter (1:10), e = jitter (1:10) , f = jitter (1:10), g = sampel (jitter (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, family = "gaussian", alpha = 0, lambda = 2.57 / 10)) Hasilnya sedikit berbeda dan menjadi jauh lebih mirip ketika Saya menggunakan lambdas yang jauh lebih tinggi untuk glmnet.
a11msp
Menarik. Koefisien tampaknya berbeda secara kasar dengan faktor 10.
tomka

Jawaban:

27

Perbedaan yang Anda amati adalah karena pembagian tambahan dengan jumlah pengamatan, N, yang digunakan GLMNET dalam fungsi tujuan mereka dan standardisasi implisit Y oleh sampel standar deviasi seperti yang ditunjukkan di bawah ini.

12NysyXβ22+λβ22/2

di mana kita menggunakan sebagai pengganti untuk , 1 / ( n - 1 ) s y s y = i ( y i - ˉ y ) 21/n1/(n1)sy

sy=i(yiy¯)2n

Dengan membedakan sehubungan dengan beta, mengatur persamaan ke nol,

XTXβXTysy+Nλβ=0

Dan pemecahan untuk beta, kami memperoleh estimasi,

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

Untuk memulihkan taksiran (dan hukumannya yang sesuai) pada metrik Y asli, GLMNET mengalikan taksiran dan dengan dan mengembalikan hasil ini kepada pengguna,sy

λunstd. =syλ

β^GL.M.NET=syβ~GL.M.NET=(XTX+Nλsayahal)-1XTy
λkamunstd.=syλ

Bandingkan solusi ini dengan derivasi standar regresi ridge.

β^=(XTX+λsayahal)-1XTy

Perhatikan bahwa diskalakan oleh faktor tambahan N. Selain itu, ketika kita menggunakan fungsi atau , penalti akan secara implisit diskalakan oleh . Dengan kata lain, ketika kita menggunakan fungsi-fungsi ini untuk mendapatkan estimasi koefisien untuk beberapa , kita secara efektif mendapatkan estimasi untuk .1 / s y λ λ = λ / s yλpredict()coef()1/syλλ=λ/sy

Berdasarkan pengamatan ini, hukuman yang digunakan dalam GLMNET perlu ditingkatkan dengan faktor .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

Hasil generalisasi dengan dimasukkannya variabel X intersep dan standar. Kami memodifikasi matriks X standar untuk memasukkan kolom yang dan matriks diagonal untuk memiliki entri nol tambahan di posisi [1,1] (yaitu tidak menghukum intersep). Anda kemudian dapat menghapus standar estimasi dengan masing-masing standar deviasi sampel (sekali lagi memastikan Anda menggunakan 1 / n saat menghitung standar deviasi).

β^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

Kode ditambahkan untuk menunjukkan X standar tanpa intersep:

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
sumber
3
+6. Selamat datang di CV dan terima kasih telah menjawab pertanyaan lama ini dengan sangat jelas.
Amoeba mengatakan Reinstate Monica
1
Seharusnya matriks identitas bukan dalam solusi , benar? ˜ βββ~
user1769197
Saya juga memperhatikan bahwa untuk bagian kedua di mana Anda mengatakan "Hasilnya menggeneralisasi untuk dimasukkannya variabel X yang mencegat dan standar"; untuk bagian ini, jika Anda mengecualikan intersep, kemudian mengikuti perhitungan yang sama, hasil glmnet menjadi berbeda dari perhitungan manual.
user1769197
Benar, saya telah memperbarui solusi dengan matriks identitas sebagai pengganti sesuai kebutuhan. Saya memeriksa solusi untuk X standar tanpa intersep dan masih mendapatkan hasil yang identik (lihat kode tambahan di atas). β
skijunkie
3

Menurut https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , ketika keluarga itu gaussian, glmnet()harus meminimalkan

(1)12nsaya=1n(ysaya-β0-xsayaTβ)2+λj=1hal(α|βj|+(1-α)βj2/2).

Ketika menggunakan glmnet(x, y, alpha=1)agar sesuai dengan laso dengan kolom dalam standar, solusi untuk penalti yang dilaporkan adalah solusi untuk meminimalkan Namun, setidaknya dalam , ketika menggunakan agar sesuai dengan regresi ridge, solusi untuk penalti yang dilaporkan adalah solusi untuk meminimalkan di mana adalah standar deviasi . Di sini, hukuman seharusnya dilaporkan sebagai .xλ

12nsaya=1n(ysaya-β0-xsayaTβ)2+λj=1hal|βj|.
glmnet_2.0-13glmnet(x, y, alpha=0)λ
12nsaya=1n(ysaya-β0-xsayaTβ)2+λ12syj=1halβj2.
syyλ/sy

Apa yang mungkin terjadi adalah bahwa fungsi pertama-tama menstandarkan ke dan kemudian meminimalkan yang secara efektif adalah untuk meminimalkan atau yang setara, untuk meminimalkan yy0

(2)12nsaya=1n(y0saya-xsayaTγ)2+ηj=1hal(α|γj|+(1-α)γj2/2),
12nsy2saya=1n(ysaya-β0-xsayaTβ)2+ηαsyj=1hal|βj|+η1-α2sy2j=1halβj2,
12nsaya=1n(ysaya-β0-xsayaTβ)2+ηsyαj=1hal|βj|+η(1-α)j=1halβj2/2.

Untuk laso ( ), skala kembali untuk melaporkan penalti karena masuk akal. Maka untuk semua , harus dilaporkan sebagai penalti untuk menjaga kesinambungan hasil di seluruh . Ini mungkin adalah penyebab masalah di atas. Ini sebagian karena menggunakan (2) untuk menyelesaikan (1). Hanya ketika atau ada beberapa kesetaraan antara masalah (1) dan (2) (yaitu, korespondensi antara di (1) dan in (2)). Untukα=1ηηsyαηsyαα=0α=1ληα(0,1), masalah (1) dan (2) adalah dua masalah pengoptimalan yang berbeda, dan tidak ada korespondensi satu-ke-satu antara di (1) dan di (2).λη

Chun Li
sumber
1
Saya tidak bisa melihat di mana jawaban Anda berbeda dari yang sebelumnya. Bisakah Anda jelaskan?
Firebug
1
@ Firebug Saya ingin menjelaskan mengapa fungsi melaporkan lambda dengan cara ini, yang tampak tidak wajar jika dilihat hanya dari sudut pandang regresi ridge, tetapi masuk akal (atau harus seperti ini) jika dilihat dari perspektif seluruh spektrum termasuk punggungan dan laso.
Chun Li