r - 为什么 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

您可以检查 R 中的 summary(my_model) 和 python 中的 my_model.summary()，您会注意到，除了 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？您可以通过在 Rstudio 上键入 logLik 并按 Enter 来访问基础 R 中使用的 stats:::logLik.lm 函数。下面的两行似乎有问题:

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 :

For ‘"lm"’ fits it is assumed that the scale has been estimated (by maximum likelihood or REML)

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

对于它的值(value)，似乎也有很多来自 statsmodels 的讨论。侧面，例如this (closed) issue关于为什么 AIC/BIC 在 R 和 statsmodels 之间不一致......

This commit in March 2002显示 Martin Maechler 将“df”(自由度/模型参数数量)属性改回 object$rank+1带有以下附加注释:
帮助页面 ?logLik.lm yield :

Note that error variance \eqn{\sigma^2} is estimated in \code{lm()} and hence counted as well.

(此消息显然在稍后的某个时间点被编辑为上面看到的版本)。
新闻文件获得(在“错误修复”下):

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

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