Comment peut-on démontrer empiriquement dans R à quelles méthodes de validation croisée l'AIC et le BIC sont équivalents?

26

Dans une question ailleurs sur ce site, plusieurs réponses ont mentionné que l'AIC est équivalent à la validation croisée avec absence de contact (LOO) et que le BIC est équivalent à la validation croisée K-fold. Existe-t-il un moyen de démontrer empiriquement cela dans R de telle sorte que les techniques impliquées dans LOO et K-fold soient clairement définies et démontrées comme équivalentes aux valeurs AIC et BIC? Un code bien commenté serait utile à cet égard. De plus, pour démontrer le BIC, veuillez utiliser le package lme4. Voir ci-dessous pour un exemple de jeu de données ...

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

Selon les commentaires précédents, ci-dessous, j'ai fourni une liste de graines de 1 à 10000 dans lesquelles AIC et BIC ne sont pas d'accord. Cela a été fait par une simple recherche parmi les graines disponibles, mais si quelqu'un pouvait fournir un moyen de générer des données qui auraient tendance à produire des réponses divergentes à partir de ces deux critères d'information, cela pourrait être particulièrement instructif.

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

En passant, j'ai pensé à commander ces graines par la mesure dans laquelle l'AIC et le BIC ne sont pas d'accord, que j'ai essayé de quantifier comme la somme des différences absolues de l'AIC et du BIC. Par exemple,

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

où ma métrique de désaccord ne s'applique raisonnablement que lorsque les observations sont notables. Par exemple,

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

Cependant, dans les cas où AIC et BIC étaient en désaccord, la valeur de désaccord calculée était toujours la même (et est fonction de la taille de l'échantillon). En repensant à la façon dont AIC et BIC sont calculés, je peux voir pourquoi cela pourrait être le cas sur le plan informatique, mais je ne sais pas pourquoi ce serait conceptuellement le cas. Si quelqu'un pouvait également élucider ce problème, je l'apprécierais.

russellpierce
la source
+1 Le code serait simple à écrire, mais je suis toujours très intéressé à voir un ensemble de données clair et illustratif.
Je ne sais pas ce que tout devrait contenir dans un ensemble de données clair et illustratif, mais j'ai tenté d'inclure un exemple d'ensemble de données.
russellpierce
Alors regardez: ce que vous avez fourni est un exemple d'un ensemble inutile, car le BIC et l'AIC donnent les mêmes résultats: 340 v.342 pour AIC et 349 v.353 pour BIC - donc bon.model gagne dans les deux cas. L'idée avec cette convergence est que certaines validations croisées sélectionneront le même modèle que son IC correspondant.
J'ai fait un balayage simple et par exemple pour la graine 76, les CI ne sont pas d'accord.
1
Wow, c'est quelque chose qui sera encore plus difficile à obtenir, je le crains; mon point général dans toute la discussion est que la convergence de ces théorèmes est trop faible pour que la différence puisse émerger de fluctuations aléatoires. (Et que cela ne fonctionne pas pour l'apprentissage automatique, mais j'espère que cela est évident.)

Réponses:

5

Dans le but de répondre partiellement à ma propre question, j'ai lu la description de Wikipedia de la validation croisée avec absence de contact

consiste à utiliser une seule observation de l'échantillon d'origine comme données de validation et les observations restantes comme données d'apprentissage. Ceci est répété de telle sorte que chaque observation dans l'échantillon est utilisée une fois comme données de validation.

Dans le code R, je soupçonne que cela signifierait quelque chose comme ça ...

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
    }

... est censé donner des valeurs en resid qui sont liées à l'AIC. En pratique, la somme des résidus carrés de chaque itération de la boucle LOO détaillée ci-dessus est un bon prédicteur de l'AIC pour les graines notables, r ^ 2 = 0,9776. Mais, ailleurs, un contributeur a suggéré que LOO devrait être asymptotiquement équivalent à l'AIC (au moins pour les modèles linéaires), donc je suis un peu déçu que r ^ 2 ne soit pas plus proche de 1. Évidemment, ce n'est pas vraiment une réponse - plus comme du code supplémentaire pour essayer d'encourager quelqu'un à essayer de fournir une meilleure réponse.

Addendum: étant donné que l'AIC et le BIC pour les modèles de taille d'échantillon fixe ne varient que d'une constante, la corrélation du BIC aux résidus au carré est la même que la corrélation de l'AIC aux résidus au carré, donc l'approche que j'ai adoptée ci-dessus semble être infructueuse.

russellpierce
la source
notez que ce sera votre réponse acceptée pour la prime (au cas où vous ne choisiriez pas de réponse, la prime sélectionnera automatiquement la réponse avec le plus de points)
robin girard
1
bien - m'attribuer la prime semble idiot - mais personne d'autre n'a soumis de réponse.
russellpierce