biglm 和 lm 之间的 AIC 不同

我一直在尝试使用 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 模型？

tl;dr在我看来，AIC 方法中有一个非常明显的错误biglm-class 对象（更具体地说，在 update 方法中），在当前（0.9-1）版本中，但是作者biglm包是一个聪明、经验丰富的人，并且biglm被广泛使用，所以也许我遗漏了一些东西。谷歌搜索"biglm AIC df.resid"，看来这已经是早在2009年就讨论过了？ https://stat.ethz.ch/pipermail/r-help/attachments/20090708/e76206c7/attachment.pl Update：软件包作者/维护者通过电子邮件报告这确实是一个错误。

这里似乎发生了一些有趣的事情。这差异在 AIC 中，模型之间的模型在整个建模框架中应该是相同的，无论使用了什么常量以及如何计算参数（因为这些常量和参数计数应该一致）within建模框架...）

原始示例：

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 而不是偏差（这两个度量的加性系数不同）...）