lmm WITHOUT random factor (lme4)

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|

lmm WITHOUT random factor (lme4)

Baugh
LMM without Random effect:

I want to run an LMM both with and without the random factor (ID). And then extract the log-lik values from the two models in order to generate a p-value.

with random factor as:    lmer(y~x+(1|ID),data)

Question: can I simply substitute a dummy var (e.g. populated by zeros) for "ID" to run the model without the random factor? when I try this R returns values that seem reasonable, but I want to be sure this is appropriate.

From other inquires in the forum it seems that substituting ("ID-1") or some variant like that, does not work.

Thanks for your thoughts.
Baugh
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Mark Difford
On Mar 17, 2011; 11:43am Baugh wrote:

>> Question: can I simply substitute a dummy var (e.g. populated by zeros) for "ID" to run the model
>> without the random factor? when I try this R returns values that seem reasonable, but I want to be sure
>> this is appropriate.

If you can fit the model using lme (and it looks like you easily can) then another check would be:

## Compare models with and without random effects
fm <- lm(Reaction ~ Days, sleepstudy)
fm1 <- lme(Reaction ~ Days, random= ~1|Subject, sleepstudy)
anova(fm1, fm)                 ## lme-fitted model must come first

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Thierry Onkelinx
Dear Mark,

You cannot compare lm() with lme() because the likelihoods are not the same. Use gls() instead of lm()

library(nlme)
data("sleepstudy", package = "lme4")
fm <- lm(Reaction ~ Days, sleepstudy)
fm0 <- gls(Reaction ~ Days, sleepstudy)
logLik(fm)
logLik(fm0)

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

Best regards,

Thierry

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

> -----Oorspronkelijk bericht-----
> Van: [hidden email]
> [mailto:[hidden email]] Namens Mark Difford
> Verzonden: donderdag 17 maart 2011 10:55
> Aan: [hidden email]
> Onderwerp: Re: [R] lmm WITHOUT random factor (lme4)
>
> On Mar 17, 2011; 11:43am Baugh wrote:
>
> >> Question: can I simply substitute a dummy var (e.g. populated by
> >> zeros) for "ID" to run the model without the random factor? when I
> >> try this R returns values that seem reasonable, but I want
> to be sure
> >> this is appropriate.
>
> If you can fit the model using lme (and it looks like you
> easily can) then another check would be:
>
> ## Compare models with and without random effects fm <-
> lm(Reaction ~ Days, sleepstudy)
> fm1 <- lme(Reaction ~ Days, random= ~1|Subject, sleepstudy)
> anova(fm1, fm)                 ## lme-fitted model must come first
>
> Regards, Mark.
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/lmm-WITHOUT-random-factor-lme4-t
> p3384054p3384072.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Mark Difford
On Mar 17, 2011; 04:29pm Thierry Onkelinx wrote:

>> You cannot compare lm() with lme() because the likelihoods are not the same. Use gls() instead of lm()

Hi Thierry,

Of course, I stand subject to correction, but unless something dramatic has changed, you can. gls() can be used if you need to accommodate a correlation structure.

The method I have outlined, i.e. anova(lme$obj, lm$obj), is detailed in Pinheiro & Bates (2000) beginning page 154. Please refer to this if you doubt me.

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Mark Difford
On Mar 17, 2011; 04:29pm Thierry Onkelinx wrote:

>> You cannot compare lm() with lme() because the likelihoods are not the same. Use gls() instead of lm()

And perhaps I should have added the following:

First para on page 155 of Pinheiro & Bates (2000) states, "The anova method can be used to compare lme and lm objects."

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Thierry Onkelinx
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Mark Difford
This post has NOT been accepted by the mailing list yet.
On Mar 18, 2011; 10:55am Thierry Onkelinx wrote:

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

Hi Thierry,

You get this error because you have not done the comparison the way I said you should, by putting the lme$obj model first in the call to anova(). This ensures that lme's anova method is invoked, rather than lm's anova method.

You really do need to be careful with this.

## From my original posting
>> ## Compare models with and without random effects
>> fm <- lm(Reaction ~ Days, sleepstudy)
>> fm1 <- lme(Reaction ~ Days, random= ~1|Subject, sleepstudy)
>> anova(fm1, fm)                 ## lme-fitted model must come first

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Mark Difford
This post has NOT been accepted by the mailing list yet.
In reply to this post by Thierry Onkelinx
On Mar 18, 2011; 10:55am Thierry Onkelinx wrote:

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

Hi Thierry,

You get this error because you have not done the comparison the way I said you should, by putting the lme$obj model first in the call to anova(). This ensures that lme's anova method is invoked, rather than lm's anova method.

You really do need to be careful with this.

## From my original posting
>> ## Compare models with and without random effects
>> fm <- lm(Reaction ~ Days, sleepstudy)
>> fm1 <- lme(Reaction ~ Days, random= ~1|Subject, sleepstudy)
>> anova(fm1, fm)                 ## lme-fitted model must come first

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Mark Difford
In reply to this post by Thierry Onkelinx
On Mar 18, 2011; 10:55am Thierry Onkelinx wrote:

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

Hi Thierry,

You get this error because you have not done the comparison the way I said you should, by putting the lme$obj model first in the call to anova(). This ensures that lme's anova method is invoked, rather than lm's anova method.

You really do need to be careful with this...

## From my original posting
>> ## Compare models with and without random effects
>> fm <- lm(Reaction ~ Days, sleepstudy)
>> fm1 <- lme(Reaction ~ Days, random= ~1|Subject, sleepstudy)
>> anova(fm1, fm)                 ## lme-fitted model must come first

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa
Reply | Threaded
Open this post in threaded view
|

Re: lmm WITHOUT random factor (lme4)

Mark Difford
Apologies to all for the multiple posting. Don't know what caused it. Maybe it _is_ time to stop using Nabble after all...

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa
Reply | Threaded
Open this post in threaded view
|

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

agalecki
In reply to this post by Thierry Onkelinx
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

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

Mark Difford
On Mar 19, 2011; 01:39am Andrzej Galecki wrote:

>> I agree with you that caution needs to be exercised. Simply because mathematically the same
>> likelihood may be defined using different constant.

Yes. But this is ensured by the implementation. If the call to anova() is made with the lm$obj first in the sequence then an error is thrown. If the call is correctly made, with the lme$obj placed first in the sequence, then the log of the likelihood of each object is calculate by nlme:::logLik.lme using the same formula [via lapply(object, logLik, REML), where logLik points to nlme:::logLik.lme].

You will note, as Andrzej Galecki has pointed out, that the logLik.lm of the lm$obj is different from logLik.lme.

##
> logLik(fm)
'log Lik.' -950.1465 (df=3)
> logLik(fm, REML=T)
'log Lik.' -946.8318 (df=3)

Regards, Mark.
Mark Difford (Ph.D.)
Research Associate
Botany Department
Nelson Mandela Metropolitan University
Port Elizabeth, South Africa