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: [EMAIL PROTECTED] ______________________________________________ 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
