K-fold atau validasi silang tahan untuk regresi ridge menggunakan R

Saya sedang mengerjakan validasi silang dari prediksi data saya dengan 200 subjek dan 1000 variabel. Saya tertarik regresi ridge karena jumlah variabel (saya ingin menggunakan) lebih besar dari jumlah sampel. Jadi saya ingin menggunakan estimator penyusutan. Berikut ini adalah contoh data:

 #random population of 200 subjects with 1000 variables 
    M <- matrix(rep(0,200*100),200,1000)
    for (i in 1:200) {
    set.seed(i)
      M[i,] <- ifelse(runif(1000)<0.5,-1,1)
    }
    rownames(M) <- 1:200

    #random yvars 
    set.seed(1234)
    u <- rnorm(1000)
    g <- as.vector(crossprod(t(M),u))
    h2 <- 0.5 
    set.seed(234)
    y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

    myd <- data.frame(y=y, M)
myd[1:10,1:10]

y X1 X2 X3 X4 X5 X6 X7 X8 X9
1   -7.443403 -1 -1  1  1 -1  1  1  1  1
2  -63.731438 -1  1  1 -1  1  1 -1  1 -1
3  -48.705165 -1  1 -1 -1  1  1 -1 -1  1
4   15.883502  1 -1 -1 -1  1 -1  1  1  1
5   19.087484 -1  1  1 -1 -1  1  1  1  1
6   44.066119  1  1 -1 -1  1  1  1  1  1
7  -26.871182  1 -1 -1 -1 -1  1 -1  1 -1
8  -63.120595 -1 -1  1  1 -1  1 -1  1  1
9   48.330940 -1 -1 -1 -1 -1 -1 -1 -1  1
10 -18.433047  1 -1 -1  1 -1 -1 -1 -1  1

Saya ingin melakukan yang berikut untuk validasi silang -

(1) membagi data menjadi dua penghentian - gunakan babak pertama sebagai latihan dan babak kedua sebagai tes

(2) validasi silang K-fold (katakanlah 10 kali lipat atau saran pada lipatan lain yang sesuai untuk kasus saya dipersilakan)

Saya hanya bisa sampel data menjadi dua (memperoleh dan menguji) dan menggunakannya:

# using holdout (50% of the data) cross validation 
training.id <- sample(1:nrow(myd), round(nrow(myd)/2,0), replace = FALSE)
test.id <- setdiff(1:nrow(myd), training.id)

 myd_train <- myd[training.id,]
 myd_test  <- myd[test.id,]   

Saya menggunakan lm.ridgedari MASSpaket R.

library(MASS)
out.ridge=lm.ridge(y~., data=myd_train, lambda=seq(0, 100,0.001))
plot(out.ridge)
select(out.ridge)

lam=0.001
abline(v=lam)

out.ridge1 =lm.ridge(y~., data=myd_train, lambda=lam)
hist(out.ridge1$coef)
    out.ridge1$ym
hist(out.ridge1$xm)

Saya punya dua pertanyaan -

(1) Bagaimana saya bisa memprediksi set tes dan menghitung akurasi (sebagai korelasi prediksi vs aktual)?

(2) Bagaimana saya bisa melakukan validasi K-fold? katakan 10 kali lipat?

Anda dapat menggunakan caret paket (sketsa , kertas ) untuk jenis hal ini, yang dapat membungkus sejumlah model pembelajaran mesin atau Anda dapat menggunakan model khusus Anda sendiri . Karena Anda tertarik pada regresi ridge di sini hanya kode kustom untuk regresi ridge, Anda mungkin ingin mengadopsi situasi Anda lebih tepat.

Untuk pemisahan data sederhana:

set.seed(107)
# stratified random split of the data
inTrain <- createDataPartition(y = myd$y, p = .5,list = FALSE)
training <- myd[ inTrain,]
testing <- myd[-inTrain,]

Untuk validasi K-fold dan tipe CV lainnya termasuk boot default

ridgeFit1 <- train(y ~ ., data = training,method = 'ridge', 
preProc = c("center", "scale"), metric = "ROC")
plot(ridgeFit1)

Berikut ini diskusi tentang cara menggunakan trainfungsi. Catatan metode punggungan tergantung pada elasticnetfungsi-fungsi paket (dan ketergantungannya pada lars, harus atau perlu diinstal). Jika tidak diinstal dalam sistem akan menanyakan apakah Anda ingin melakukannya.

