GLMのRのAIC / AICc / BIC式
Rがglm()
モデルオブジェクトの統計AIC、AICc(修正されたAIC)、およびBICを計算する方法を理解していることを確認しようとしています(revoScaleR::rxGlm()
オブジェクト、特にデフォルトでは使用できないAICcで同じ計算を実行できるようにするため))
私はこれらが次のように定義されていることを理解していました:
let p
=モデルパラメータの数
let n
=データポイントの数
AIC = deviance + 2p
AICc = AIC + (2p^2 + 2p)/(n-p-1)
BIC = deviance + 2p.log(n)
そこで、これらの番号を複製して、対応するR関数呼び出しと比較しようとしました。それは機能しませんでした:
library(AICcmodavg) # for the AICc() function
data(mtcars)
glm_a1 <- glm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb,
data = mtcars,
family = gaussian(link = "identity"),
trace = TRUE)
summary(glm_a1)
n <- nrow(glm_a1$data) # 32 p <- glm_a1$rank # 11
dev <- glm_a1$deviance# 147.49
my_AIC <- dev + 2 * p
my_AICc <- my_AIC + (2 * p^2 + 2 * p)/(n - p - 1)
my_BIC <- dev + 2 * p * log(n)
AIC(glm_a1) # 163.71
my_AIC # 169.49
AICc(glm_a1) # 180.13 (from AICcmodavg package)
my_AICc # 182.69
BIC(glm_a1) # 181.30
my_BIC # 223.74
を使用debug(AIC)
すると、計算が異なることがわかります。これは12個のパラメーターに基づいています(推定された分散/スケールパラメーター用に1つ余分ですか?)。また、対数尤度は、を使用しlogLik()
て取得されます-69.85
。これは、モデルの逸脱度が-2*-69.85 = 139.71
(そうではない)であることを私に示唆します。
誰かが私が間違ったことを知っていますか?ありがとうございました。
回答
でextractAIC
マニュアルページ
- Lは、適合の尤度およびedfと同等の自由度(つまり、通常のパラメトリックモデルのパラメーターの数)です。
- 一般化線形モデル(つまり、lm、aov、およびglmの場合)の場合、-2log Lは、deviance(fit)によって計算される逸脱度です。
- k = 2は従来のAICに対応し、k = log(n)を使用すると、代わりにBIC(ベイズIC)が提供されます。
したがって、
@ user20650のコメントと入力で次の議論を編集します
glm_a1$ranks
ガウス族で使用される近似分散を考慮せずに、近似パラメーターの数を返します。?glm
州逸脱度:最大の対数尤度の2倍を引いた定数まで。賢明な場合、飽和モデルの逸脱度がゼロになるように定数が選択されます。
それが理由です
-2*logLik(glm_a1) - deviance(glm_a1) = 7.78 > 0
summary(glm_a1)
次の行は、Dispersion parameter for gaussian family taken to be 7.023544
-2対数尤度と逸脱度のおよその差を返します。
library(AICcmodavg)
#> Warning: package 'AICcmodavg' was built under R version 3.6.2
#> Warning: no function found corresponding to methods exports from 'raster' for:
#> 'wkt'
data(mtcars)
glm_a1 <- glm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb,
data = mtcars,
family = gaussian(link = "identity"),
trace = TRUE)
#> Deviance = 147.4944 Iterations - 1
#> Deviance = 147.4944 Iterations - 2
(loglik <- logLik(glm_a1))
#> 'log Lik.' -69.85491 (df=12)
# thus the degrees of freedom r uses are 12 instead of 11
n <- attributes(loglik)$nobs # following user20650 recommendation
p <- attributes(loglik)$df # following user20650 recommendation
dev <- -2*as.numeric(loglik)
my_AIC <- dev + 2 * p
my_AICc <- my_AIC + (2 * p^2 + 2 * p)/(n - p - 1)
my_BIC <- dev + p * log(n)
BIC(glm_a1)
#> [1] 181.2986
my_BIC
#> [1] 181.2986
AIC(glm_a1)
#> [1] 163.7098
my_AIC
#> [1] 163.7098
AICc(glm_a1)
#> [1] 180.1309
my_AICc
#> [1] 180.1309
(逸脱度の「一定まで」の差を調整する)のrxGlm()
処理と一致するオブジェクトのこれらの量を計算する関数glm()
:
wrc_information_criteria <- function(rx_glm) # an object created by rxGlm()
{
# add 1 to parameter count for cases where the GLM scale parameter needs to be estimated (notably Gamma/gaussian)
extra_parameter_flag <- case_when(
rx_glm$family$family == "gaussian" ~ 1,
rx_glm$family$family == "Gamma" ~ 1,
rx_glm$family$family == "poisson" ~ 0,
rx_glm$family$family == "binomial" ~ 0,
TRUE ~ 999999999
)
n <- rx_glm$nValidObs p <- rx_glm$rank + extra_parameter_flag
dev <- rx_glm$deviance cat("\n") cat("n :", n, "\n") cat("p :", p, "\n") cat("deviance:", dev, "\n") AIC <- dev + 2 * p AICc <- AIC + (2 * p^2 + 2 * p)/(n - p - 1) BIC <- dev + p * log(n) # make a constant adjustment to AIC/AICc/BIC to give consistency with R's built in AIC/BIC functions applied to glm objects # can do this because rxGlm() supplies AIC already (consistent with R/glm()) - as long as computeAIC = TRUE in the function call deviance_constant_adjustment <- rx_glm$aic[1] - AIC
AIC <- AIC + deviance_constant_adjustment
AICc <- AICc + deviance_constant_adjustment
BIC <- BIC + deviance_constant_adjustment
cat("\n")
cat("AIC: ", AIC , "\n")
cat("AICc:", AICc, "\n")
cat("BIC: ", BIC , "\n")
}
それをテストしましょう...
data(mtcars)
glm_a1 <- glm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb,
data = mtcars,
family = gaussian(link = "identity"),
trace = TRUE)
glm_b1 <- rxGlm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb,
data = mtcars,
family = gaussian(link = "identity"),
verbose = 1,
computeAIC = TRUE)
AIC(glm_a1)
AICc(glm_a1)
BIC(glm_a1)
wrc_information_criteria(glm_b1) # gives same results for glm_b1 as I got for glm_a1