[R-sig-ME] Variance explained by random factor
Douglas Bates
bates at stat.wisc.edu
Tue Aug 26 15:35:31 CEST 2008
On Tue, Aug 26, 2008 at 12:29 AM, Ken Beath <ken at kjbeath.com.au> wrote:
> 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.
Yesterday I was looking at what is done in the glm function. Several
glm families (the non-"quasi" ones) have a member function to evaluate
the AIC then use the AIC value to determine the log-likelihood. I'm
not sure how that will extend to the generalized linear mixed model.
I think I will be able to reconstruct a log-likelihood with the
appropriate scaling factors for the binomial, poisson, gaussian and
Gamma families. I don't know about the inverse Gaussian and my guess
is that returning NA for the log-likelihood of the quasi families is a
sensible approach. Ben mentioned the QAIC definition which may be
reasonable.
I was considering these questions because we wanted to check the
results from the AGQ method developed by Bin Dai as a Google Summer of
Code project. In version 0.999375-26, which is now on CRAN, one can
specify nAGQ > 1 in glmer and nlmer but only for models with a single
grouping factor for the random effects.
>> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list