我一直在尝试使用 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在我看来biglm，在当前 (0.9-1) 版本中，用于 -class 对象的 AIC 方法（更具体地说，在更新方法中）存在一个非常明显的错误，但该biglm包的作者是聪明，经验丰富的家伙，并且biglm被广泛使用，所以也许我错过了一些东西。谷歌搜索"biglm AIC df.resid"，这似乎早在 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 而不是偏差（这两个度量的相加系数不同）...）
我玩过这个。我不确定，但我认为AIC包使用的公式biglm是：
2 * (n.parameters + obs.added - 1) + deviance(a)
其中obs_added是 中的观察chunk2数加上 中的观察数chunk3：
obs.added <- dim(chunk2)[1] + dim(chunk3)[1]
并且是由（其中是误差项）n.parameters返回的估计系数的数量，并且是模型的偏差。summary(a) + 1+1deviance(a)a
####################################################
data(trees)
ff <- log(Volume)~log(Girth)+log(Height)
n.parm <- 4
chunk1<-trees[1:10,]
chunk2<-trees[11:20,]
chunk3<-trees[21:31,]
obs.added <- dim(chunk2)[1] + dim(chunk3)[1]
library(biglm)
a <- biglm(ff,chunk1)
a <- update(a,chunk2)
a <- update(a,chunk3)
AIC(a)
summary(a)
deviance(a)
2 * (n.parm + obs.added - 1) + deviance(a)
round(AIC(a), 5) == round(2 * (n.parm + obs.added - 1) + deviance(a), 5)
# [1] TRUE
####################################################
由于我不能 100% 确定我的答案是正确的，因此您可以使用下面的代码，看看您是否可以找到 AIC 的建议公式不起作用的场景。如果我发现任何此类情况，我将尝试根据需要修改下面的代码和上面的公式。
#########################################################
# Generate some data
n <- 118 # number of observations
B0 <- 2 # intercept
B1 <- -1.5 # slope 1
B2 <- 0.4 # slope 2
B3 <- 2.0 # slope 3
B4 <- -0.8 # slope 4
sigma2 <- 5 # residual variance
x1 <- round(runif(n, -5 , 5), digits = 3) # covariate 1
x2 <- round(runif(n, 10 , 20), digits = 3) # covariate 2
x3 <- round(runif(n, 2 , 8), digits = 3) # covariate 3
x4 <- round(runif(n, 12 , 15), digits = 3) # covariate 4
eps <- rnorm(n, mean = 0, sd = sqrt(sigma2)) # error
y <- B0 + B1 * x1 + B2 * x2 + B3 * x3 + B4 * x4 + eps # dependent variable
my.data <- data.frame(y, x1, x2, x3, x4)
# analyze data with linear regression
model.1 <- lm(my.data$y ~ my.data$x1 + my.data$x2 + my.data$x3 + my.data$x4)
summary(model.1)
AIC(model.1)
n.parms <- length(model.1$coefficients) + 1
my.AIC <- 2 * n.parms - 2 * as.numeric(logLik(model.1))
my.AIC
#########################################################
ff0 <- y ~ 1
ff1 <- y ~ x1
ff2 <- y ~ x1 + x2
ff3 <- y ~ x1 + x2 + x3
ff4 <- y ~ x1 + x2 + x3 + x4
n.parm0 <- 2
n.parm1 <- 3
n.parm2 <- 4
n.parm3 <- 5
n.parm4 <- 6
n.chunks <- 5
chunk1<-my.data[ 1:round(((nrow(my.data)/n.chunks)*1)+0),]
chunk2<-my.data[round(((nrow(my.data)/n.chunks)*1)+1):round(((nrow(my.data)/n.chunks)*2)+0),]
chunk3<-my.data[round(((nrow(my.data)/n.chunks)*2)+1):round(((nrow(my.data)/n.chunks)*3)+0),]
chunk4<-my.data[round(((nrow(my.data)/n.chunks)*3)+1):round(((nrow(my.data)/n.chunks)*4)+0),]
chunk5<-my.data[round(((nrow(my.data)/n.chunks)*4)+1):nrow(my.data),]
obs.added <- dim(chunk2)[1] + dim(chunk3)[1] + dim(chunk4)[1] + dim(chunk5)[1]
# check division of data
foo <- list()
foo[[1]] <- chunk1
foo[[2]] <- chunk2
foo[[3]] <- chunk3
foo[[4]] <- chunk4
foo[[5]] <- chunk5
all.data.foo <- do.call(rbind, foo)
all.equal(my.data, all.data.foo)
####################################################
library(biglm)
####################################################
a0 <- biglm(ff0, chunk1)
a0 <- update(a0, chunk2)
a0 <- update(a0, chunk3)
a0 <- update(a0, chunk4)
a0 <- update(a0, chunk5)
AIC(a0)
summary(a0)
deviance(a0)
print(a0)
2 * (n.parm0 + obs.added - 1) + deviance(a0)
round(AIC(a0), 5) == round(2 * (n.parm0 + obs.added - 1) + deviance(a0), 5)
####################################################
a1 <- biglm(ff1, chunk1)
a1 <- update(a1, chunk2)
a1 <- update(a1, chunk3)
a1 <- update(a1, chunk4)
a1 <- update(a1, chunk5)
AIC(a1)
summary(a1)
deviance(a1)
print(a1)
2 * (n.parm1 + obs.added - 1) + deviance(a1)
round(AIC(a1), 5) == round(2 * (n.parm1 + obs.added - 1) + deviance(a1), 5)
####################################################
a2 <- biglm(ff2, chunk1)
a2 <- update(a2, chunk2)
a2 <- update(a2, chunk3)
a2 <- update(a2, chunk4)
a2 <- update(a2, chunk5)
AIC(a2)
summary(a2)
deviance(a2)
print(a2)
2 * (n.parm2 + obs.added - 1) + deviance(a2)
round(AIC(a2), 5) == round(2 * (n.parm2 + obs.added - 1) + deviance(a2), 5)
####################################################
a3 <- biglm(ff3, chunk1)
a3 <- update(a3, chunk2)
a3 <- update(a3, chunk3)
a3 <- update(a3, chunk4)
a3 <- update(a3, chunk5)
AIC(a3)
summary(a3)
deviance(a3)
print(a3)
2 * (n.parm3 + obs.added - 1) + deviance(a3)
round(AIC(a3), 5) == round(2 * (n.parm3 + obs.added - 1) + deviance(a3), 5)
####################################################
a4 <- biglm(ff4, chunk1)
a4 <- update(a4, chunk2)
a4 <- update(a4, chunk3)
a4 <- update(a4, chunk4)
a4 <- update(a4, chunk5)
AIC(a4)
summary(a4)
deviance(a4)
print(a4)
2 * (n.parm4 + obs.added - 1) + deviance(a4)
round(AIC(a4), 5) == round(2 * (n.parm4 + obs.added - 1) + deviance(a4), 5)
####################################################
编辑
我建议biglm使用以下等式AIC：
2 * (n.parameters + obs.added - 1) + deviance(a)
Ben Bolker 指出，用于的等式biglm是AIC：
deviance(object) + k * (object$n - object$df.resid)
Ben 还确定biglm没有更新残差 df 的第一个值。
鉴于这些新信息，我现在看到这两个方程是等价的。
首先，将两个方程限制为以下，这是它们唯一不同的地方：
(n.parameters + obs.added - 1) # mine
(object$n - object$df.resid) # Ben's
重新排列我的，以说明我在参数数量上加 1，然后减 1：
((n.parameters-1) + obs.added) = ((4-1) + obs.added) = (3 + 21) = 24
现在把我的方程变成 Ben 的方程：
我3的是一样的：
(number of observations in chunk1 - object$df.resid) = (10 - 7) = 3
给予：
((number of obs in chunk1 - object$df.resid) + obs.added) = ((10-7) + 21)
或者：
(3 + 21) = 24
改编：
((number of obs in chunk1 + obs.added) - object$df.resid) = ((10 + 21) - 7)
或者：
(31 - 7) = 24
和：
((number of observations in chunk1 + obs.added) - object$df.resid)
是相同的：
(total number of observations - object$df.resid)
这与以下内容相同：
(object$n - object$df.resid) = (31 - 7) = 24
看来我提出的方程确实是biglm用于的方程AIC，只是以不同的形式表示。
当然，我之所以能够意识到这一点，是因为 Ben 提供了关键代码和错误的关键解释。