Walidacja krzyżowa K-hold lub hold-out dla regresji grzbietu za pomocą R.

Pracuję nad weryfikacją krzyżową prognoz moich danych z 200 podmiotami i 1000 zmiennymi. Interesuje mnie regresja grzbietu, ponieważ liczba zmiennych (chcę użyć) jest większa niż liczba próbek. Więc chcę użyć estymatorów skurczu. Oto przykładowe dane:

 #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

Chciałbym wykonać następujące czynności w celu weryfikacji krzyżowej -

(1) podziel dane na dwie części - wykorzystaj pierwszą połowę jako trening, a drugą połowę jako test

(2) Walidacja krzyżowa K-krotnie (mile widziane jest 10-krotnie lub sugestia przy każdym innym odpowiednim foldem dla mojej sprawy)

Mogę po prostu próbkować dane na dwie części (zbieranie i testowanie) i używać ich:

# 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,]   

Korzystam lm.ridgez MASSpakietu 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)

Mam dwa pytania -

(1) Jak mogę przewidzieć zestaw testów i obliczyć dokładność (jako korelację prognozy z rzeczywistością)?

(2) Jak mogę przeprowadzić walidację K-fold? powiedz 10-krotnie?

Możesz użyć caret pakietu (winiety , papier ) do tego rodzaju rzeczy, które mogą owijać wiele modeli uczenia maszynowego lub możesz używać własnych modeli niestandardowych . Ponieważ interesuje Cię regresja kalenicy, tutaj są tylko niestandardowe kody regresji kalenicy, możesz chcieć bardziej precyzyjnie dostosować się do swojej sytuacji.

Dla prostego podziału danych:

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

Do weryfikacji K-fold i innego rodzaju CV, w tym domyślnego rozruchu

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

Oto dyskusja na temat korzystania z trainfunkcji. Uwaga: metoda ridge zależy od elasticnetfunkcji pakietu (i jej zależności od lars, powinna lub powinna zostać zainstalowana). Jeśli nie jest zainstalowany w systemie, zapyta, czy chcesz to zrobić.

rodzaj zastosowanego ponownego próbkowania, domyślnie używany jest prosty pasek startowy. Aby zmodyfikować metodę ponownego próbkowania, używana jest funkcja trainControl

Metoda opcjonalna kontroluje typ ponownego próbkowania i domyślnie „rozruch”. Inna metoda, „repeatcv”, służy do określania powtarzanej krzyżowej weryfikacji K-krotności (a argument repeats kontroluje liczbę powtórzeń). K jest kontrolowany przez argument liczby i domyślnie wynosi 10.

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

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

plot(ridgefit)

Do prognoz:

plsClasses <- predict(ridgeFit, newdata = testing)

Jest to rozszerzenie sugestii Franka w komentarzach. Dr Harrel, proszę o poprawienie, jeśli się mylę (doceniam poprawki).

Twoje dane:

#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)

Zainstaluj rmspakiet i załaduj go.

require(rms)

ols funkcja służy do szacowania modelu liniowego przy użyciu zwykłych najmniejszych kwadratów, w których można określić okres kary.

Jak sugerowałem poniżej w komentarzach, dodałem petracefunkcję. Ta funkcja śledzi AIC i BIC vs. Kara.

# 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="+")))

Ważna uwaga Nie mogłem użyć wszystkich 1000 zmiennych, ponieważ program narzeka, jeśli liczba zmiennych przekracza 100. y~.Nie działało też oznaczenie formuły typu. Zobacz więc powyższy sposób tworzenia tego samego obiektu formułyfrm

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)

„W przypadku zwykłego niezapenalizowanego dopasowania od lrm lub ols oraz do wektora lub listy kar pasuje do szeregu modeli logistycznych lub liniowych z zastosowaniem karanego oszacowania maksymalnego prawdopodobieństwa i oszczędza efektywne stopnie swobody, Akaike Information Criterion (AIC), Schwarz Bayesian Kryterium informacyjne (BIC) oraz AIC poprawione przez Hurvicha i Tsai (AIC_c). Opcjonalnie pentrace może użyć funkcji nlminb do rozwiązania optymalnego współczynnika kary lub kombinacji czynników karających różne rodzaje terminów w modelu. ” z rmsinstrukcji pakietu.

calibratefunkcja służy do kalibracji modelu ponownego próbkowania i używa ładowania początkowego lub walidacji krzyżowej w celu uzyskania skorygowanych uprzedzeń (przeregulowanych) oszacowań prognozowanych i obserwowanych wartości w oparciu o podzbiór prognoz na przedziały. validateFunkcja jest resampling walidację modelu regresji, z lub bez delecji w tył zmiennej obniżającego. B = liczba powtórzeń. Dla metody = „crossvalidation” oznacza liczbę grup pominiętych obserwacji

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

Możesz użyć Predictfunkcji do obliczenia przewidywanych wartości i granic ufności. Nie jestem pewien, czy to działa w sytuacji testowej.

Pakiet R glmnet( winieta ) ma funkcję otoki, która wykonuje dokładnie to, co chcesz, o nazwie cv.glmnet( doc ). Właśnie go użyłem wczoraj, działa jak sen.