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

[Ken Beath]

On 26/08/2008, at 10:30 AM, <brandon at brandoninvergo.com> wrote:
>
> 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?
>

In addition to the problem with the missing constants in one of the  
likelihoods, there is also a scale parameter in the likelihood for  
Gamma so twice the LL difference is no longer distributed as a  
chisquare.

Ken


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