Variabilidade nos resultados do cv.glmnet

18

Estou usando cv.glmnetpara encontrar preditores. A configuração que eu uso é a seguinte:

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,nfolds=cvfold)
bestlambda<-lassoResults$lambda.min

results<-predict(lassoResults,s=bestlambda,type="coefficients")

choicePred<-rownames(results)[which(results !=0)]

Garantir que os resultados sejam reprodutíveis set.seed(1). Os resultados são altamente variáveis. Corri exatamente o mesmo código 100 para ver como os resultados eram variáveis. Nas corridas 98/100, sempre havia um preditor em particular selecionado (às vezes por conta própria); outros preditores foram selecionados (o coeficiente era diferente de zero) geralmente 50/100 vezes.

Por isso, diz-me que cada vez que a validação cruzada estiver em execução, provavelmente irá selecionar um melhor lambda diferente, porque a randomização inicial das dobras é importante. Outros viram esse problema ( resultados do CV.glmnet ), mas não há uma solução sugerida.

Eu estou pensando que talvez aquele que aparece 98/100 provavelmente esteja bastante correlacionado com todos os outros? Os resultados não estabilizar, se eu LOOCV apenas correr ( ), mas estou curioso por que eles são tão variáveis quando \ text {nfold} <n .fold-size=nnfold<n

user4673
fonte
1
Para ser claro, você quer dizer que você set.seed(1)corre uma vez cv.glmnet()100 vezes? Essa não é uma ótima metodologia para reprodutibilidade; melhor para a set.seed()direita antes de cada corrida, ou mantenha as dobras constantes entre as corridas. Cada uma das chamadas para cv.glmnet()está chamando sample()N vezes. Portanto, se o comprimento dos seus dados mudar, a reprodubilidade será alterada.
smci 24/02

Respostas:

14

O ponto aqui é que nas cv.glmnetdobras K ("partes") são escolhidas aleatoriamente.

Na validação cruzada de dobras K, o conjunto de dados é dividido em partes partes são usadas para prever a parte K-ésima (isso é feito vezes, usando uma parte diferente a cada vez). Isso é feito para todas as lambdas e é a que fornece o menor erro de validação cruzada.K - 1 K KKK1KKlambda.min

É por isso que, quando você usa os resultados não mudam: cada grupo é composto de um, portanto, não há muita escolha para os gruposKnfolds=nK

No cv.glmnet()manual de referência:

Observe também que os resultados do cv.glmnet são aleatórios, pois as dobras são selecionadas aleatoriamente. Os usuários podem reduzir essa aleatoriedade executando o cv.glmnet várias vezes e calculando a média das curvas de erro.

### cycle for doing 100 cross validations
### and take the average of the mean error curves
### initialize vector for final data.frame with Mean Standard Errors
MSEs <- NULL
for (i in 1:100){
                 cv <- cv.glmnet(y, x, alpha=alpha, nfolds=k)  
                 MSEs <- cbind(MSEs, cv$cvm)
             }
  rownames(MSEs) <- cv$lambda
  lambda.min <- as.numeric(names(which.min(rowMeans(MSEs))))

MSEs é o quadro de dados que contém todos os erros para todas as lambdas (para as 100 execuções), lambda.miné sua lambda com erro médio mínimo.

Alice
fonte
O que mais me preocupa é que a seleção de n realmente parece importar às vezes. Devo confiar em resultados que podem ser tão variáveis? Ou devo classificá-lo como incompleto, mesmo que eu o execute várias vezes?
user4673
1
Dependendo do tamanho da amostra, você deve escolher n para ter pelo menos 10 observações por grupo. Portanto, é melhor diminuir o n padrão (= 10) se você tiver um tamanho de amostra menor que 100. Dito isso, consulte a resposta editada com o trecho de código: com isso para o loop, você pode repetir o cv.glmnet 100 vezes e calcular a média dos curvas de erro. Experimente algumas vezes e você verá que o lambda.min não muda.
Alice
2
Eu gosto de como você fez isso. Eu tenho o mesmo loop, mas com uma exceção no final: eu olho para a frequência com que diferentes recursos aparecem em oposição ao MSE mais baixo de todas as iterações. Eu escolho um ponto de corte arbitrário (ou seja, mostro 50/100 iterações) e uso esses recursos. Curioso contraste entre as duas abordagens.
user4673
1
Este rownames (MPEs) <- cv lambda no meu caso é mais do que MPEs (I assumir o seu devido à convergência ...)lambdaerror,sincecv
user4581
Como o usuário4581 observou, esta função pode falhar devido à variabilidade no comprimento de cv.glmnet(...)$lambda. Meus correções alternativas isto: stats.stackexchange.com/a/173895/19676
Max Ghenis
9

