Obtaining GAMM intervals

Question

I'm using a GAMM and want to extract the intervals for the model with the Variance-Covariance Matrix but am having problems doing so

Coding:

fishdata <- read.csv("http://dl.dropbox.com/s/4w0utkqdhqribl4/fishdata.csv", 
header=T)
attach(fishdata)
require(mgcv)

gammdl <- gamm(inlandfao ~ s(marinefao), correlation = corAR1(form = ~1 | year), 
family=poisson())

summary(gammdl$gam)
intervals(gammdl$lme)

but the final line of coding returns,

Error in intervals.lme(gammdl$lme) : cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance

I don't see why this error message appears.

I'm trying to replicate what Simon Wood does on page 316 of Generalized Additive Models: An Introduction with R using my data.

Answer 1

This is really an issue of statistical computation, not R usage. Essentially, the error means that whilst the computations might have converged to a fitted model, there are problems with that model. What this boils down to is that the Hessian computed at the fitted values of the model is non-positive definite, and hence the maximum likelihood estimate of the covariance matrix is not available. Often this occurs where the log-likelihood function becomes flat, no further progress in optimisation can occur, and arises from the fact that the model is most likely too complex for the given data.

You could try fitting the model with and without the AR(1) and compare them using a generalised likelihood ratio test (via anova()). I would guess that the AR(1) nested within year is redundant - i.e. the parameter $\phi$ is effectively 0 - and it is this that is causing the error.

Actually, now it also occurs to me that you are trying to model the variable inlandfao as a smooth function of itself s(inlandfao). There is something wrong here - you'd expect a perfect fit without needing a smooth function. You need to fix that; you can't model the response as a function of itself.