为什么 R 产生不正确的 AIC 和 BIC

Question

我已经用谷歌搜索了这个并找不到解决方案。

R 似乎在 AIC/BIC 计算方面存在问题。它会产生错误的结果。一个简单的例子如下所示：

link = 'https://gist.githubusercontent.com/seankross/a412dfbd88b3db70b74b/raw/5f23f993cd87c283ce766e7ac6b329ee7cc2e1d1/mtcars.csv'
df = read.csv(link, row.names = 'model')
form = 'mpg ~ disp + hp + wt + qsec + gear'
my_model = lm(form, data = df)
summary(my_model)
cat('AIC:',AIC(my_model),'\tBIC:',AIC(my_model, k = log(nrow(df))))

AIC: 157.4512   BIC: 167.7113

在 python 中做完全相同的事情，我得到：

import pandas as pd
from statsmodels.formula.api import ols as lm

link = 'https://gist.githubusercontent.com/seankross/a412dfbd88b3db70b74b/raw/5f23f993cd87c283ce766e7ac6b329ee7cc2e1d1/mtcars.csv'
df = pd.read_csv(link, index_col='model')
form = 'mpg ~ disp + hp + wt + qsec + gear'
my_model = lm(form, df).fit()
my_model.summary()
print(f'AIC: {my_model.aic:.4f}\tBIC: {my_model.bic:.4f}')

AIC: 155.4512   BIC: 164.2456

您可以summary(my_model)在 R 和my_model.summary()python 中检查，您会注意到这两个模型在所有方面都完全相同，除了 AIC 和 BIC。

我决定在 R 中手动计算它：

p = length(coef(my_model)) # number of predictors INCLUDING the Intercept ie 6
s = sqrt(sum(resid(my_model)^2)/nrow(df)) #sqrt(sigma(my_model)^2 * (nrow(df) - p)/nrow(df))
logl = -2* sum(dnorm(df$mpg, fitted(my_model),s, log = TRUE))

c(aic = logl + 2*p, bic = logl + log(nrow(df))*p)
 aic      bic 
155.4512 164.2456 

这与python产生的结果相匹配。

深入挖掘，我注意到 AIC 确实使用了该logLik功能。这就是问题出现的地方：在乘以之前给出logLik(my_model)与logl上面显示的完全相同的结果，-2但df给出的是 7 而不是 6。

如果我强制排名以使其成为 6，我会得到正确的结果，即：

my_model$rank = my_model$rank - 1
cat('AIC:',AIC(my_model),'\tBIC:',AIC(my_model, k = log(nrow(df))))

AIC: 155.4512   BIC: 164.2456

为什么 R 将预测变量的数量加 1？您可以logLik通过stats:::logLik.lm在 Rstudio 上键入并按 Enter来访问在基础 R 中使用的函数。下面的两行似乎有问题：

function (object, REML = FALSE, ...) 
{
    ...
    p <- object$rank
    ...
    attr(val, "df") <- p + 1 # This line here. Why does R ADD 1?
    ...
}

Answer 1

这显然是一个深思熟虑的选择：R 计算估计参数集中的尺度参数。来自?logLik.lm：

对于 '"lm"' 拟合，假设已经估计了比例（通过最大似然或 REML）

（另见此处，@MrFlick 在评论中指出）。这种歧义（以及，归一化常数是否包含在对数似然中：在 R 中，它们是）在跨平台比较结果之前总是必须检查，有时甚至跨同一平台内的过程或函数。

对于它的价值，似乎也有很多关于这一点的讨论statsmodels，例如这个（封闭的）问题关于为什么 R 和 statsmodels 之间的 AIC/BIC 不一致......

2002 年 3 月的这个提交显示 Martin Maechler 将“df”（自由度/模型参数数量）属性改回object$rank+1了以下附加注释：

帮助页面?logLik.lm获得：

请注意，误差方差 \eqn{\sigma^2} 是在 \code{lm()} 中估计的，因此也被计算在内。

（此消息显然在稍后的某个时间点被编辑为上面看到的版本）。

新闻文件获得（在“BUG FIXES”下）：

 o    logLik.lm() now uses  "df = p + 1" again (`+ sigma'!).

我很难再做比这更远的考古学（即大概是基于这里的消息p+1，最初使用了推算，然后有人改用了它，p而 MM 在 2002 年改回了它），因为函数四处移动（这个文件于 2001 年创建，因此较早的版本将更难找到）。我在2002 年 2 月或 3 月的r-devel 邮件列表存档中没有找到任何关于此的讨论......