Validação cruzada K-fold ou hold-out para regressão de crista usando R

Estou trabalhando na validação cruzada da previsão dos meus dados com 200 indivíduos e 1000 variáveis. Estou interessado em regressão de cume porque o número de variáveis ​​(eu quero usar) é maior que o número de amostra. Então, eu quero usar estimadores de encolhimento. A seguir, são compostos dados de exemplo:

 #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

Gostaria de fazer o seguinte para validação cruzada -

(1) dividir os dados em duas paradas - use a primeira metade como treinamento e a segunda metade como teste

(2) validação cruzada com dobra em K (digamos 10 vezes ou sugestão em qualquer outra dobra apropriada para o meu caso)

Posso simplesmente amostrar os dados em dois (ganho e teste) e usá-los:

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

Estou usando lm.ridgedo MASSpacote 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)

Eu tenho duas perguntas -

(1) Como posso prever o conjunto de testes e calcular a precisão (como correlação entre o previsto e o real)?

(2) Como posso executar a validação de dobras em K? digamos 10 vezes?

Você pode usar caret pacotes (vinhetas , papel ) para esse tipo de coisa, que pode agrupar vários modelos de aprendizado de máquina ou usar seus próprios modelos personalizados . Como você está interessado em regressão de cume, aqui estão apenas códigos personalizados para regressão de cume, convém adotar sua situação com mais precisão.

Para uma simples divisão de dados:

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

Para validação K-fold e outro tipo de CV, incluindo inicialização padrão

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

Aqui está a discussão sobre como usar a trainfunção. Observe que o método cume depende das elasticnetfunções do pacote (e sua dependência lars, deve ou precisa ser instalada). Se não estiver instalado no sistema, perguntará se você deseja fazê-lo.

o tipo de reamostragem usada, O bootstrap simples é usado por padrão. Para modificar o método de reamostragem, uma função trainControl é usada

O método option controla o tipo de reamostragem e o padrão é "boot". Outro método, "repeatcv", é usado para especificar a validação cruzada repetida em dobras K (e o argumento repetido controla o número de repetições). K é controlado pelo argumento numérico e o padrão é 10.

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

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

plot(ridgefit)

Para previsões:

plsClasses <- predict(ridgeFit, newdata = testing)

Esta é uma extensão da sugestão de Frank nos comentários. Dr. Harrel, por favor, corrija se estou errado (aprecie as correções).

Seus dados:

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

Instale o rmspacote e carregue-o.

require(rms)

ols A função é usada para Estimativa de modelo linear usando mínimos quadrados ordinários, onde é possível especificar o termo da penalidade.

Como sugerido abaixo nos comentários, adicionei a petracefunção. Esta função rastreia AIC e BIC vs Penalidade.

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

Nota importante Eu não poderia usar todas as 1000 variáveis, pois o programa reclama se o número de variáveis ​​exceder 100. Também a y~.designação da fórmula do tipo não funcionou. Portanto, veja acima a maneira de fazer o mesmo objeto de fórmula de criaçãofrm

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)

"Para um ajuste comum e não penalizado de lrm ou ols e para um vetor ou lista de multas, se encaixa em uma série de modelos logísticos ou lineares usando a estimativa de máxima verossimilhança penalizada e salva os graus efetivos de liberdade", Schwarz Bayesian Critério de Informação (BIC) e AIC corrigido de Hurvich e Tsai (AIC_c). Opcionalmente, o pentrace pode usar a função nlminb para resolver o fator de penalidade ideal ou a combinação de fatores que penalizam diferentes tipos de termos no modelo ". do rmsmanual do pacote.

calibrateA função destina-se a reamostragem de calibração de modelos e usa inicialização ou validação cruzada para obter estimativas corrigidas de viés (correção de ajuste excessivo) dos valores previstos versus os observados, com base na definição de subconjuntos em intervalos. A validatefunção faz a reamostragem da validação de um modelo de regressão, com ou sem exclusão de variável de etapa para trás. B = número de repetições. Para method = "crossvalidation", é o número de grupos de observações omitidas

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

Você pode usar a Predictfunção para calcular valores previstos e limites de confiança. Não tenho certeza se isso funciona na situação de teste.

O pacote R glmnet( vinheta ) possui uma função de wrapper que faz exatamente o que você deseja, chamada cv.glmnet( doc ). Eu usei ontem ontem, funciona como um sonho.