[R] polynomials REML and ML in nlme

[Spencer Graves]

  Am I correct that changing the parameterization should NOT change the 
estimates of the variance components, because you are minimizing 
essentially the same objective function over the same subspace? The only 
thing that changes is the logdet(X'WX) term mentioned by I White? 
Moreover, letting X = QR, and using the fact that det(AB) = 
det(A)*det(B) if they are both square, we get det(R'Q'WQR) = 
det(Q'WQ)*det(R)^2. Thus, the change in parameterization affects only R, 
not Q, which means that it can't affect det(Q'WQ).

Is this accurate?

As a simple sanity check, I ran a very simple model with the same fixed 
effects in different parameterizations, and got the same estimates for 
the variance components under REEL but different 
"log-restricted-likelihood" (see below).

Comments?
spencer graves
############################
fit2 <- lme(y~i+I(i^2), random=~1|a, data=DF)
fit.p <- lme(y~poly(i, 2), random=~1|a, data=DF)

fit2
fit.p
Linear mixed-effects model
Fixed: y ~ i + I(i^2)
Data: DF
log-restricted-likelihood: -9.969998
Random effects:
Groups Name Variance Std.Dev.
a (Intercept) 1.4212 1.1921
Residual 1.9928 1.4117
# of obs: 6, groups: a, 2

Fixed effects:
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) -0.70889 2.67847 3 -0.2647 0.8084
i 0.19403 1.65425 3 0.1173 0.9140
I(i^2) -0.01081 0.23104 3 -0.0468 0.9656
 > fit.p
Linear mixed-effects model
Fixed: y ~ poly(i, 2)
Data: DF
log-restricted-likelihood: -6.728954
Random effects:
Groups Name Variance Std.Dev.
a (Intercept) 1.4212 1.1921
Residual 1.9928 1.4117
# of obs: 6, groups: a, 2

Fixed effects:
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) -0.193736 1.021136 3 -0.1897 0.8616
poly(i, 2)1 0.495131 1.454805 3 0.3403 0.7560
poly(i, 2)2 -0.066053 1.411677 3 -0.0468 0.9656
 >


I M S White wrote:

>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
>
>______________________________________________
>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
>  
>