[R] nlme: Computing REML likelihood value from ML likelihood value

[CHAOUCH Aziz]

Dear All,

I have a question concerning the computation of the value of the Restricted Maximum Likelihood (REML) function evaluated at a given set of parameter estimates from the Maximum likelihood (ML) value. Following the book of Fitzmaurice, Laird and Ware (2004) "Applied Longitudinal Analysis" pp101, the REML likelihood can be computed by multiplying the ML likleihood by the square root of the generalized variance of the fixed effects. When working on the log scale, this means:

log.REML = 0.5*(log(det(W))) + log.ML

where W is the covariance matrix of the fixed effects estimates.

However when using the function lme in the nlme package, the REML likelihood value does not agree with this definition. As an example consider the following code:

library(nlme)
data(Glucose)
model <- lme(conc~Time, data=Glucose, random=~Time|Subject, na.action=na.omit)
LL.ml <- logLik(model,REML=F) #=688.588
LL.reml <- logLik(model,REML=T) #=692.654
LL.reml2 <- as.numeric(0.5*log(det(model$varFix)) + LL.ml) #=694.489

Here LL.reml2 is computed according to the above formula and it does not agree with the REML value as computed in lme. Does anybody know why is this so and how I may evaluate log.REML at a set of parameter values using log.ML?

Regards,

Aziz Chaouch