AIC biglm和LM之间不同

Question

我一直在试图用biglm到大型数据集（约6000万线）运行线性回归。我想用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安全使用的同一数据集的比较biglm模型生成的工商行政管理机关？

Answer 1

文艺青年最爱的它看起来对我来说，没有在AIC方法biglm类对象一个很明显的错误（更具体地说，在更新方法），在目前（0.9-1）版本，但biglm包的作者是一个聪明，有经验的家伙，biglm被广泛使用，所以也许我失去了一些东西。谷歌搜索"biglm AIC df.resid"，似乎这已经discussed way back in 2009?更新：通过电子邮件，这确实是一个bug包的作者/维护者报告。

有趣的东西似乎是怎么回事。模型之间的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 *数似然，而不是越轨行为（这两项措施，在其添加剂系数不同）...）

Answer 2

我跟这个有点发挥各地。我不能肯定，但我认为由包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（其中+1是误差项）返回的估计系数的数量，并且deviance(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)

本Bolker指出，方程biglm用来AIC是：

deviance(object) + k * (object$n - object$df.resid)

本还确定，biglm未更新的残余DF的第一个值。

由于新的信息，我现在看到的是两个方程是等价的。 首先，限制两个方程以下，这是它们的不同的唯一的地方：

(n.parameters  + obs.added - 1) # mine

(object$n    - object$df.resid) # Ben's

重新安排矿占我加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使用公式。

当然，我才能够实现这一点，因为本提供两个关键代码和错误的关键解释。