Hello Thierry,

Based on the code below, it looks like you do not need to worry about

likelihoods from lm() and gls(). They are defined the same way. I agree

with you that caution needs to be exercised. Simply because

mathematically the same likelihood may be defined using different constant.

Regards,

Have a good weekend

Andrzej Galecki

library(nlme)

data("sleepstudy", package = "lme4")

lmx <- lm(Reaction ~ Days, sleepstudy)

glsx<- gls(Reaction ~ Days, sleepstudy)

lmex <- lme(Reaction ~ Days, random= ~1|Subject, sleepstudy)

# LogLik values are the same

logLik(lmx) #

'log Lik.' -950.1465 (df=3)

logLik(glsx,REML=FALSE) # 'log Lik.'

-950.1465 (df=3)

# REML values are the same

logLik(lmx, REML=TRUE) # 'log Lik.'

-946.8318 (df=3)

logLik(glsx, REML=TRUE) # 'log Lik.'

-946.8318 (df=3)

# Two anovas give the same results

# (NOTE: Do not use anova for testing the need of random intercept!

Reference distribution other than chi-squre should be used)

anova(lmex, lmx)

anova(lmex, glsx)

# Model df AIC BIC logLik Test L.Ratio p-value

# lmex 1 4 1794.465 1807.192 -893.2325

# lmx 2 3 1899.664 1909.209 -946.8318 1 vs 2 107.1986 <.0001

On 3/18/2011 4:55 AM, ONKELINX, Thierry wrote:

> Dear Mark,

>

> I'm cc'ing this to the mixed models list to get some input from other experts. For them a link to the entire thread:

http://r.789695.n4.nabble.com/lmm-WITHOUT-random-factor-lme4-tp3384054p3384823.html>

> My comment was based on what I have read in Zuur et al. (2009).

>

> What worries me is that the loglikelihood of a lm() model and the equivalent gls() model is different. Although both models should be mathematically identical. Assuming that the loglikelihood is calculated on the same way within a package, I therefore have more confidence in comparing two models from the same package, thus gls() versus lme().

> Furthermore, I get an error when doing an anova between a lm() and a lme() model.

>

> Best regards,

>

> Thierry

>

>> library(nlme)

>> data("sleepstudy", package = "lme4")

>> fm<- lm(Reaction ~ Days, sleepstudy)

>> fm0<- gls(Reaction ~ Days, sleepstudy)

>> logLik(fm)

> 'log Lik.' -950.1465 (df=3)

>> logLik(fm0)

> 'log Lik.' -946.8318 (df=3)

>> fm1<- lme(Reaction ~ Days, random= ~1|Subject, sleepstudy)

>> anova(fm, fm1)

> Error in anova.lmlist(object, ...) :

> models were not all fitted to the same size of dataset

>> anova(fm0, fm1)

> Model df AIC BIC logLik Test L.Ratio p-value

> fm0 1 3 1899.664 1909.209 -946.8318

> fm1 2 4 1794.465 1807.192 -893.2325 1 vs 2 107.1986<.0001

>> sessionInfo()

> R version 2.12.2 (2011-02-25)

> Platform: i386-pc-mingw32/i386 (32-bit)

>

> locale:

> [1] LC_COLLATE=Dutch_Belgium.1252 LC_CTYPE=Dutch_Belgium.1252

> [3] LC_MONETARY=Dutch_Belgium.1252 LC_NUMERIC=C

> [5] LC_TIME=Dutch_Belgium.1252

>

> attached base packages:

> [1] stats graphics grDevices utils datasets methods base

>

> other attached packages:

> [1] nlme_3.1-98

>

> loaded via a namespace (and not attached):

> [1] grid_2.12.2 lattice_0.19-19 tools_2.12.2

>

> ----------------------------------------------------------------------------

> ir. Thierry Onkelinx

> Instituut voor natuur- en bosonderzoek

> team Biometrie& Kwaliteitszorg

> Gaverstraat 4

> 9500 Geraardsbergen

> Belgium

>

> Research Institute for Nature and Forest

> team Biometrics& Quality Assurance

> Gaverstraat 4

> 9500 Geraardsbergen

> Belgium

>

> tel. + 32 54/436 185

>

[hidden email]
> www.inbo.be

>

> To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.

> ~ Sir Ronald Aylmer Fisher

>

> The plural of anecdote is not data.

> ~ Roger Brinner

>

> The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.

> ~ John Tukey

>

>

> _______________________________________________

>

[hidden email] mailing list

>

https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>

>

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.