[R-sig-ME] [R] lmm WITHOUT random factor (lme4)

[Andrzej Galecki]

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
> Thierry.Onkelinx at inbo.be
> 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
>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>