[R] Fitting binomial lmer-model, high deviance and low logLik
Ken Beath
kbeath at efs.mq.edu.au
Thu Dec 15 12:44:54 CET 2005
Try using method="AGQ" to use the adaptive Gaussian quadrature
method. This will generally give a more accurate result than PQL. If
this doesn't give a more meaningful result, then it may be your data.
Within each mother are the outcomes all identical ? This will give
the random effects model a lot of problems.
Ken
> From: Ivar Herfindal <ivar.herfindal at bio.ntnu.no>
> Subject: [R] Fitting binomial lmer-model, high deviance and low logLik
> To: r-help at stat.math.ethz.ch
> Message-ID: <439FF524.1040403 at bio.ntnu.no>
> Content-Type: text/plain; charset=us-ascii; format=flowed
>
> Hello
>
> I have a problem when fitting a mixed generalised linear model with
> the
> lmer-function in the Matrix package, version 0.98-7. I have a respons
> variable (sfox) that is 1 or 0, whether a roe deer fawn is killed
> or not
> by red fox. This is expected to be related to e.g. the density of red
> fox (roefoxratio) or other variables. In addition, we account for
> family
> effects by adding the mother (fam) of the fawns as random factor. I
> want
> to use AIC to select the best model (if no other model selection
> criterias are suggested).
>
> the syntax looks like this:
>> mod <- lmer(sfox ~ roefoxratio + (1|fam), data=manu2,
>> family=binomial)
>
> The output looks ok, except that the deviance is extremely high
> (1.798e+308).
>
>> mod
> Generalized linear mixed model fit using PQL
> Formula: sfox ~ roefoxratio + (1 | fam)
> Data: manu2
> Family: binomial(logit link)
> AIC BIC logLik deviance
> 1.797693e+308 1.797693e+308 -8.988466e+307 1.797693e+308
> Random effects:
> Groups Name Variance Std.Dev.
> fam (Intercept) 17.149 4.1412
> # of obs: 128, groups: fam, 58
>
> Estimated scale (compare to 1) 0.5940245
>
> Fixed effects:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -2.60841 1.06110 -2.45820 0.01396 *
> roefoxratio 0.51677 0.63866 0.80915 0.41843
>
> I suspect this may be due to a local maximum in the ML-fitting, since:
>
>> mod at logLik
> 'log Lik.' -8.988466e+307 (df=4)
>
> However,
>
>> mod at deviance
> ML REML
> 295.4233 295.4562
>
> So, my first question is what this second deviance value represent. I
> have tried to figure out from the lmer-syntax
> (https://svn.r-project.org/R-packages/trunk/Matrix/R/lmer.R)
> but I must admit I have problems with this.
>
> Second, if the very high deviance is due to local maximum, is there a
> general procedure to overcome this problem? I have tried to alter the
> tolerance in the control-parameters. However, I need a very high
> tolerance value in order to get a more reasonable deviance, e.g.
>
>> mod <- lmer(sfox ~ roefoxratio + (1|fam), data=manu2,
> family=binomial,
> control=list(tolerance=sqrt(sqrt(sqrt(sqrt(.Machine$double.eps))))))
>> mod
> Generalized linear mixed model fit using PQL
> Formula: sfox ~ roefoxratio + (1 | fam)
> Data: manu2
> Family: binomial(logit link)
> AIC BIC logLik deviance
> 130.2166 141.6247 -61.10829 122.2166
> Random effects:
> Groups Name Variance Std.Dev.
> fam (Intercept) 15.457 3.9316
> # of obs: 128, groups: fam, 58
>
> Estimated scale (compare to 1) 0.5954664
>
> Fixed effects:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -2.55690 0.98895 -2.58548 0.009724 **
> roefoxratio 0.50968 0.59810 0.85216 0.394127
>
> The tolerance value in this model represent 0.1051 on my machine. Does
> anyone have an advice how to handle such problems? I find the
> tolerance
> needed to achieve reasonable deviances rather high, and makes me
> not too
> confident about the estimates and the model. Using the other methods,
> ("Laplace" or "AGQ") did not help.
>
> My system is windows 2000,
>> version
> _
> platform i386-pc-mingw32
> arch i386
> os mingw32
> system i386, mingw32
> status
> major 2
> minor 2.0
> year 2005
> month 10
> day 06
> svn rev 35749
> language R
>
> Thanks
>
> Ivar Herfindal
>
> By the way, great thanks to all persons contributing to this package
> (and other), it makes my research more easy (and fun).
>
More information about the R-help
mailing list