[R] lme correlation structure error

[Douglas Bates]

Wayne Jones <JonesW at kssg.com> writes:

> I am trying to follow an example of modelling a serial correlation structure
> in the textbook "Mixed Effects Model in S and Splus". 
> However, I am getting some very odd results. Here is what I am trying to
> run: 
> 
> 
> library(nlme)
> data(Ovary)
> fm1<-lme(follicles~sin(2*pi*Time)+cos(2*pi*Time),data=Ovary,random=pdDiag(~s
> in(2*pi*Time)))
> 
> ### The example is fine up to here with all parameter estimates being
> identical to that in the book. 
> 
> 
> fm2<-update(fm1,correlation=corAR1())
> 
> #### The parameters of fm2 are different to that in the book and 
> 
> plot(ACF(fm2))
> 
> ##### signifies that serial correlation still exists in the residuals.
> 
> 
> 
> 
> ###### When I try and run this (which runs fine in the book)
> 
> fm5<-update(fm1,corr=corARMA(p=1,q=1))
> 
> #### I get the following error message 
> 
> #Error in "coef<-.corARMA"(*tmp*, value = c(62.3428455941166,
> 62.3428517930051 : 
> #       Coefficient matrix not invertible
> 
> 
> I have tried running the example on R for windows versions 1.7.1 and 1.8.1
> with the same results. 
> 
> Can anyone shed light on this?? Perhaps another R-user could kindly run this
> example and see if they get the same results??

I believe it is the optimizer being used by lme that causes the
problem.  In the initial iterations the nlm function will sometimes
take very large steps in the parameter space and end up stuck in
places where the likelihood is very flat.  Use 

control = list(msVerbose = TRUE)

in the call to lme to see where the parameter vector is being
directed.

The examples in our book were done in S-PLUS (version 3.4, I think)
which uses different optimizer code, and that optimizer code does not
appear to be freely available.  (It's from the PORT library from Bell
Labs.  It can be downloaded but you have to click through an agreement
panel and give your name, email address, etc.)