K-fache oder Hold-out-Kreuzvalidierung für die Gratregression mit R.

Ich arbeite an einer Kreuzvalidierung der Vorhersage meiner Daten mit 200 Probanden und 1000 Variablen. Ich bin an einer Ridge-Regression interessiert, da die Anzahl der Variablen (die ich verwenden möchte) größer ist als die Anzahl der Stichproben. Ich möchte also Schrumpfungsschätzer verwenden. Die folgenden Beispieldaten bestehen aus:

 #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

Ich möchte Folgendes zur Kreuzvalidierung tun -

(1) Teilen Sie die Daten in zwei Pausen auf - verwenden Sie die erste Hälfte als Training und die zweite Hälfte als Test

(2) K-fache Kreuzvalidierung (z. B. 10-fach oder Vorschlag für eine andere geeignete Falte für meinen Fall sind willkommen)

Ich kann die Daten einfach in zwei Teile zerlegen (gewinnen und testen) und sie verwenden:

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

Ich benutze lm.ridgevon MASSR-Paket.

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)

Ich habe zwei Fragen -

(1) Wie kann ich den Testsatz vorhersagen und die Genauigkeit berechnen (als Korrelation zwischen vorhergesagt und tatsächlich)?

(2) Wie kann ich eine K-fache Validierung durchführen? 10-fach sagen?

Sie können für diese Art von Dingen ein caret Paket (Vignnetten , Papier ) verwenden, das eine Reihe von Modellen für maschinelles Lernen umschließen kann , oder Sie können Ihre eigenen benutzerdefinierten Modelle verwenden . Da Sie an einer Ridge-Regression interessiert sind, handelt es sich hier nur um benutzerdefinierte Codes für die Ridge-Regression. Vielleicht möchten Sie diese genauer auf Ihre Situation anwenden.

Für eine einfache Aufteilung der Daten:

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

Für die K-Fold-Validierung und andere Arten von Lebensläufen, einschließlich Standardstart

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

Hier finden Sie eine Diskussion zur Verwendung der trainFunktion. Beachten Sie, dass die Ridge-Methode von den Paketfunktionen abhängt elasticnet(und von der Abhängigkeit lars, ob sie installiert werden muss oder muss). Wenn nicht im System installiert, werden Sie gefragt, ob Sie dies möchten.

Die Art des verwendeten Resamplings. Standardmäßig wird der einfache Bootstrap verwendet. Um die Resampling-Methode zu ändern, wird eine trainControl-Funktion verwendet

Die Optionsmethode steuert die Art des Resamplings und ist standardmäßig "boot". Eine andere Methode, "repeatcv", wird verwendet, um die wiederholte K-fache Kreuzvalidierung anzugeben (und das Argument "Wiederholungen" steuert die Anzahl der Wiederholungen). K wird durch das Zahlenargument gesteuert und ist standardmäßig 10.

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

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

plot(ridgefit)

Für Vorhersagen:

plsClasses <- predict(ridgeFit, newdata = testing)

Dies ist eine Erweiterung des Vorschlags von Frank in den Kommentaren. Dr. Harrel, bitte korrigieren Sie, wenn ich falsch liege (schätzen Sie Korrekturen).

Deine Daten:

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

Installieren Sie das rmsPaket und laden Sie es.

require(rms)

ols Die Funktion wird für die lineare Modellschätzung unter Verwendung gewöhnlicher kleinster Quadrate verwendet, in denen der Strafbegriff angegeben werden kann.

Wie unten in den Kommentaren vorgeschlagen, habe ich eine petraceFunktion hinzugefügt . Diese Funktion verfolgt AIC und BIC gegen 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="+")))

Wichtiger Hinweis Ich konnte nicht alle 1000 Variablen verwenden, da sich das Programm beschwert, wenn die Anzahl der Variablen 100 überschreitet. Auch die y~.Typformelbezeichnung funktionierte nicht. Siehe oben, wie Sie dasselbe Formelobjekt erstellenfrm

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)

"Für eine gewöhnliche ungestrafte Anpassung von lrm oder ols und für einen Vektor oder eine Liste von Strafen passt eine Reihe von logistischen oder linearen Modellen unter Verwendung der bestraften Maximum-Likelihood-Schätzung an und speichert die effektiven Freiheitsgrade, Akaike Information Criterion (AIC), Schwarz Bayesian Information Criterion (BIC) und Hurvich und Tsais korrigierter AIC (AIC_c). Optional kann Pentrace die Funktion nlminb verwenden, um nach dem optimalen Straffaktor oder einer Kombination von Faktoren zu suchen, die verschiedene Arten von Begriffen im Modell benachteiligen. " aus dem rmsPakethandbuch.

calibrateDie Funktion dient zur Resampling-Modellkalibrierung und verwendet Bootstrapping oder Kreuzvalidierung, um vorspannungskorrigierte (überanpassungskorrigierte) Schätzungen von vorhergesagten vs. beobachteten Werten basierend auf Teilmengenvorhersagen in Intervallen zu erhalten. Die validateFunktion führt eine erneute Abtastvalidierung eines Regressionsmodells mit oder ohne Rückwärts-Step-Down-Variablenlöschung durch. B = Anzahl der Wiederholungen. Für method = "crossvalidation" ist die Anzahl der Gruppen ausgelassener Beobachtungen

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

Mit der PredictFunktion können Sie vorhergesagte Werte und Konfidenzgrenzen berechnen. Ich bin mir nicht sicher, ob dies in einer Testsituation funktioniert.

Das R-Paket glmnet( Vignette ) verfügt über eine Wrapper-Funktion namens cv.glmnet( doc ) , die genau das tut, was Sie wollen . Ich habe es erst gestern benutzt, es funktioniert wie ein Traum.