Zmienność wyników cv.glmnet

Używam cv.glmnetdo znajdowania predyktorów. Konfiguracja, której używam jest następująca:

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

Aby upewnić się, że wyniki są powtarzalne ja set.seed(1). Wyniki są bardzo zmienne. Uruchomiłem dokładnie ten sam kod 100, aby zobaczyć, jak zmienne były wyniki. W biegach 98/100 zawsze wybierano jeden konkretny predyktor (czasem tylko sam); inne predyktory zostały wybrane (współczynnik był niezerowy) zwykle 50/100 razy.

Mówi mi więc, że za każdym razem, gdy przeprowadzana jest walidacja krzyżowa, prawdopodobnie wybierze inną najlepszą lambda, ponieważ początkowa randomizacja fałdów ma znaczenie. Inni widzieli ten problem ( wyniki CV.glmnet ), ale nie ma sugerowanego rozwiązania.

Myślę, że może ten, który pokazuje 98/100, jest prawdopodobnie dość mocno skorelowany ze wszystkimi innymi? Wyniki należy ustabilizować jeśli po prostu uruchomić LOOCV ( $\text{fold-size} = n$ ), ale jestem ciekaw, dlaczego są one tak zmiennej kiedy $\text{nfold} < n$ .

Chodzi o to, że w cv.glmnetK fałdy („części”) są wybierane losowo.

W krzyżowej weryfikacji K-fałdów zestaw danych jest podzielony na części , a części K - 1 służą do przewidywania części K (odbywa się to razy K , za każdym razem przy użyciu innej części K ). Odbywa się to dla wszystkich lambd, i to ten, który daje najmniejszy błąd weryfikacji krzyżowej.$K$$K-1$$K$$K$lambda.min

Dlatego przy użyciu wyniki się nie zmieniają: każda grupa składa się z jednej, więc nie ma dużego wyboru dla grup $nfolds = n$$K$

Z cv.glmnet()podręcznika referencyjnego:

Zauważ również, że wyniki cv.glmnet są losowe, ponieważ fałdy są wybierane losowo. Użytkownicy mogą zmniejszyć tę losowość, uruchamiając wiele razy cv.glmnet i uśredniając krzywe błędów.

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

MSE to ramka danych zawierająca wszystkie błędy dla wszystkich lambd (dla 100 przebiegów), lambda.minto twoja lambda z minimalnym średnim błędem.

$\lambda$$\alpha$$\alpha$$\lambda$$\alpha$

$\alpha$$\lambda$

Następnie dla każdego predyktora otrzymuję:

• średni współczynnik
• odchylenie standardowe
• Podsumowanie 5 liczb (mediana, kwartyle, min i maks)
• procent razy różni się od zera (tzn. ma wpływ)

W ten sposób uzyskuję dość solidny opis działania predyktora. Gdy masz już rozkłady dla współczynników, nie możesz uruchomić żadnych statystyk, które Twoim zdaniem warto uzyskać CI, wartości p itp., Ale jeszcze tego nie zbadałem.

Ta metoda może być używana z mniej więcej dowolną metodą selekcji, o której mogę myśleć.

Dodam inne rozwiązanie, które obsługuje błąd w @ Alice z powodu brakujących lambd, ale nie wymaga dodatkowych pakietów, takich jak @Max Ghenis. Podziękowania należą się wszystkim pozostałym odpowiedziom - każdy robi użyteczne punkty!

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)

Odpowiedź Alicji działa dobrze w większości przypadków, ale czasami pomija błędy z powodu cv.glmnet$lambdaczasami zwracających wyniki o różnej długości, np .:

Błąd w nazwach <- (tmp, wartość = c (0.135739830284452, 0.12368107787663,: długość „dimnames” [1] nie jest równa zasięgowi tablicy).

OptimLambdaponiżej powinno działać w ogólnym przypadku, a także jest szybsze, wykorzystując mclapplydo równoległego przetwarzania i unikania pętli.

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

You can control the randomness if you explicitly set foldid. Here an example for 5-fold CV

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

Now run cv.glmnet with these foldids.

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

You will get the same results each time.