Skip to content

Likelihood ratios

12 messages · Michael Lawrence, fengsj at mail.utexas.edu, Ben Bolker +2 more

Original
Michael Lawrence
# 
Hi folks,

I have 2 lmer fits, one (fit1) nested in the other (fit2), and I'd
like to compute the likelihood ratio comparing the models so I can say
something like "there is X times more evidence for fit1 than for fit2"
(as in Glover & Dixon, 2004, www.ncbi.nlm.nih.gov/pubmed/15732688).

I know I can use anova(fit1,fit2) to obtain a null-hypothesis
significance test of the fits, and I suspect the output also contains
the information I need to make my evidentiary statement, but I'm not
confident of what I'm doing here. Is it correct that the reported
value of chi-square from anova() is simply the D of the likelihood
ratio test (http://en.wikipedia.org/wiki/Likelihood_ratio_test)? If
so, does it sound right that I can simply derive the
complexity-corrected likelihood ratio as:

LR = exp( -2 * anova( fit1 , fit2 )$Chisq[2] )

?


Mike
Michael Lawrence
# 
oops, I guess that should be:

LR = exp( anova( fit1 , fit2 )$Chisq[2] / -2 )
On Tue, Jun 1, 2010 at 1:28 PM, Mike Lawrence <Mike.Lawrence at dal.ca> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Michael Lawrence
    

    

      
      #
    

  


  

      
After posting this, I thought to contact Pete Dixon himself and indeed
it seems he already coded the functions to obtain a likelihood ratio
comparing two lmer models:

AIC_lmer = function(x){
	require(lme4)
	print(formula(attr(x,"call")))
	summary(x)@AICtab
}

LR_lmer = function(m0,m1){
	exp((AIC_lmer(m0)[[1]]-AIC_lmer(m1)[[1]])/2)
}

#example usage:
LR_lmer( my_fit1 , my_fit2 )

        
On Tue, Jun 1, 2010 at 1:50 PM, Mike Lawrence <Mike.Lawrence at dal.ca> wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Michael Lawrence
    

    

      
      #
    

  


  

      
(sorry to spam!)

Pete Dixon also mentioned: "I think maybe you also have to do "REML=F"
in the lmer fits as well."

        
On Tue, Jun 1, 2010 at 3:08 PM, Mike Lawrence <Mike.Lawrence at dal.ca> wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      fengsj at mail.utexas.edu
    

    

      
      #
    

  


  

      
Is there someting similar for glmer models?
Thanks!


Quoting Mike Lawrence <Mike.Lawrence at dal.ca>:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
Yes, or

exp(logLik(fit2,REML=FALSE)-logLik(fit1,REML=FALSE))


   What is a "complexity-corrected" likelihood ratio?
   I'm very nervous about mixing in the AIC penalty terms with the rest
of the likelihood ratio -- and, as I think has been pointed out, you
have to be careful to set REML=FALSE ...

        
fengsj at mail.utexas.edu wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Michael Lawrence
    

    

      
      #
    

  


  

      
By "complexity-corrected", I simply meant a likelihood ratio that
takes into account potential differences in the number of parameters
between models, penalizing the more complex model. In the Glover &
Dixon paper I cited in my original post, they compute a likelihood
ratio based on ANOVA sums of squares, then apply a complexity
correction factor (they show formulae for both AIC and BIC) to yield a
final ratio that I've been calling a "complexity-corrected likelihood
ratio". Sorry for causing confusion by using my idiosyncratic
nomenclature!

        
On Tue, Jun 1, 2010 at 8:25 PM, Ben Bolker <bolker at ufl.edu> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      Shujuan Feng
    

    

      
      #
    

  


  

      
Hi all,



Is there a function similar to"lmList" for GLM?  I'd like to fit 
fixed-effect GLMs within each level of the random factors to take a look at 
the distribution of estimated parameters.



Thanks

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
Not quite as convenient as lmList, but something like

  L <- lapply(split(yourdata,yourdata$random_factor),glm,
              family=...)
  sapply(L,coef)

  should work fine.

        
Shujuan Feng wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      David Atkins
    

    

      
      #
    

  


  

      
Ben et al.--

Methinks lmList() in lme4 package takes a family argument:

Arguments

formula 	a linear formula object of the form y ~ x1+...+xn | g. In the 
formula object, y represents the response, x1,...,xn the covariates, and 
g the grouping factor specifying the partitioning of the data according 
to which different lm fits should be performed.

data 	a data frame in which to interpret the variables named in object.

family 	an optional family specification for a generalized linear model.

[snip]

So, I don't think a workaround is needed.

cheers, Dave

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      fengsj at mail.utexas.edu
    

    

      
      #
    

  


  

      
I use lmList to fit GLM (logit, with 3 contious predictors) within  
each leavel of my random factor. I can get results but with some  
warnings. I think the reason for the warnings is because GLM(logit)  
can not be used or estimated within some levels. (The dependents are  
all 0 within some levels. There are some NA for some estimated  
coefficients in the results).
Is there a way to remove those lm objects with NA values from the  
lmList results? I cannot use plot(),intervals() for the lmList  
results(I guess it is because of those NAs )
Thanks!




Quoting David Atkins <datkins at u.washington.edu>:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      fengsj at mail.utexas.edu
    

    

      
      #
    

  


  

      
I just read about:? No definition of residuals is completely  
satisfactory for binary data (in GLM)?.  Does anyone know how to check  
linearity for GLM(M) with the binary data?
Thanks


Quoting fengsj at mail.utexas.edu: