R中nlme包中的gls函数出错答案

【问题标题】：Error with gls function in nlme package in RR中nlme包中的gls函数出错
【发布时间】：2011-07-15 11:09:50
【问题描述】：

我不断收到这样的错误：

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)

使用nlme 中的gls 函数。

前一个例子是模型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。

罗兰

【问题讨论】：

我是否正确地阅读了您想从 10 次观察中估计具有 2 个系数加上另外 3 个系数的 ARMA 过程的 gls 模型？

标签： r

【解决方案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) 的结果。）

【讨论】：