jenis resampling yang digunakan, Bootstrap sederhana digunakan secara default. Untuk memodifikasi metode resampling, fungsi trainControl digunakan

Metode opsi mengontrol jenis resampling dan standar untuk "boot". Metode lain, "repeatcv", digunakan untuk menentukan validasi K-fold cross-validation (dan argumen yang diulang mengontrol jumlah pengulangan). K dikendalikan oleh argumen angka dan default ke 10.

 ctrl <- trainControl(method = "repeatedcv", repeats = 5)

 ridgeFit <- train(y ~ ., data = training,method = 'ridge',
preProc = c("center", "scale"),trControl = ctrl, metric = "ROC")

plot(ridgefit)

Untuk prediksi:

plsClasses <- predict(ridgeFit, newdata = testing)

Ini adalah perpanjangan dari saran oleh Frank di komentar. Harrel tolong perbaiki jika saya salah (menghargai koreksi).

Data Anda:

#random population of 200 subjects with 1000 variables 
    M <- matrix(rep(0,200*100),200,1000)
    for (i in 1:200) {
    set.seed(i)
      M[i,] <- ifelse(runif(1000)<0.5,-1,1)
    }
    rownames(M) <- 1:200

    #random yvars 
    set.seed(1234)
    u <- rnorm(1000)
    g <- as.vector(crossprod(t(M),u))
    h2 <- 0.5 
    set.seed(234)
    y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

    myd <- data.frame(y=y, M)

Instal rmspaket dan muat.

require(rms)

ols fungsi digunakan untuk Estimasi Model Linier Menggunakan Kuadrat Terkecil Biasa di mana dapat menentukan jangka waktu penalti.

Seperti yang disarankan di bawah ini di komentar saya menambahkan petracefungsi. Fungsi ini melacak AIC dan BIC vs Penalty.

# using holdout (50% of the data) cross validation 
training.id <- sample(1:nrow(myd), round(nrow(myd)/2,0), replace = FALSE)
test.id <- setdiff(1:nrow(myd), training.id)

 myd_train <- myd[training.id,]
 myd_test  <- myd[test.id,] 

frm <- as.formula(paste("y~",paste(names(myd_train)[2:100],collapse="+")))

Catatan penting Saya tidak bisa menggunakan semua 1000 variabel karena program mengeluh jika jumlah variabel melebihi 100. Juga y~.ketik penunjukan rumus tidak berfungsi. Jadi lihat di atas cara melakukan objek rumus yang sama membuatfrm

f <- ols(frm, data = myd_train, method="qr", x=TRUE, y=TRUE)


p <- pentrace(f, seq(.2,1,by=.05))

Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x),  : 
'data' must be of a vector type, was 'NULL'

 plot(p)

"Untuk kecocokan biasa yang tidak diloloskan dari lrm atau ols dan untuk vektor atau daftar penalti, cocok dengan serangkaian model logistik atau linier menggunakan estimasi kemungkinan maksimum yang dihukum, dan menyimpan derajat kebebasan efektif, Akaike Information Criterion (AIC), Schwarz Bayesian Kriteria Informasi (BIC), dan Hurvich dan AAI Tsai yang dikoreksi (AIC_c). Secara opsional pentrace dapat menggunakan fungsi nlminb untuk menyelesaikan faktor penalti optimal atau kombinasi faktor yang menghukum berbagai jenis istilah dalam model. " dari rmsmanual paket.

calibratefungsinya adalah untuk Resampling Model Calibration dan Menggunakan bootstrap atau cross-validation untuk mendapatkan estimasi bias-koreksi (overfitting-koreksi) dari nilai-nilai yang diprediksi vs yang diobservasi berdasarkan subset prediksi menjadi interval. The validatefungsi tersebut resampling validasi dari model regresi, dengan atau tanpa mundur step-down variabel penghapusan. B = jumlah pengulangan. Untuk metode = "crossvalidation", adalah jumlah kelompok pengamatan yang dihilangkan

cal <- calibrate(f, method = "cross validation", B=20)  
plot(cal)

Anda dapat menggunakan Predictfungsi untuk menghitung nilai yang diprediksi dan batas kepercayaan. Saya tidak yakin saya ini bekerja dalam situasi ujian.

Paket R glmnet( sketsa ) memiliki fungsi wrapper yang melakukan persis apa yang Anda inginkan, disebut cv.glmnet( doc ). Saya hanya menggunakannya kemarin, itu berfungsi seperti mimpi.