r - biglm 和 lm 之间的 AIC 不同

Question

我一直在尝试使用 biglm 在大型数据集(约 60,000,000 行)上运行线性回归。我想使用 AIC 进行模型选择。但是我发现在较小的数据集上使用 biglm 时，biglm 返回的 AIC 变量与 lm 返回的变量不同。这甚至适用于 biglm 帮助中的示例。

data(trees)
ff<-log(Volume)~log(Girth)+log(Height)

chunk1<-trees[1:10,]
chunk2<-trees[11:20,]
chunk3<-trees[21:31,]

library(biglm)
a <- biglm(ff,chunk1)
a <- update(a,chunk2)
a <- update(a,chunk3)

AIC(a)#48.18546

a_lm <- lm(ff, trees)
AIC(a_lm)#-62.71125

有人可以解释一下这里发生了什么吗？使用 biglm 生成的 AIC 是否可以安全地用于比较同一数据集上的 biglm 模型？

Answer 1

tl;博士 在我看来，biglm 的 AIC 方法中存在一个非常明显的错误。 -class 对象(更具体地说，在更新方法中)，在当前(0.9-1)版本中，但作者biglm package 是一个聪明、有经验的人，而且 biglm被广泛使用，所以也许我错过了一些东西。谷歌搜索 "biglm AIC df.resid" ，看来这已经discussed way back in 2009? 更新 :包作者/维护者通过电子邮件报告这确实是一个错误。

这里似乎发生了一些有趣的事情。 差异 在模型之间的 AIC 中，跨建模框架应该是相同的，无论使用了哪些常量以及计算参数(因为这些常量和参数计数在建模框架内应该是一致的......)

原始示例:

data(trees)
ff <- log(Volume)~log(Girth)+log(Height)
chunk1<-trees[1:10,]
chunk2<-trees[11:20,]
chunk3<-trees[21:31,]
library(biglm)
a <- biglm(ff,chunk1)
a <- update(a,chunk2)
a <- update(a,chunk3)
a_lm <- lm(ff, trees)

现在拟合一个简化模型:

ff2 <- log(Volume)~log(Girth)    
a2 <- biglm(ff2, chunk1)
a2 <- update(a2, chunk2)
a2 <- update(a2 ,chunk3)
a2_lm <- lm(ff2,trees)

现在比较 AIC 值:

AIC(a)-AIC(a2)
## [1] 1.80222

AIC(a_lm)-AIC(a2_lm)
## [1] -20.50022

检查我们没有搞砸什么:

all.equal(coef(a),coef(a_lm))  ## TRUE
all.equal(coef(a2),coef(a2_lm))  ## TRUE

看看引擎盖下:

biglm:::AIC.biglm
## function (object, ..., k = 2) 
##    deviance(object) + k * (object$n - object$df.resid)

原则上这是正确的公式(观察数减去残差 df 应该是拟合的参数数)，但深入研究，它看起来像 $df.resid对象的组件未正确更新:

a$n  ## 31, correct
a$df.resid  ## 7, only valid before updating!

看着 biglm:::update.biglm ，我要补充

object$df.resid <- object$df.resid + NROW(mm)

在读取的行之前或之后

object$n <- object$n + NROW(mm)

...

这对我来说似乎是一个相当明显的错误，所以也许我遗漏了一些明显的东西，或者它已经被修复了。

一个简单的解决方法是将您自己的 AIC 函数定义为

AIC.biglm <- function (object, ..., k = 2) {
    deviance(object) + k * length(coef(object))
}

AIC(a)-AIC(a2)  ## matches results from lm()

(尽管请注意 AIC(a_lm) 仍然不等于 AIC(a) ，因为 stats:::AIC.default() 使用 2*log-likelihood 而不是偏差(这两个度量的附加系数不同)...)