2

我正在尝试根据 AIC 统计使用 R 进行模型选择。在比较带或不带加权的线性模型时，我在 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 估计要低得多。这种情况类似于在带有加权非线性回归 (nls) 的相关问题 AIC中提出的问题。

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

4

1 回答 1

1

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

在 的情况下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我相信in使用的公式的推导m2非常简单和正确，但我不太确定m4. logLik.nls()通过阅读有关（示例 1、示例 2 ）的其他一些线程，似乎对 nls 估计的正确方法存在一些混淆。总而言之，我认为AIC是正确的m2；我无法验证加权nls模型的数学，并且在这种情况下倾向于再次使用该公式（但用加权残差m2替换计算），或者（也许更好）不用于模型devianceAICnls

于 2018-01-13T03:39:26.847 回答