Como demonstrar empiricamente em R quais métodos de validação cruzada são equivalentes o AIC e o BIC?

26

Em uma pergunta em outro lugar neste site, várias respostas mencionaram que o AIC é equivalente à validação cruzada de não uso único (LOO) e que o BIC é equivalente à validação cruzada de dobras em K Existe uma maneira de demonstrar isso empiricamente em R, de modo que as técnicas envolvidas na LOO e na dobra K sejam esclarecidas e demonstradas como equivalentes aos valores da AIC e da BIC? Um código bem comentado seria útil nesse sentido. Além disso, na demonstração do BIC, use o pacote lme4. Veja abaixo um exemplo de conjunto de dados ...

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

Pelos comentários anteriores, abaixo, forneci uma lista de sementes de 1 a 10000, nas quais a AIC e a BIC discordam. Isso foi feito por uma pesquisa simples nas sementes disponíveis, mas se alguém pudesse fornecer uma maneira de gerar dados que tenderiam a produzir respostas divergentes a partir desses dois critérios de informação, isso seria particularmente informativo.

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

Como um aparte, pensei em encomendar essas sementes na medida em que a AIC e a BIC discordam, que tentei quantificar como a soma das diferenças absolutas da AIC e da BIC. Por exemplo,

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

onde minha métrica de desacordo se aplica apenas razoavelmente quando as observações são notáveis. Por exemplo,

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

No entanto, nos casos em que a AIC e a BIC discordaram, o valor calculado da discordância sempre foi o mesmo (e é uma função do tamanho da amostra). Olhando para trás, como o AIC e o BIC são calculados, posso ver por que isso pode ser o caso computacionalmente, mas não sei por que seria o caso conceitualmente. Se alguém puder esclarecer esse problema também, eu agradeceria.

russellpierce
fonte
+1 O código seria simples de escrever, ainda estou muito interessado em ver um conjunto de dados claro e ilustrativo.
Não tenho certeza do que todos precisariam em um conjunto de dados claro e ilustrativo, mas tentei incluir um conjunto de dados de amostra.
russellpierce
Então veja: o que você forneceu é um exemplo de um conjunto inútil, porque o BIC e o AIC fornecem os mesmos resultados: 340 v. 342 para o AIC e 349 v. 353 para o BIC - tão bom. O modelo vence em ambos os casos. A ideia toda com essa convergência é que determinada validação cruzada selecione o mesmo modelo que seu IC correspondente.
Fiz uma varredura simples e, por exemplo, para a semente 76, os CIs discordam.
11
Uau, isso é algo que será ainda mais difícil de obter, eu tenho medo; meu argumento geral em toda a discussão é que a convergência desses teoremas é muito fraca, de modo que a diferença possa emergir de flutuações aleatórias. (E que não está trabalhando para a aprendizagem de máquina, mas espero que isso é óbvio.)

Respostas:

5

Em uma tentativa de responder parcialmente à minha própria pergunta, li a descrição da Wikipedia sobre validação cruzada de deixar um fora

envolve o uso de uma única observação da amostra original como dados de validação e as demais observações como dados de treinamento. Isso é repetido de modo que cada observação na amostra seja usada uma vez como dados de validação.

No código R, suspeito que isso significaria algo assim ...

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
    }

... deve gerar valores em residentes relacionados ao AIC. Na prática, a soma dos resíduos quadráticos de cada iteração do loop LOO detalhada acima é um bom preditor da AIC para as sementes notáveis. R ^ 2 = 0,9776. Mas, em outro lugar , um colaborador sugeriu que a LOO deveria ser assintoticamente equivalente à AIC (pelo menos para modelos lineares), então estou um pouco decepcionado por r ^ 2 não estar mais perto de 1. Obviamente, isso não é realmente uma resposta - mais como código adicional para tentar incentivar alguém a tentar fornecer uma resposta melhor.

Adendo: Como o AIC e o BIC para modelos de tamanho fixo de amostra variam apenas uma constante, a correlação do BIC com os resíduos ao quadrado é a mesma que a correlação do AIC com os resíduos ao quadrado, portanto, a abordagem adotada acima parece infrutífera.

russellpierce
fonte
Note-se que esta será a sua resposta aceito pela graça (no caso de você não escolher uma resposta a recompensa selecionar automaticamente a resposta com mais pontos)
Girard robin
11
bem - conceder a recompensa a mim parece bobo - mas ninguém mais enviou uma resposta.
russellpierce