r - R 中 nlme 包中的 gls 函数出错

Question

我不断收到这样的错误:

Error in `coef<-.corARMA`(`*tmp*`, value = c(18.3113452983211, -1.56626248550284,  :
  Coefficient matrix not invertible

或者像这样:

Error in gls(archlogfl ~ co2, correlation = corARMA(p = 3)) : false convergence (8)

与 gls函数在 nlme .

前一个例子是模型 gls(archlogflfornma~nma,correlation=corARMA(p=3))哪里 archlogflfornma是

[1] 2.611840 2.618454 2.503317 2.305531 2.180464 2.185764 2.221760 2.211320

和 nma是

[1] 138 139 142 148 150 134 137 135

您可以在后者中看到模型，以及 archlogfl是

[1] 2.611840 2.618454 2.503317 2.305531 2.180464 2.185764 2.221760 2.211320
[9] 2.105556 2.176747

和 co2是

[1]  597.5778  917.9308 1101.0430  679.7803  886.5347  597.0668  873.4995 
[8]  816.3483 1427.0190  423.8917

我有 R 2.13.1。

罗兰

Answer 1

@GavinSimpson 上面的评论是，尝试从 10 个观察中估计具有 5 个参数的模型很有希望，这是正确的。一般的经验法则是，您的数据点数至少应是参数的 10 倍，这适用于标准的固定效应/回归参数。 (通常方差结构参数，如 AR 参数甚至比回归参数更难/需要更多的数据来估计。)

也就是说，在一个完美的世界中，人们甚至可以希望从过度拟合的模型中估计参数。让我们来探讨一下会发生什么:

archlogfl <- c(2.611840,2.618454,2.503317,
               2.305531,2.180464,2.185764,2.221760,2.211320,
               2.105556,2.176747)

co2 <- c(597.5778,917.9308,1101.0430,679.7803,
         886.5347,597.0668,873.4995,
         816.3483,1427.0190,423.8917)

看看数据，

plot(archlogfl~co2,type="b")
library(nlme)
g0 <- gls(archlogfl~co2)
plot(ACF(g0),alpha=0.05)

这是残差的自相关函数，具有 95% 的置信区间(请注意，这些是曲线置信区间，因此我们预计在任何情况下都会有大约 1/20 点落在这些边界之外)。

所以这里确实有一些(图形)证据表明存在一些自相关。我们将拟合具有详细输出的 AR(1) 模型(要了解估计这些参数的规模，您可能需要在 Pinheiro 和 Bates 2000 中深入研究:打印输出中显示的是参数，摘要中打印的是约束值......

g1 <- gls(archlogfl ~co2,correlation=corARMA(p=1),
    control=glsControl(msVerbose=TRUE))

让我们看看我们拟合了 AR1 后还剩下什么:

plot(ACF(g1,resType="normalized"),alpha=0.05)

现在适合 AR(2):

g2 <- gls(archlogfl ~co2,correlation=corARMA(p=2),
    control=glsControl(msVerbose=TRUE))

plot(ACF(g2,resType="normalized"),alpha=0.05)

正如您正确说明的那样，尝试转到 AR(3) 失败。

gls(archlogfl ~co2,correlation=corARMA(p=3))

您可以调整容差、起始条件等，但我认为这不会有太大帮助。

gls(archlogfl ~co2,correlation=corARMA(p=3,value=c(0.9,-0.5,0)),
    control=glsControl(tolerance=1e-4,msVerbose=TRUE),verbose=TRUE)

如果我绝对渴望获得这些值，我会编写自己的广义最小二乘函数，从头开始构建 AR(3) 相关矩阵，并尝试使用一些更强大的优化器运行它，但我真的必须有一个努力工作的好理由......

另一种选择是使用 arima以适应 gls 的残差或 lm无自相关拟合: arima(residuals(g0),c(3,0,0)) . (你可以看到，如果你用 arima(residuals(g0),c(2,0,0)) 这样做，答案接近(但不完全等于)来自 gls 和 corARMA(p=2) 的结果。)