Jak można empirycznie wykazać w R, jakie metody walidacji krzyżowej są równoważne z AIC i BIC?

26

W pytaniu w innym miejscu na tej stronie, w kilku odpowiedziach wspomniano, że AIC jest równoważny walidacji krzyżowej z pominięciem jednego (LOO) i że BIC jest równoważny krzyżowej walidacji K-krotnie. Czy istnieje sposób empirycznego zademonstrowania tego w R, aby techniki zastosowane w LOO i K-fold zostały wyjaśnione i wykazano, że są równoważne wartościom AIC i BIC? Dobrze skomentowany kod byłby pomocny w tym względzie. Ponadto, demonstrując BIC, skorzystaj z pakietu lme4. Poniżej znajduje się przykładowy zestaw danych ...

library(lme4) #for the BIC function

generate.data <- function(seed)
{
    set.seed(seed) #Set a seed so the results are consistent (I hope)
    a <- rnorm(60) #predictor
    b <- rnorm(60) #predictor
    c <- rnorm(60) #predictor
    y <- rnorm(60)*3.5+a+b #the outcome is really a function of predictor a and b but not predictor c
    data <- data.frame(y,a,b,c) 
    return(data)    
}

data <- generate.data(76)
good.model <- lm(y ~ a+b,data=data)
bad.model <- lm(y ~ a+b+c,data=data)
AIC(good.model)
BIC(logLik(good.model))
AIC(bad.model)
BIC(logLik(bad.model))

We wcześniejszych komentarzach poniżej przedstawiłem listę nasion od 1 do 10000, w których AIC i BIC się nie zgadzają. Dokonano tego poprzez proste przeszukanie dostępnych nasion, ale jeśli ktoś mógłby zapewnić sposób generowania danych, który miałby tendencję do uzyskiwania rozbieżnych odpowiedzi na podstawie tych dwóch kryteriów informacyjnych, może to być szczególnie pouczające.

notable.seeds <- read.csv("http://student.ucr.edu/~rpier001/res.csv")$seed

Nawiasem mówiąc, zastanawiałem się nad zamówieniem tych nasion na podstawie tego, w jakim stopniu AIC i BIC się nie zgadzają, co próbowałem określić jako sumę absolutnych różnic AIC i BIC. Na przykład,

AICDiff <- AIC(bad.model) - AIC(good.model) 
BICDiff <- BIC(logLik(bad.model)) - BIC(logLik(good.model))
disagreement <- sum(abs(c(AICDiff,BICDiff)))

gdzie moje nieporozumienie ma sens tylko wtedy, gdy obserwacje są godne uwagi. Na przykład,

are.diff <- sum(sign(c(AICDiff,BICDiff)))
notable <- ifelse(are.diff == 0 & AICDiff != 0,TRUE,FALSE)

Jednak w przypadkach, w których AIC i BIC nie zgadzały się, obliczona wartość niezgodności była zawsze taka sama (i jest funkcją wielkości próbki). Patrząc wstecz na sposób obliczania AIC i BIC, widzę, dlaczego tak jest w przypadku obliczeń, ale nie jestem pewien, dlaczego tak jest w przypadku koncepcji. Jeśli ktoś mógłby to wyjaśnić, doceniłbym to.

russellpierce
źródło
+1 Kod byłby prosty do napisania, ale nadal bardzo interesuje mnie przejrzysty, ilustracyjny zestaw danych.
Nie jestem pewien, co wszystko powinno być w jasnym i ilustracyjnym zestawie danych, ale podjąłem próbę włączenia przykładowego zestawu danych.
russellpierce
Spójrz więc: to, co podałeś, jest przykładem bezużytecznego zestawu, ponieważ BIC i AIC dają te same wyniki: 340 v. 342 dla AIC i 349 v. 353 dla BIC - tak dobry. Model wygrywa w obu przypadkach. Cały pomysł z tą konwergencją polega na tym, że pewna walidacja krzyżowa wybierze ten sam model co odpowiadający jej układ scalony.
Zrobiłem proste skanowanie i na przykład dla seed 76, układy scalone nie zgadzają się.
1
Wow, obawiam się, że będzie to jeszcze trudniejsze. moim ogólnym punktem w całej dyskusji jest to, że zbieżność tych twierdzeń jest zbyt słaba, więc różnica może wynikać z przypadkowych fluktuacji. (I że nie działa przy uczeniu maszynowym, ale mam nadzieję, że to oczywiste.)

Odpowiedzi:

5

Próbując częściowo odpowiedzieć na moje pytanie, przeczytałem opis Wikipedii dotyczący weryfikacji krzyżowej z pominięciem jednego z nich

obejmuje wykorzystanie pojedynczej obserwacji z oryginalnej próbki jako danych walidacyjnych, a pozostałe obserwacje jako danych treningowych. Jest to powtarzane tak, że każda obserwacja w próbce jest wykorzystywana raz jako dane walidacyjne.

Podejrzewam, że w kodzie R oznaczałoby to coś takiego ...

resid <- rep(NA, Nobs) 
for (lcv in 1:Nobs)
    {
        data.loo <- data[-lcv,] #drop the data point that will be used for validation
        loo.model <- lm(y ~ a+b,data=data.loo) #construct a model without that data point
            resid[lcv] <- data[lcv,"y"] - (coef(loo.model)[1] + coef(loo.model)[2]*data[lcv,"a"]+coef(loo.model)[3]*data[lcv,"b"]) #compare the observed value to the value predicted by the loo model for each possible observation, and store that value
    }

... ma dawać wartości w pozostałościach, które są związane z AIC. W praktyce suma kwadratów reszt z każdej iteracji opisanej powyżej pętli LOO jest dobrym predyktorem AIC dla nasion godnych uwagi, r ^ 2 = .9776. Ale gdzie indziej autor sugerował, że LOO powinno być asymptotycznie równoważne z AIC (przynajmniej dla modeli liniowych), więc jestem trochę rozczarowany, że r ^ 2 nie jest bliżej 1. Oczywiście nie jest to tak naprawdę odpowiedź - bardziej jak dodatkowy kod, aby zachęcić kogoś do udzielenia lepszej odpowiedzi.

Dodatek: Ponieważ AIC i BIC dla modeli o stałej wielkości próby różnią się tylko o stałą, korelacja BIC z kwadratowymi resztami jest taka sama jak korelacja AIC z kwadratowymi resztami, więc podejście, które podjąłem powyżej, wydaje się bezowocne.

russellpierce
źródło
zauważ, że to będzie twoja zaakceptowana odpowiedź na nagrodę (w przypadku, gdy nie wybierzesz odpowiedzi, nagroda automatycznie wybierze odpowiedź z największą liczbą punktów)
robin girard
1
dobrze - przyznanie nagrody sobie wydaje się głupie - ale nikt inny nie udzielił odpowiedzi.
russellpierce