Duncan, I think that the problem in your comparison is locatted in the method of estimation of fixed-effects: in REML they are not estimated maximizing a joint likelihood function including fixed-effects and variance- covariance components as ML does. So when it comes down to different fixed-effects specifications, they can't be compared directly by means of the log-likelihoods.
Best Regards, Alex. > > -----Original Message----- > > From: [EMAIL PROTECTED] > > [mailto:[EMAIL PROTECTED] On Behalf Of dgoliche > > Sent: Tuesday, February 01, 2005 10:54 AM > > To: [email protected] > > Subject: [R] polynomials REML and ML in nlme > > > > Hello everyone, > > > > I hope this is a fair enough question, but I don�t have > > access to a copy > > of Bates and Pinheiro. It is probably quite obvious but the > > answer might > > be of general interest. > > If I fit a fixed effect with an added quadratic term and then do it as > > an orthogonal polynomial using maximum likelihood I get the expected > > result- they have the same logLik. > > > > mod2a<-lme(wthole~nplants+I(nplants^2),data=d3,random=~1|field > > /subplot,m > > ethod="ML") > > mod2b<-lme(wthole~poly(nplants,2),data=d3,random=~1|field/subp > > lot,method > > ="ML") > > > > > anova(mod2a,mod2b) > > Model df AIC BIC logLik > > mod2a 1 6 6698.231 6723.869 -3343.116 > > mod2b 2 6 6698.231 6723.869 -3343.116 > > > > However if I fit the two models by REML they are not considered to be > > the same and I get warned. > > > > > anova(mod2a.REML,mod2b.REML) > > Model df AIC BIC logLik > > mod2a.REML 1 6 6680.616 6706.219 -3334.308 > > mod2b.REML 2 6 6666.915 6692.518 -3327.457 > > > > Warning message: > > Fitted objects with different fixed effects. REML comparisons are not > > meaningful. in: anova.lme(mod2a.REML, mod2b.REML) > > > > Well yes, I suppose that�s right, they are not the same fixed effects. > > But why does REML give them such different Log likelihoods? And what > > should I do if I want to compare a larger set of models. For > > example the > > following, admittedly overparameterised model, can be fitted > > (slowly) by > > either method > > > > >mod4<-lme(wthole~poly(nplants,2),data=d3,random=~poly(nplants > > ,2)|field/ > > subplot,method="ML") > > > > But this doesn�t work by either method� > > > > > > > mod4<-lme(wthole~nplants+I(nplants^2),data=d3,random=~nplants+ > > I(nplants^ > > 2)|field/subplot,method="ML") > > > > Error in solve.default(pdMatrix(a, fact = TRUE)) : > > system is computationally singular: reciprocal > > condition number > > = 2.58558e-045 > > > > Thanks in advance for clearing up my confusion. The gentler the > > explanation, the more useful it would be as far as I am concerned as I > > am not a statistician and have to admit I am not at all clear how REML > > actually works. > > > > Duncan > > > > Dr Duncan Golicher > > Ecologia y Sistematica Terrestre > > Conservaci�n de la Biodiversidad > > El Colegio de la Frontera Sur > > San Cristobal de Las Casas, Chiapas, Mexico > > Tel. 967 1883 ext 9814 > > [EMAIL PROTECTED] > > > > > > [[alternative HTML version deleted]] > > > > ______________________________________________ > > [email protected] mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide! > > http://www.R-project.org/posting-guide.html > > > > ______________________________________________ > [email protected] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > __________________________________________________________________________ AntiPop-up UOL - � gr�tis! [[alternative HTML version deleted]] ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
