Variablity dalam hasil cv.glmnet

18

Saya menggunakan cv.glmnetuntuk menemukan prediktor. Pengaturan yang saya gunakan adalah sebagai berikut:

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,nfolds=cvfold)
bestlambda<-lassoResults$lambda.min

results<-predict(lassoResults,s=bestlambda,type="coefficients")

choicePred<-rownames(results)[which(results !=0)]

Untuk membuat yakin hasilnya direproduksi saya set.seed(1). Hasilnya sangat bervariasi. Saya menjalankan kode yang sama persis 100 untuk melihat bagaimana variabel hasilnya. Dalam menjalankan 98/100 memiliki satu prediktor tertentu selalu dipilih (kadang-kadang hanya pada itu sendiri); prediktor lain dipilih (co-efisien bukan nol) biasanya 50/100 kali.

Jadi dikatakan kepada saya bahwa setiap kali validasi silang berjalan itu mungkin akan memilih lambda terbaik yang berbeda, karena pengacakan awal dari lipatan penting. Orang lain telah melihat masalah ini ( hasil CV.glmnet ) tetapi tidak ada solusi yang disarankan.

Saya berpikir bahwa mungkin yang muncul 98/100 mungkin sangat berkorelasi dengan yang lain? Hasilnya memang stabil jika saya hanya menjalankan LOOCV ( fold-size=n ), tapi saya ingin tahu mengapa mereka sangat variabel ketika nfold<n .

pengguna4673
sumber
1
Untuk lebih jelasnya, apakah maksud Anda Anda set.seed(1)menjalankan cv.glmnet()100 kali? Itu bukan metodologi yang bagus untuk reproduksibilitas; lebih baik ke set.seed()kanan sebelum setiap lari, atau jaga agar lipatannya konstan di seluruh lintasan. Setiap panggilan Anda untuk cv.glmnet()menelepon sample()N kali. Jadi jika panjang data Anda pernah berubah, maka reprodubilitas akan berubah.
smci

Jawaban:

14

Intinya di sini adalah bahwa dalam cv.glmnetlipatan K ("bagian") diambil secara acak.

Dalam validasi silang K-folds, dataset dibagi dalam bagian , dan bagian K - 1 digunakan untuk memprediksi bagian K-th (ini dilakukan kali K , menggunakan bagian K yang berbeda setiap kali). Ini dilakukan untuk semua lambda, dan yang memberikan kesalahan validasi silang terkecil.KK1KKlambda.min

Inilah sebabnya ketika Anda menggunakan hasilnya tidak berubah: setiap grup terbuat dari satu, jadi tidak ada banyak pilihan untuk grup K.nfolds=nK

Dari cv.glmnet()manual referensi:

Perhatikan juga bahwa hasil cv.glmnet adalah acak, karena lipatannya dipilih secara acak. Pengguna dapat mengurangi keacakan ini dengan menjalankan cv.glmnet berkali-kali, dan rata-rata kurva kesalahan.

### cycle for doing 100 cross validations
### and take the average of the mean error curves
### initialize vector for final data.frame with Mean Standard Errors
MSEs <- NULL
for (i in 1:100){
                 cv <- cv.glmnet(y, x, alpha=alpha, nfolds=k)  
                 MSEs <- cbind(MSEs, cv$cvm)
             }
  rownames(MSEs) <- cv$lambda
  lambda.min <- as.numeric(names(which.min(rowMeans(MSEs))))

UMK adalah kerangka data yang berisi semua kesalahan untuk semua lambdas (untuk 100 berjalan), lambda.minadalah lambda Anda dengan kesalahan rata-rata minimum.

Alice
sumber
Hal yang paling saya khawatirkan adalah bahwa pemilihan n benar-benar tampaknya penting kadang-kadang. Haruskah saya memercayai hasil yang bisa sangat bervariasi? Atau haruskah saya menganggapnya samar meskipun saya menjalankannya berkali-kali?
user4673
1
Bergantung pada ukuran sampel Anda harus memilih n sehingga Anda memiliki setidaknya 10 pengamatan per kelompok. Jadi lebih baik untuk mengurangi n standar (= 10) jika Anda memiliki ukuran sampel lebih kecil dari 100. Ini mengatakan, lihat jawaban yang diedit dengan potongan kode: dengan ini untuk loop Anda dapat mengulangi cv.glmnet 100 kali dan rata-rata kurva kesalahan. Cobalah beberapa kali dan Anda akan melihat bahwa lambda.min tidak akan berubah.
Alice
2
Saya suka bagaimana Anda melakukannya. Saya memiliki loop yang sama tetapi dengan satu pengecualian di akhir: Saya melihat seberapa sering fitur yang berbeda muncul sebagai lawan dari MSE terendah dari semua iterasi. Saya memilih titik potong sewenang-wenang (yaitu menunjukkan iterasi 50/100) dan menggunakan fitur-fitur itu. Penasaran kontras kedua pendekatan itu.
user4673
1
lambdaerror,sincecv
Seperti yang dicatat oleh user4581, fungsi ini bisa gagal karena variabilitas panjangnya cv.glmnet(...)$lambda. Alternatif saya memperbaiki ini: stats.stackexchange.com/a/173895/19676
Max Ghenis
9

