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;drbiglm現在の (0.9-1) バージョンでは、-class オブジェクトの AIC メソッド (より具体的には update メソッド) にかなり明白なバグがあるように見えますが、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)は と等しくないことに注意してください。これは、逸脱ではなく 2*対数尤度を使用するAIC(a)ためです(これら 2 つの尺度は加法係数が異なります) ...)stats:::AIC.default()

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)

biglmBen Bolker は、方程式が for を使用すると次のようになると指摘しましたAIC。

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

Ben はまたbiglm、残差 df の最初の値を更新していないと判断しました。

その新しい情報を考えると、2 つの方程式が同等であることがわかります。
まず、2 つの式を次のように制限します。ここが唯一の違いです。

(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

次に、私の方程式をベンの方程式に変形します。

私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であり、別の形式で表現されているようです。

もちろん、ベンが重要なコードとエラーの重要な説明の両方を提供してくれたからこそ、私はこれに気付くことができました。