[R] polynomials REML and ML in nlme

[Berton Gunter]

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
>