AIC/AICc/BIC GLM 的 R 公式

Question

我正在尝试检查我是否了解 R 如何计算 glm() 模型对象的统计 AIC、AICc（校正后的 AIC）和 BIC（以便我可以对 revoScaleR::rxGlm() 对象 - 特别是 AICc，默认情况下不可用）

我了解到这些定义如下：

让p =模型参数的数量

令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 个参数（一个额外的估计 dispersion/scale 参数？）。此外，使用 logLik() 获得对数似然，它返回一个数字 -69.85，这对我来说表明模型偏差将是 -2*-69.85 = 139.71（事实并非如此）。

有谁知道我做错了什么吗？谢谢。

Answer 1

在 extractAIC 手册中 page

其中：

L 是拟合的可能性，edf 是等效的自由度（即，通常参数模型的参数数量）。
对于广义线性模型（即，对于 lm、aov 和 glm），-2log L 是偏差，由 deviance(fit) 计算得出。
k = 2对应于传统的AIC，使用k = log(n)代替提供BIC（贝叶斯IC）。

因此

根据 @user20650

的评论和输入中的讨论进行编辑

glm_a1$ranks returns 不考虑高斯族中使用的拟合方差的拟合参数的数量。
?glm 状态

deviance: up to a constant, minus twice the maximized log-likelihood. Where sensible, the constant is chosen so that a saturated model has deviance zero.

这就是为什么 -2*logLik(glm_a1) - deviance(glm_a1) = 7.78 > 0
summary(glm_a1) returns 下一行 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

Answer 2

计算 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

AIC/AICc/BIC GLM 的 R 公式

AIC/AICc/BIC Formula in R for GLM

r

glm

根据 @user20650