R 中的 AIC：使用加权数据时手动值与内部值的差异
AIC in R: differences in manual vs. internal value when using weighted data
我正在尝试使用 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 差异更为明显。以下是米氏动力学数据的示例代码：
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 值相同。由于手动 AIC 估计要低得多，加权 (m4) 会出现问题。这种情况类似于相关问题 AIC with weighted nonlinear regression (nls) 中提出的问题。
我之前提到过 GraphPad Prism，对于上面给出的模型和数据集，当使用加权时，它显示出较低的 AIC。那么我的问题是，为什么在对数据进行加权时，R 中的内部 AIC 估计与手动 AIC 估计存在如此差异（非线性模型与线性模型的结果不同）？最终，我应该认为内部 AIC 值或手动值更正确，还是我使用了错误的公式？
您看到的差异是由于在加权模型的手动计算中使用了未加权的对数似然公式。例如，您可以通过以下调整复制 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
我相信 logLik 在 m2 中使用的公式的推导是非常直接和正确的，但我不确定 m4。通过阅读其他一些关于 logLik.nls() (example 1, example 2) 的帖子，似乎对 nls 估计的正确方法存在一些困惑。总而言之，我认为 AIC 对 m2 是正确的；我无法验证加权 nls 模型的数学运算，在这种情况下我会倾向于再次使用 m2 公式（但用加权残差替换 deviance 计算），或者（也许更好）不要将 AIC 用于 nls 模型
我正在尝试使用 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 差异更为明显。以下是米氏动力学数据的示例代码：
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 值相同。由于手动 AIC 估计要低得多，加权 (m4) 会出现问题。这种情况类似于相关问题 AIC with weighted nonlinear regression (nls) 中提出的问题。
我之前提到过 GraphPad Prism，对于上面给出的模型和数据集，当使用加权时，它显示出较低的 AIC。那么我的问题是，为什么在对数据进行加权时，R 中的内部 AIC 估计与手动 AIC 估计存在如此差异（非线性模型与线性模型的结果不同）？最终，我应该认为内部 AIC 值或手动值更正确，还是我使用了错误的公式？
您看到的差异是由于在加权模型的手动计算中使用了未加权的对数似然公式。例如，您可以通过以下调整复制 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
我相信 logLik 在 m2 中使用的公式的推导是非常直接和正确的，但我不确定 m4。通过阅读其他一些关于 logLik.nls() (example 1, example 2) 的帖子，似乎对 nls 估计的正确方法存在一些困惑。总而言之，我认为 AIC 对 m2 是正确的；我无法验证加权 nls 模型的数学运算，在这种情况下我会倾向于再次使用 m2 公式（但用加权残差替换 deviance 计算），或者（也许更好）不要将 AIC 用于 nls 模型