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 glmnet
untuk mendapatkan hasil yang sama seperti dengan ridge?
r
ridge-regression
glmnet
John
sumber
sumber
Jawaban:
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.
di mana kita menggunakan sebagai pengganti untuk , 1 / ( n - 1 ) s y s y = ∑ i ( y i - ˉ y ) 21/n 1/(n−1) sy
Dengan membedakan sehubungan dengan beta, mengatur persamaan ke nol,
Dan pemecahan untuk beta, kami memperoleh estimasi,
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λ
Bandingkan solusi ini dengan derivasi standar regresi ridge.
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λ 1 / dtky λ∗ λ = λ∗/ sy
predict()
coef()
Berdasarkan pengamatan ini, hukuman yang digunakan dalam GLMNET perlu ditingkatkan dengan faktor .sy/ N
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).
Kode ditambahkan untuk menunjukkan X standar tanpa intersep:
sumber
Menurut https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , ketika keluarga itu
gaussian
,glmnet()
harus meminimalkanKetika menggunakanx λ
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 .glmnet_2.0-13
glmnet(x, y, alpha=0)
Apa yang mungkin terjadi adalah bahwa fungsi pertama-tama menstandarkan ke dan kemudian meminimalkan yang secara efektif adalah untuk meminimalkan atau yang setara, untuk meminimalkany y0
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).λ η
sumber