[R-sig-ME] Variance explained by random factor

[brandon at brandoninvergo.com]

On Wed, 13 Aug 2008 16:17:48 +0100, "Renwick, A. R." <a.renwick at abdn.ac.uk>
wrote:
> 
> I am currently trying to run a lmer model with poisson distrubution.  I
> tested the model with a model without the random effect and it inferred
> that I should include the random effect:
> 
> ma1<-glm(RoundedOverlap~sess+breedfem,family=poisson,data=Male)
>
mixed<-lmer(RoundedOverlap~sess+breedfem+sess:breedfem+(1|Site),family=poisson,data=Male)
> 
> #test to see if sig difference between glm and glmm
> as.numeric(2*(logLik(mixed)-logLik(ma)))
> #99.16136
> pchisq(99.16136,1,lower=FALSE)
> #2.327441e-23  so should use a GLMM
> 
> However,the model output that I get states that the variance explained by
> the random factor is 0:
> 
> 
> Generalized linear mixed model fit by the Laplace approximation
> Formula: RoundedOverlap ~ sess + breedfem + sess:breedfem + (1 | Site)
>    Data: Male
>    AIC   BIC logLik deviance
>  109.9 127.2 -45.93    91.86
> Random effects:
>  Groups Name        Variance Std.Dev.
>  Site   (Intercept)  0        0
> Number of obs: 51, groups: Site, 14
> 
> I would really appreciate if somebody could help me understand why the
> variance is 0.
> Many thanks,
> Anna


I'm having a similar problem as Anna, however the responses in this thread
don't seem to apply to my problem. In my case, when I fit a GLMM to my data
(family = Gamma), the variance of the random effect is effectively zero
(1e-12). I also fit a GLM to the data and this model has practically the
same estimates and t-values for all the parameters. However, when I check
for significant differences between the two models via the log-likelihood
ratio test, the result is highly significant. I understand that a model
without the random effects term will be interpreted as a simpler model and
thus it will have a lower log-likelihood value but I don't understand how
the addition of a random effects term with such a small variance can cause
such a large increase in the log-likelihood of the model. Am I missing
something obvious here?

Thanks for your help!
-brandon


Here's my output:

> mixed.model <- lmer(dev.time ~ sex*temp + (1|sleeve.in.temp),
data=data.clean, family=Gamma(link="log"), method="ML")
> summary(mixed.model)
Generalized linear mixed model fit using Laplace 
Formula: dev.time ~ sex * temp + (1 | sleeve.in.temp) 
   Data: data.clean 
 Family: Gamma(log link)
   AIC   BIC logLik deviance
 29.28 87.32 -3.639    7.277
Random effects:
 Groups         Name        Variance   Std.Dev.  
 sleeve.in.temp (Intercept) 2.5683e-12 1.6026e-06
 Residual                   5.1366e-03 7.1670e-02
number of obs: 1446, groups: sleeve.in.temp, 54

Fixed effects:
             Estimate Std. Error t value
(Intercept)  3.659500   0.005993   610.6
sexm        -0.088450   0.008956    -9.9
temp21      -0.207295   0.008132   -25.5
temp23      -0.433156   0.008363   -51.8
temp25      -0.612696   0.008491   -72.2
temp27      -0.755628   0.008297   -91.1
sexm:temp21  0.008189   0.012000     0.7
sexm:temp23 -0.003357   0.012243    -0.3
sexm:temp25 -0.014887   0.012321    -1.2
sexm:temp27  0.010674   0.012405     0.9

Correlation of Fixed Effects:
            (Intr) sexm   temp21 temp23 temp25 temp27 sxm:21 sxm:23 sxm:25
sexm        -0.669                                                 
temp21      -0.737  0.493                                          
temp23      -0.717  0.480  0.528                                   
temp25      -0.706  0.472  0.520  0.506                            
temp27      -0.722  0.483  0.532  0.518  0.510                     
sexm:temp21  0.499 -0.746 -0.678 -0.358 -0.353 -0.361              
sexm:temp23  0.490 -0.731 -0.361 -0.683 -0.346 -0.354  0.546       
sexm:temp25  0.486 -0.727 -0.358 -0.349 -0.689 -0.351  0.542  0.532
sexm:temp27  0.483 -0.722 -0.356 -0.346 -0.341 -0.669  0.539  0.528 0.525


> glm.model <- glm(dev.time ~ sex*temp, data=data.clean,
family=Gamma(link="log"), na.action=na.omit)
> summary(glm.model)

Call:
glm(formula = dev.time ~ sex * temp, family = Gamma(link = "log"), 
    data = data.clean, na.action = na.omit)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.217280  -0.050703  -0.004718   0.043783   0.292775  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.659429   0.006014 608.474   <2e-16 ***
sexm        -0.088440   0.008987  -9.841   <2e-16 ***
temp21      -0.207203   0.008161 -25.391   <2e-16 ***
temp23      -0.433163   0.008392 -51.617   <2e-16 ***
temp25      -0.612562   0.008520 -71.895   <2e-16 ***
temp27      -0.755615   0.008326 -90.752   <2e-16 ***
sexm:temp21  0.008232   0.012041   0.684    0.494    
sexm:temp23 -0.003237   0.012285  -0.264    0.792    
sexm:temp25 -0.015077   0.012363  -1.220    0.223    
sexm:temp27  0.010813   0.012448   0.869    0.385    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1
‘ ’ 1 

(Dispersion parameter for Gamma family taken to be 0.005172238)

    Null deviance: 114.2564  on 1445  degrees of freedom
Residual deviance:   7.2771  on 1436  degrees of freedom
AIC: 5762.4

Number of Fisher Scoring iterations: 3
 

> as.numeric(-2*(logLik(glm.model)-logLik(mixed.model)))
[1] 5733.084
> pchisq(5733.084,1,lower=FALSE)
[1] 0

I checked the model's deviance as mentioned my Douglas Bates in this thread
and I got this:
> mixed.model at deviance
      ML     REML 
7.277151       NA