Как можно эмпирически продемонстрировать в R, каким методам перекрестной проверки AIC и BIC эквивалентны?

26

В вопросе, приведенном в другом месте на этом сайте, в нескольких ответах упоминалось, что AIC эквивалентна перекрестной проверке с пропуском (LOO) и что BIC эквивалентна перекрестной проверке в K-кратном размере. Есть ли способ эмпирически продемонстрировать это в R так, чтобы методы, задействованные в LOO и K-кратном прояснении, стали ясными и продемонстрировали, что они эквивалентны значениям AIC и BIC? Хорошо прокомментированный код будет полезен в этом отношении. Кроме того, при демонстрации BIC используйте пакет lme4. Ниже приведен пример набора данных ...

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

В предыдущих комментариях ниже я привел список семян от 1 до 10000, в которых AIC и BIC не согласны. Это было сделано с помощью простого поиска по доступным семенам, но если кто-то может предоставить способ генерирования данных, которые будут давать разные ответы из этих двух информационных критериев, это может быть особенно информативным.

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

Кроме того, я подумал о том, чтобы упорядочить эти семена по степени разногласий AIC и BIC, которые я попытался определить количественно как сумму абсолютных различий между AIC и BIC. Например,

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

где моя метрика несогласия применима только в разумных случаях, когда наблюдения заметны. Например,

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

Однако в случаях, когда AIC и BIC не соглашались друг с другом, рассчитанное значение разногласий всегда было одинаковым (и является функцией размера выборки). Оглядываясь назад на то, как рассчитываются AIC и BIC, я понимаю, почему это может иметь место в вычислительном отношении, но я не уверен, почему это было бы концептуально. Если бы кто-то мог объяснить и эту проблему, я был бы признателен.

russellpierce
источник
+1 Код было бы просто написать, но мне очень интересно увидеть четкий, иллюстративный набор данных.
Я не уверен, что все должно быть в четком и иллюстративном наборе данных, но я попытался включить примерный набор данных.
Расселпирс
Посмотрите: то, что вы предоставили, является примером бесполезного набора, потому что BIC и AIC дают одинаковые результаты: 340 против 342 для AIC и 349 против 353 для BIC - так что good.model выигрывает в обоих случаях. Вся идея этой конвергенции состоит в том, что определенная перекрестная проверка выберет ту же модель, что и соответствующая ей IC.
Я сделал простое сканирование, и, например, для семян 76 микросхемы не согласны.
1
Вау, это то, что будет еще сложнее получить, я боюсь; Мой общий смысл всей дискуссии состоит в том, что сходимость этих теорем слишком слабая, поэтому разница может возникнуть из-за случайных флуктуаций. (И это не работает для машинного обучения, но я надеюсь, что это очевидно.)

Ответы:

5

В попытке частично ответить на мой собственный вопрос, я прочитал описание Википедии о перекрестной проверке «оставь один раз»

включает в себя использование одного наблюдения из исходной выборки в качестве данных проверки, а остальные наблюдения - в качестве обучающих данных. Это повторяется так, что каждое наблюдение в выборке используется один раз в качестве данных проверки.

В коде R, я подозреваю, что это будет означать что-то вроде этого ...

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
    }

... должен давать значения в остатке, который связан с AIC. На практике сумма квадратов невязок от каждой итерации цикла LOO, описанного выше, является хорошим предиктором AIC для notable.seeds, r ^ 2 = .9776. Но в другом месте один из авторов предположил, что LOO должен быть асимптотически эквивалентен AIC (по крайней мере, для линейных моделей), поэтому я немного разочарован тем, что r ^ 2 не ближе к 1. Очевидно, это не совсем ответ - больше похоже на дополнительный код, чтобы побудить кого-то попытаться дать лучший ответ.

Приложение: Поскольку AIC и BIC для моделей с фиксированным размером выборки изменяются только на константу, корреляция BIC с квадратом остатков такая же, как и с AIC с квадратом остатков, поэтому подход, который я использовал выше, кажется бесплодным.

russellpierce
источник
обратите внимание, что это будет ваш принятый ответ за вознаграждение (в случае, если вы не выберете ответ, вознаграждение автоматически выберет ответ с наибольшим количеством очков)
Робин Жирар
1
хорошо - вручать награду мне кажется глупым - но никто другой не представил ответ.
russellpierce