Ultimamente, enfrentei o mesmo problema. Tentei repetir o CV várias vezes, como 100, 200, 1000 no meu conjunto de dados, tentando encontrar o melhor e (estou usando uma rede elástica). Mas, mesmo se eu criar 3 testes cv, cada um com 1000 iterações que calculam a média dos mínimos MSEs para cada , recebo três casais diferentes ( , ).α α λ αλααλα

Não tocarei no problema aqui, mas decidi que minha melhor solução não é calcular a média dos min MSEs, mas extrair os coeficientes para cada iteração best e depois tratá-los como uma distribuição de valores (uma variável aleatória).λαλ

Então, para cada preditor, recebo:

  • coeficiente médio
  • desvio padrão
  • Resumo de 5 números (mediana, quartis, mínimo e máximo)
  • A porcentagem de vezes é diferente de zero (ou seja, tem influência)

Dessa forma, recebo uma descrição bastante sólida do efeito do preditor. Depois de ter distribuições para os coeficientes, você pode executar qualquer coisa estatística que considere valer para obter IC, valores de p, etc ... mas ainda não investiguei isso.

Este método pode ser usado com mais ou menos qualquer método de seleção que eu possa pensar.

Bakaburg
fonte
4
Você pode postar seu código aqui, por favor?
Rbm #
Sim, você pode postar seu código aqui?
smci
4

Vou adicionar outra solução, que lida com o bug no @ Alice devido à falta de lambdas, mas não requer pacotes extras como o @Max Ghenis. Agradecemos a todas as outras respostas - todo mundo faz pontos úteis!

lambdas = NULL
for (i in 1:n)
{
    fit <- cv.glmnet(xs,ys)
    errors = data.frame(fit$lambda,fit$cvm)
    lambdas <- rbind(lambdas,errors)
}
# take mean cvm for each lambda
lambdas <- aggregate(lambdas[, 2], list(lambdas$fit.lambda), mean)

# select the best one
bestindex = which(lambdas[2]==min(lambdas[2]))
bestlambda = lambdas[bestindex,1]

# and now run glmnet once more with it
fit <- glmnet(xy,ys,lambda=bestlambda)
Sideshow Bob
fonte
3

A resposta de Alice funciona bem na maioria dos casos, mas às vezes ocorre erros devido às cv.glmnet$lambdavezes retornando resultados de diferentes tamanhos, por exemplo:

Erro nos nomes de usuários <- (tmp, valor = c (0.135739830284452, 0.12368107787663,: comprimento de 'dimnames' [1] diferente da extensão da matriz).

OptimLambdaabaixo deve funcionar no caso geral e também é mais rápido, aproveitando o mclapplyprocessamento paralelo e evitando loops.

Lambdas <- function(...) {
  cv <- cv.glmnet(...)
  return(data.table(cvm=cv$cvm, lambda=cv$lambda))
}

OptimLambda <- function(k, ...) {
  # Returns optimal lambda for glmnet.
  #
  # Args:
  #   k: # times to loop through cv.glmnet.
  #   ...: Other args passed to cv.glmnet.
  #
  # Returns:
  #   Lambda associated with minimum average CV error over runs.
  #
  # Example:
  #   OptimLambda(k=100, y=y, x=x, alpha=alpha, nfolds=k)
  #
  require(parallel)
  require(data.table)
  MSEs <- data.table(rbind.fill(mclapply(seq(k), function(dummy) Lambdas(...))))
  return(MSEs[, list(mean.cvm=mean(cvm)), lambda][order(mean.cvm)][1]$lambda)
}
Max Ghenis
fonte
1

Você pode controlar a aleatoriedade se definir explicitamente foldid. Aqui está um exemplo de CV 5 vezes

library(caret)
set.seed(284)
flds <- createFolds(responseDiffs, k = cvfold, list = TRUE, returnTrain = FALSE)
foldids = rep(1,length(responseDiffs))
foldids[flds$Fold2] = 2
foldids[flds$Fold3] = 3
foldids[flds$Fold4] = 4
foldids[flds$Fold5] = 5

Agora execute o cv.glmnet com esses foldids.

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,foldid = foldids)

Você obterá os mesmos resultados todas as vezes.

Brigitte
fonte