[R] polynomials REML and ML in nlme

[I M S White]

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