Wie kann man in R empirisch nachweisen, welchen Kreuzvalidierungsmethoden der AIC und der BIC gleichwertig sind?

26

In einer Frage an anderer Stelle auf dieser Website wurde in mehreren Antworten darauf hingewiesen, dass die AIC der LOO-Kreuzvalidierung und die BIC der K-fachen Kreuzvalidierung entspricht. Gibt es eine Möglichkeit, dies in R empirisch zu demonstrieren, sodass die mit LOO und K-fach verbundenen Techniken klargestellt werden und den AIC- und BIC-Werten entsprechen? Gut kommentierter Code wäre in dieser Hinsicht hilfreich. Verwenden Sie zur Demonstration des BIC außerdem das Paket lme4. Unten finden Sie einen Beispieldatensatz ...

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

In früheren Kommentaren habe ich unten eine Liste von Samen von 1 bis 10000 angegeben, in denen AIC und BIC nicht übereinstimmen. Dies geschah durch eine einfache Suche in den verfügbaren Samen, aber wenn jemand eine Möglichkeit zur Generierung von Daten bieten könnte, die dazu neigen, unterschiedliche Antworten auf diese beiden Informationskriterien zu liefern, könnte dies besonders informativ sein.

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

Nebenbei habe ich darüber nachgedacht, diese Samen nach dem Ausmaß zu ordnen, in dem AIC und BIC nicht übereinstimmen. Ich habe versucht, sie als Summe der absoluten Unterschiede zwischen AIC und BIC zu quantifizieren. Beispielsweise,

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

wo meine Meinungsverschiedenheit nur dann angemessen ist, wenn die Beobachtungen bemerkenswert sind. Beispielsweise,

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

In Fällen, in denen AIC und BIC nicht übereinstimmten, war der berechnete Nichtübereinstimmungswert immer der gleiche (und ist eine Funktion der Stichprobengröße). Wenn ich zurückblicke, wie AIC und BIC berechnet werden, kann ich sehen, warum dies rechnerisch der Fall sein könnte, aber ich bin nicht sicher, warum dies konzeptionell der Fall wäre. Wenn jemand auch dieses Problem klären könnte, würde ich es begrüßen.

russellpierce
quelle
+1 Der Code wäre einfach zu schreiben, trotzdem bin ich sehr daran interessiert, einen klaren, anschaulichen Datensatz zu sehen.
Ich bin nicht sicher, was alles in einem übersichtlichen und anschaulichen Datensatz enthalten sein muss, aber ich habe versucht, einen Beispieldatensatz aufzunehmen.
Russellpierce
Schauen Sie: Was Sie angegeben haben, ist ein Beispiel für einen unbrauchbaren Satz, da BIC und AIC die gleichen Ergebnisse liefern: 340 gegenüber 342 für AIC und 349 gegenüber 353 für BIC - also gewinnt good.model in beiden Fällen. Die ganze Idee mit dieser Konvergenz ist, dass eine bestimmte Kreuzvalidierung dasselbe Modell wie sein entsprechender IC auswählt.
Ich habe ein einfaches Scannen durchgeführt und zum Beispiel für Seed 76 stimmen die ICs nicht überein.
1
Wow, das ist etwas, das noch schwerer zu bekommen sein wird, fürchte ich; Mein allgemeiner Punkt in der gesamten Diskussion ist, dass die Konvergenz dieser Theoreme zu schwach ist, so dass der Unterschied aus zufälligen Fluktuationen hervorgehen kann. (Und dass es nicht für maschinelles Lernen funktioniert, aber ich hoffe, dass dies offensichtlich ist.)

Antworten:

5

In dem Versuch, meine eigene Frage teilweise zu beantworten, las ich die Beschreibung der einmaligen Kreuzvalidierung in Wikipedia

Dabei wird eine einzelne Beobachtung aus der Originalprobe als Validierungsdaten und die verbleibenden Beobachtungen als Trainingsdaten verwendet. Dies wird so wiederholt, dass jede Beobachtung in der Probe einmal als Validierungsdaten verwendet wird.

Im R-Code vermute ich, dass das so etwas bedeuten würde ...

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
    }

... soll Werte in Rückständen ergeben, die mit dem AIC zusammenhängen. In der Praxis ist die Summe der quadrierten Residuen aus jeder Iteration der LOO-Schleife, die oben beschrieben wurde, ein guter Prädiktor für den AIC für die notable.seeds, r ^ 2 = .9776. Aber an anderer Stelle schlug ein Mitwirkender vor, dass LOO asymptotisch der AIC entsprechen sollte (zumindest für lineare Modelle), daher bin ich ein wenig enttäuscht, dass r ^ 2 nicht näher bei 1 liegt. Offensichtlich ist dies keine wirkliche Antwort. eher wie zusätzlicher Code, um jemanden zu ermutigen, eine bessere Antwort zu geben.

Nachtrag: Da AIC und BIC für Modelle mit fester Stichprobengröße nur um eine Konstante variieren, ist die Korrelation von BIC zu quadratischen Residuen dieselbe wie die Korrelation von AIC zu quadratischen Residuen, so dass der Ansatz, den ich oben gewählt habe, erfolglos erscheint.

russellpierce
quelle
Beachten Sie, dass dies Ihre akzeptierte Antwort für das Kopfgeld ist (falls Sie keine Antwort wählen, wählt das Kopfgeld automatisch die Antwort mit den meisten Punkten)
Robin Girard
1
Gut - die Belohnung an mich selbst zu vergeben, scheint dumm - aber niemand anderes hat eine Antwort eingereicht.
Russellpierce