R中的AIC：使用加权数据时手动值与内部值的差异

Posted 2023-03-12

【中文标题】R中的AIC：使用加权数据时手动值与内部值的差异

【英文标题】：AIC in R: differences in manual vs. internal value when using weighted data

【发布时间】：2018-06-22 09:12:31

【问题描述】：

我正在尝试使用 R 进行基于 AIC 统计的模型选择。在比较带或不带加权的线性模型时，我在 R 中的代码告诉我，与不加权相比，加权更可取，并且这些结果在其他软件 (Graphpad Prism) 中得到了证实。我有使用来自标准曲线的真实数据的示例代码：

#Linear Curve Fitting
a <- c(0.137, 0.412, 1.23, 3.7, 11.1 ,33.3)
b <- c(0.00198, 0.00359, 0.00816, 0.0220, 0.0582, 0.184)
m1 <- lm(b ~ poly(a,1))
m2 <- lm(b ~ poly(a,1), weight=1/a)
n1 <- 6 #Number of observations
k1 <- 2 #Number of parameters

当我使用 R 中的内部函数或通过手动计算计算 AIC 时：

AIC = n + n log 2π + n log(RSS/n) + 2(k + 1) 有 n 个观测值和 k 参数

我得到了非加权模型的等效 AIC 值。当我分析加权的影响时，手动 AIC 值较低，但最终结果是内部和手动 AIC 都表明优先考虑加权。

> AIC(m1); n1+(n1*log(2*pi))+n1*(log(deviance(m1)/n1))+(2*(k1+1))
[1] -54.83171
[1] -54.83171
> AIC(m2); n1+(n1*log(2*pi))+n1*(log(deviance(m2)/n1))+(2*(k1+1))
[1] -64.57691
[1] -69.13025

当我尝试使用非线性模型进行相同的分析时，内部函数和手动计算之间的 AIC 差异更加深刻。以下是示例 Michaelis-Menten 动力学数据的代码：

c <- c(0.5, 1, 5, 10, 30, 100, 300)
d <- c(3, 5, 20, 50, 75, 200, 250)
m3 <- nls(d ~ (V * c)/(K + c), start=list(V=10, K=1))
m4 <- nls(d ~ (V * c)/(K + c), start=list(V=10, K=1), weight=1/d^2)
n2 <- 7
k2 <- 2

按照前两个模型的说明计算 AIC：

> AIC(m3); n2+(n2*log(2*pi))+n2*(log(deviance(m3)/n2))+(2*(k2+1))
[1] 58.48839
[1] 58.48839
> AIC(m4); n2+(n2*log(2*pi))+n2*(log(deviance(m4)/n2))+(2*(k2+1))
[1] 320.7105
[1] 0.1538546

与线性示例类似，当数据未加权 (m3) 时，内部 AIC 和手动 AIC 值相同。权重 (m4) 会出现问题，因为手动 AIC 估计要低得多。这种情况类似于在相关问题AIC with weighted nonlinear regression (nls) 中提出的问题。

我之前提到过 GraphPad Prism，对于上面给出的模型和数据集，当使用加权时，它的 AIC 较低。那么我的问题是，为什么在对数据进行加权时，R 中的内部与手动 AIC 估计存在如此差异（非线性模型与线性模型的结果不同）？最终，我应该认为内部 AIC 值还是手动值更正确，还是我使用了错误的公式？

【问题讨论】：

查看另一个问题中的 cmets，您的权重总和是否为 1？就“正确”而言，您链接的问题中有一个交叉验证的链接表明AIC is not a good metric for nonlinear models，因此统计上的正确答案可能是“都不正确”。

【参考方案1】：

您看到的差异是由于在加权模型的手动计算中使用了未加权对数似然公式。例如，您可以通过以下调整复制 m2 和 m4 的 AIC 结果：

对于m2，您只需从计算中减去sum(log(m2$weights))：

AIC(m2); n1+(n1*log(2*pi))+n1*(log(deviance(m2)/n1))+(2*(k1+1)) - sum(log(m2$weights))
[1] -64.57691
[1] -64.57691

对于m4，您必须将deviance 调用与加权残差计算交换，并从结果中减去n2 * sum(log(m4$weights))：

AIC(m4); n2+(n2*log(2*pi))+n2*(log(sum(m4$weights * m4$m$resid()^2)/n2))+(2*(k2+1)) - n2 * sum(log(m4$weights))
[1] 320.7105
[1] 320.7105

我相信logLikin m2 使用的公式的推导是非常直接和正确的，但我不确定m4。从阅读有关logLik.nls()（example 1，example 2）的其他一些线程来看，似乎对 nls 估计的正确方法有些混淆。总而言之，我相信AIC 对m2 是正确的；我无法验证加权 nls 模型的数学运算，并且在这种情况下倾向于再次使用 m2 公式（但用加权残差替换 deviance 计算），或者（可能更好）不使用 @987654342 @ 为nls 模型

【讨论】：

感谢丹尼尔森的回复！您对加权 AIC 方程的渲染与我正在研究的另一个模型的内部 AIC 匹配。