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)

使用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 版本。

Roland

- Roland Sookias

---

我理解你的意思是想要从仅有10个观测值中，估计具有2个系数及另外3个ARMA过程系数的gls模型。 - Gavin Simpson

Answer 1

@GavinSimpson的上面评论是正确的，尝试从10个观测值中估计具有5个参数的模型，是非常令人希望的。通常的经验法则是，您应该至少有10倍于参数数量的数据点，而这只适用于标准固定效应/回归参数。（通常方差结构参数如自回归参数甚至比回归参数更难/需要更多数据来估计。）

话虽如此，在理想情况下，我们可以希望甚至从过度拟合的模型中估计参数。不过，我们来探究一下会发生什么：

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))进行此操作，则答案接近于（但不完全等于）使用corARMA(p=2)的gls结果。）

- Ben Bolker