λααλα

αλ

Lalu, untuk setiap prediktor saya dapatkan:

  • koefisien rata-rata
  • standar deviasi
  • 5 ringkasan angka (median, kuartil, min, dan maks)
  • Persentase kali berbeda dari nol (mis. memiliki pengaruh)

Dengan cara ini saya mendapatkan deskripsi yang cukup kuat tentang efek prediksi. Setelah Anda memiliki distribusi untuk koefisien, daripada Anda dapat menjalankan hal-hal statistik yang Anda pikir layak untuk mendapatkan CI, nilai p, dll ... tapi saya belum menyelidiki ini.

Metode ini dapat digunakan dengan lebih atau kurang metode pemilihan apa pun yang dapat saya pikirkan.

Bakaburg
sumber
4
Bisakah Anda memposting kode Anda di sini?
rbm
Ya, bisakah Anda memposting kode Anda di sini?
smci
4

Saya akan menambahkan solusi lain, yang menangani bug di @ Alice's karena lambda yang hilang, tetapi tidak memerlukan paket tambahan seperti @Max Ghenis. Terima kasih terutang pada semua jawaban lain - semua orang membuat poin berguna!

lambdas = NULL
for (i in 1:n)
{
    fit <- cv.glmnet(xs,ys)
    errors = data.frame(fit$lambda,fit$cvm)
    lambdas <- rbind(lambdas,errors)
}
# take mean cvm for each lambda
lambdas <- aggregate(lambdas[, 2], list(lambdas$fit.lambda), mean)

# select the best one
bestindex = which(lambdas[2]==min(lambdas[2]))
bestlambda = lambdas[bestindex,1]

# and now run glmnet once more with it
fit <- glmnet(xy,ys,lambda=bestlambda)
Tontonan Bob
sumber
3

Sebagian besar jawaban Alice berfungsi dengan baik, tetapi terkadang kesalahan keluar karena cv.glmnet$lambdaterkadang mengembalikan hasil dengan panjang yang berbeda, misalnya:

Kesalahan dalam rownames <- (tmp, value = c (0.135739830284452, 0.12368107787663,: panjang 'dimnames' [1] tidak sama dengan tingkat array.

OptimLambdadi bawah ini akan berfungsi dalam kasus umum, dan juga lebih cepat dengan memanfaatkan mclapplypemrosesan paralel dan menghindari loop.

Lambdas <- function(...) {
  cv <- cv.glmnet(...)
  return(data.table(cvm=cv$cvm, lambda=cv$lambda))
}

OptimLambda <- function(k, ...) {
  # Returns optimal lambda for glmnet.
  #
  # Args:
  #   k: # times to loop through cv.glmnet.
  #   ...: Other args passed to cv.glmnet.
  #
  # Returns:
  #   Lambda associated with minimum average CV error over runs.
  #
  # Example:
  #   OptimLambda(k=100, y=y, x=x, alpha=alpha, nfolds=k)
  #
  require(parallel)
  require(data.table)
  MSEs <- data.table(rbind.fill(mclapply(seq(k), function(dummy) Lambdas(...))))
  return(MSEs[, list(mean.cvm=mean(cvm)), lambda][order(mean.cvm)][1]$lambda)
}
Max Ghenis
sumber
1

Anda dapat mengontrol keacakan jika Anda secara eksplisit mengatur foldid. Berikut ini contoh untuk CV 5 kali lipat

library(caret)
set.seed(284)
flds <- createFolds(responseDiffs, k = cvfold, list = TRUE, returnTrain = FALSE)
foldids = rep(1,length(responseDiffs))
foldids[flds$Fold2] = 2
foldids[flds$Fold3] = 3
foldids[flds$Fold4] = 4
foldids[flds$Fold5] = 5

Sekarang jalankan cv.glmnet dengan lipatan ini.

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,foldid = foldids)

Anda akan mendapatkan hasil yang sama setiap kali.

Brigitte
sumber