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/subplot,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 dgoliche@sclc.ecosur.mx [[alternative HTML version deleted]]
Duncan: Bates and Pinheiro do show explicitly and in detail that the restricted log likelihood depends on the parameterization: "An important difference between the likelihood function and the restricted likelihood is that the former is invariant to one-to-one reparametrizations of the fixed effects ... while the latter is not." (p.76) So that answers both your questions: With different fixed effect parameterizations you get different restricted log likelihoods even for the same random effects specifications. So you cannot compare models and you see differences of the sort you saw. That's just the way REML works. I believe the the basic issue can be summarized as: the restricted log likelihood is a nonlinear function of the fixed effects specification. Hence linear transformations (reparametrizations) of the fixed effects change the log likelihood. But if you want more than that, get the book (or another reference of your choice). You could also try googling "REML". All the above is subject to correction by Doug Bates, of course. -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA "The business of the statistician is to catalyze the scientific learning process." - George E. P. Box> -----Original Message----- > From: r-help-bounces at stat.math.ethz.ch > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of dgoliche > Sent: Tuesday, February 01, 2005 10:54 AM > To: r-help at stat.math.ethz.ch > 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 > dgoliche at sclc.ecosur.mx > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html >
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: r-help-bounces@stat.math.ethz.ch > > [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of dgoliche > > Sent: Tuesday, February 01, 2005 10:54 AM > > To: r-help@stat.math.ethz.ch > > 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 > > dgoliche@sclc.ecosur.mx > > > > > > [[alternative HTML version deleted]] > > > > ______________________________________________ > > R-help@stat.math.ethz.ch mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide! > > http://www.R-project.org/posting-guide.html > > > > ______________________________________________ > R-help@stat.math.ethz.ch 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]]
The REML loglikelihood includes a term -(1/2)logdet(X'WX) where X is the design matrix for the fixed effects and W is the inverse covariance matrix for the observations. Under reparametrisation, X becomes XM with M a non-singular matrix, and the REML loglikelihood changes by logdet(M). On Tue, 1 Feb 2005, dgoliche wrote:> 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/subplot,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 >=====================================I.White University of Edinburgh Ashworth Laboratories, West Mains Road Edinburgh EH9 3JT Tel: 0131 650 5490 Fax: 0131 650 6564 E-mail: iwhite at staffmail.ed.ac.uk