Validation croisée K-fold ou hold-out pour la régression de crête à l'aide de R

Je travaille sur la validation croisée de la prédiction de mes données avec 200 sujets et 1000 variables. Je suis intéressé par la régression des crêtes car le nombre de variables (que je veux utiliser) est supérieur au nombre d'échantillons. Je veux donc utiliser des estimateurs de retrait. Voici des exemples de données:

 #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

Je voudrais faire ce qui suit pour la validation croisée -

(1) diviser les données en deux - utilisez la première moitié comme formation et la seconde moitié comme test

(2) Validation croisée K-fold (disons 10 fois ou une suggestion sur tout autre pli approprié pour mon cas est la bienvenue)

Je peux simplement échantillonner les données en deux (gagner et tester) et les utiliser:

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

J'utilise lm.ridgedepuis le MASSpackage 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)

J'ai deux questions -

(1) Comment puis-je prédire l'ensemble de test et calculer la précision (en tant que corrélation entre le prévu et le réel)?

(2) Comment puis-je effectuer la validation K-fold? dites 10 fois?

Vous pouvez utiliser un caret package (vignettes , papier ) pour ce type de choses, qui peut envelopper un certain nombre de modèles d'apprentissage automatique ou vous pouvez utiliser vos propres modèles personnalisés . Comme vous êtes intéressé par la régression des crêtes, voici juste des codes personnalisés pour la régression des crêtes, vous voudrez peut-être vous adapter plus précisément à votre situation.

Pour une répartition simple des données:

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

Pour la validation K-fold et d'autres types de CV, y compris le démarrage par défaut

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

Voici une discussion sur la façon d'utiliser la trainfonction. Notez que la méthode ridge dépend des fonctions du package elasticnet(et de sa dépendance sur lars, devrait ou doit être installée). S'il n'est pas installé dans le système, il vous demandera si vous souhaitez le faire.

le type de rééchantillonnage utilisé, le bootstrap simple est utilisé par défaut.Pour modifier la méthode de rééchantillonnage, une fonction trainControl est utilisée

La méthode d'option contrôle le type de rééchantillonnage et par défaut "boot". Une autre méthode, "repeatcv", est utilisée pour spécifier la validation croisée multipliée par K (et l'argument se répète contrôle le nombre de répétitions). K est contrôlé par l'argument nombre et par défaut à 10.

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

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

plot(ridgefit)

Pour les prédictions:

plsClasses <- predict(ridgeFit, newdata = testing)

Il s'agit d'une extension de la suggestion de Frank dans les commentaires. Dr Harrel, veuillez corriger si je me trompe (veuillez apprécier les corrections).

Vos données:

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

Installez le rmspackage et chargez-le.

require(rms)

ols La fonction est utilisée pour l'estimation du modèle linéaire en utilisant les moindres carrés ordinaires où peut spécifier le terme de pénalité.

Comme suggéré ci-dessous dans les commentaires, j'ai ajouté une petracefonction. Cette fonction trace AIC et 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="+")))

Remarque importante Je n'ai pas pu utiliser les 1000 variables car le programme se plaint si le nombre de variables dépasse 100. De plus, la y~.désignation de formule de type n'a pas fonctionné. Donc, voyez ci-dessus la façon de faire la même chose en créant un objet de formulefrm

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)

"Pour un ajustement ordinaire non pénalisé à partir de lrm ou ols et pour un vecteur ou une liste de pénalités, correspond à une série de modèles logistiques ou linéaires utilisant une estimation du maximum de vraisemblance pénalisée et enregistre les degrés de liberté effectifs, Akaike Information Criterion (AIC), Schwarz Bayesian Information Criterion (BIC) et AIC corrigé de Hurvich et Tsai (AIC_c). En option, pentrace peut utiliser la fonction nlminb pour résoudre le facteur de pénalité optimal ou la combinaison de facteurs pénalisant différents types de termes dans le modèle. " du rmsmanuel du paquet.

calibrateLa fonction est pour l'étalonnage du modèle de rééchantillonnage et utilise le bootstrap ou la validation croisée pour obtenir des estimations corrigées du biais (sur-ajustées) des valeurs prédites par rapport aux valeurs observées sur la base d'un sous-ensemble des prédictions dans les intervalles. La validatefonction effectue un rééchantillonnage de validation d'un modèle de régression, avec ou sans suppression de variable descendante descendante. B = nombre de répétitions. Pour method = "crossvalidation", est le nombre de groupes d'observations omises

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

Vous pouvez utiliser la Predictfonction pour calculer les valeurs prédites et les limites de confiance. Je ne suis pas sûr que cela fonctionne en situation de test.

Le package R glmnet( vignette ) a une fonction wrapper qui fait exactement ce que vous voulez, appelée cv.glmnet( doc ). Je l'ai utilisé hier, ça marche comme un rêve.