Error with gls function in nlme package in R

Question

I keep getting an error like this:

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

or like this:

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

with the gls function in nlme.

The former example was with the model gls(archlogflfornma~nma,correlation=corARMA(p=3)) where archlogflfornma is

[1] 2.611840 2.618454 2.503317 2.305531 2.180464 2.185764 2.221760 2.211320

and nma is

[1] 138 139 142 148 150 134 137 135

You can see the model in the latter, and archlogfl is

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

and co2 is

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

I have R 2.13.1.

Roland

Answer 1

@GavinSimpson's comment above, that trying to estimate a model with 5 parameters from 10 observations is very hopeful, is correct. The general rule of thumb is that you should have at least 10 times as many data points as parameters, and that's for standard fixed effect/regression parameters. (Generally variance structure parameters such as AR parameters are even a bit harder/require a bit more data than regression parameters to estimate.)

That said, in a perfect world one could hope to estimate parameters even from overfitted models. Let's just explore what happens though:

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)

Take a look at the data,

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

This is an autocorrelation function of the residuals, with 95% confidence intervals (note that these are curvewise confidence intervals, so we would expect about 1/20 points to fall outside these boundaries in any case).

So there is indeed some (graphical) evidence for some autocorrelation here. We'll fit an AR(1) model, with verbose output (to understand the scale on which these parameters are estimated, you'll probably have to dig around in Pinheiro and Bates 2000: what's presented in the printout are the unconstrained values of the parameters, what's printed in the summaries are the constrained values ...

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

Let's see what's left after we fit AR1:

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

Now fit AR(2):

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

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

As you correctly state, trying to go to AR(3) fails.

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

You can play with tolerances, starting conditions, etc., but I don't think it's going to help much.

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

If I were absolutely desperate to get these values I would code my own generalized least-squares function, constructing the AR(3) correlation matrix from scratch, and try to run it with some slightly more robust optimizer, but I would really have to have a good reason to work that hard ...

Another alternative would be to use arima to fit to the residuals from a gls or lm fit without autocorrelation: arima(residuals(g0),c(3,0,0)). (You can see that if you do this with arima(residuals(g0),c(2,0,0)) the answers are close to (but not quite equal to) the results from gls with corARMA(p=2).)