如何从 lm_robust 对象中获取 AIC答案

【问题标题】：How to get AIC from lm_robust object如何从 lm_robust 对象中获取 AIC
【发布时间】：2020-01-01 22:13:28
【问题描述】：

如何从lm_robust 对象（包estimatr）中获取AIC？我使用lm_robust 是因为我想使用稳健的估计器来计算SE。与lm 函数不同，当您运行摘要函数时不提供AIC，并且在lm_robust 对象上运行AIC 函数会产生错误。下面是我尝试运行的那种模型的玩具示例。

library(estimatr)

fake_data<-data.frame(outcome=rnorm(100,3.65,1),
                      pred1=rnorm(100,15,7),
                      pred2=as.factor(sample(1:5, 100, replace = T)))

mod1<-lm_robust(outcome~pred1+pred2,data=fake_data)
AIC(mod1)

错误信息如下所示：

> AIC(mod1)
Error in UseMethod("logLik") : 
  no applicable method for 'logLik' applied to an object of class "lm_robust"

【问题讨论】：

我不是 100% 相信稳健回归甚至具有明确定义的对数似然（即，可以在计算似然比测试中按原样使用的数字, AIC 值, 等等等等) ... ???下面的解决方案给出了数字，但我担心它们可能不是有意义的数字。例如见stats.stackexchange.com/questions/180451/…（未答复）
@BenBolker 实际上，即使我给出了答案，我也同意你的看法：)。因为，AIC 和 BIC 标准主要用于时间序列分析。我在lm 中没有看到任何对AIC 感兴趣的人。（如果我错了，如果在学术界存在使用AIC 的人，请纠正我）。底线，如果问题是计算方法，下面是正确的。但是，如果问题是lm_robust 中的 AIC 的必要性...我不太确定。
生态学和进化领域的很多人（我的研究领域）将 AIC 用于线性模型（尽管对于 GLM 可能稍微更频繁一些）
@BenBolker 哇...我学到了一些东西。谢谢你。实际上，我主要使用 AIC 来选择时间序列中的滞后顺序等来确定更适合的模型。当谈到线性模型之间的比较时，因为它很常见，我一直认为应该看看 R-square。那么是什么让 AIC 成为线性模型所必需的，而不是考虑 R 方？
AIC 将非常类似于 adjusted R2 或 Mallows 的 Cp 统计量，即它考虑了模型的复杂性。 This CV answer 建议 AIC 具有更好的模型选择属性……“??????2????????????中的惩罚足以使??????2?????? ?????? 人口的无偏估计量 ??????2 当没有回归量实际上属于模型时......但是，??????2???????????? 是不是最佳模型选择器。”

标签： r regression lm

【解决方案1】：

这里有一个解决方法

mod1 = lm_robust(outcome ~ pred1 + pred2, data = fake_data)

#Create any fitted model using 'lm' as a placeholder
mod2 = with(list(x = rnorm(10), y = rnorm(10)), lm(y ~ x))

#Copy values in `mod2` from `mod1`
mod2[names(mod2)] = mod1[names(mod2)]

#Calculate residuals in `mod2`
mod2$residuals = mod2$fitted.values - fake_data$outcome

AIC(mod2)
#[1] 326.6092

【讨论】：

但是有没有办法为 lm_robust 获得 AIC？对于此分析，我想使用稳健的标准误差估计器，这就是我使用 lm_robust 的原因

【解决方案2】：

如果一定要lm_robust，可以选择自己计算如下，

AIC的公式，

AIC = 2*k + n [Ln( 2(pi) RSS/n ) + 1]

# n : Number of observation
# k : All variables including all distinct factors and constant
# RSS : Residual Sum of Square

如果我们针对您的情况将其应用于R，

# Note that, I take k=7 since you have, 5 factors + 1 continuous and 1 constant

AIC_calculated <- 2*7 + 100* (log( 2*pi* (1-mod1$r.squared)*mod1$tss/100 ) + 1)


[1] 332.2865

这与lm 和glm 输出相同。

mod2<-lm(outcome~pred1+pred2,data=fake_data)

> AIC(mod2)
[1] 332.2865

最后，当然，您可以将这个计算放入一个函数中，以便随时调用，只需在其中提供lm_robust 模型，而无需为任何给定数据设置N 和k 参数，例如，

myAIC <- function(data) {

    2*(data$k+1) + data$N * (log(2*pi* (1-data$r.squared)*data$tss/data$N ) + 1)

}

> myAIC(mod1)
[1] 332.2865

注意：由于在数据帧中运行sample() 函数时的播种差异，结果可能在您的计算机中显示不同。

【讨论】：