Wayne Jones <[EMAIL PROTECTED]> 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.) ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html