[R-sig-ME] quasi-binomial family in lme4
Jarrod Hadfield
j.hadfield at ed.ac.uk
Wed Nov 10 10:53:32 CET 2010
Hi,
It seems worrying to me that the Herd:Id estimate and the Id estimate
are exactly the same. I get a (slightly) different estimate with
version 0.999375-35 on Linux, but again the estimates are identical.
MCMCglmm fitting the same model gives a very different answer.
Cheers,
Jarrod
> summary(gm2)
Generalized linear mixed model fit by the Laplace approximation
Formula: incidence ~ period + (1 | Id/herd)
Data: cbpp
AIC BIC logLik deviance
102.2 114.4 -45.11 90.21
Random effects:
Groups Name Variance Std.Dev.
herd:Id (Intercept) 0.29654 0.54456
Id (Intercept) 0.29654 0.54456
Number of obs: 56, groups: herd:Id, 56; Id, 56
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1129 0.2477 4.492 7.04e-06 ***
period2 -1.1969 0.4185 -2.860 0.004233 **
period3 -1.4223 0.4381 -3.246 0.001169 **
period4 -2.0049 0.5292 -3.788 0.000152 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Correlation of Fixed Effects:
(Intr) perid2 perid3
period2 -0.592
period3 -0.565 0.335
period4 -0.468 0.277 0.265
Using a parameter expanded prior in order to improve mixing as the
herd variance approaches zero:
prior<-list(R=list(V=1, nu=0), G=list(G1=list(V=1, nu=1, alpha.mu=0,
alpha.V=1000)))
mcmc2<-MCMCglmm(incidence~period, random=~herd, family="poisson", data=cbpp)
> summary(mcmc2)
Iterations = 12991
Thinning interval = 3001
Sample size = 1000
DIC: 175.9590
G-structure: ~herd
post.mean l-95% CI u-95% CI eff.samp
herd 0.01433 1.067e-16 0.0713 61.58
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.7347 0.1564 1.421 222.9
Location effects: incidence ~ period
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 1.0749 0.4584 1.5987 844.5 0.002 **
period2 -1.2054 -2.0862 -0.3556 771.8 0.004 **
period3 -1.4568 -2.4362 -0.4839 674.8 0.004 **
period4 -2.0584 -3.1851 -0.9430 446.5 <0.001 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Quoting Gosselin Frederic <frederic.gosselin at cemagref.fr>:
> Dear Colleague,
>
> good question:
>
> the commands are under R2.5.1:
> ************************************************************************
> library(lme4)
> cbppbis<-cbind.data.frame(cbpp,Id=as.factor(1:dim(cbpp)[1]))
> gm1 <- lmer(incidence ~ period + (1 | herd), family = quasipoisson,
> data = cbpp)
> summary(gm1)
>
> gm2 <- lmer(incidence ~ period + (1 | Id/herd), family = poisson,
> data = cbppbis)
> summary(gm2)
> **************************************************************************
> (fope it is declared in the good order for the random effects in gm2)
>
> The results do show a slight discrepancy between both methods:
>
> *************************************************************************
>> summary(gm1)
> Generalized linear mixed model fit using Laplace
> Formula: incidence ~ period + (1 | herd)
> Data: cbpp
> Family: quasipoisson(log link)
> AIC BIC logLik deviance
> 112.2 122.3 -51.11 102.2
> Random effects:
> Groups Name Variance Std.Dev.
> herd (Intercept) 0.35085 0.59233
> Residual 1.40470 1.18520
> number of obs: 56, groups: herd, 15
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 1.2812 0.2200 5.824
> period2 -1.1240 0.3315 -3.391
> period3 -1.3203 0.3579 -3.689
> period4 -1.9477 0.4808 -4.051
>
> Correlation of Fixed Effects:
> (Intr) perid2 perid3
> period2 -0.339
> period3 -0.314 0.219
> period4 -0.233 0.163 0.151
>
>
>
>
>> summary(gm2)
> Generalized linear mixed model fit using Laplace
> Formula: incidence ~ period + (1 | Id/herd)
> Data: cbppbis
> Family: poisson(log link)
> AIC BIC logLik deviance
> 102.2 114.4 -45.11 90.21
> Random effects:
> Groups Name Variance Std.Dev.
> herd:Id (Intercept) 0.29608 0.54413
> Id (Intercept) 0.29608 0.54413
> number of obs: 56, groups: herd:Id, 56; Id, 56
>
> Estimated scale (compare to 1 ) 0.9249959
>
> Fixed effects:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) 1.1149 0.2476 4.503 6.69e-06 ***
> period2 -1.2013 0.4184 -2.871 0.004089 **
> period3 -1.4224 0.4378 -3.249 0.001159 **
> period4 -2.0089 0.5294 -3.795 0.000148 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Correlation of Fixed Effects:
> (Intr) perid2 perid3
> period2 -0.592
> period3 -0.565 0.335
> period4 -0.468 0.277 0.264
> ************************************************************************
>
> The estimates and the standard errors are not exactly the same -
> which might not be unlogical given that the relationship between
> variance and mean is not the same in both models. They however are
> not very far one from the other.
>
> Of course, this needs further investigation.
>
> I remember of a paper that motivated the use of the quasi-poisson
> method based on empirical relationships between the residual
> variance and the mean:
>
> Ver Hoef J.M. et Boveng P.L., 2007, Quasi-poisson vs. negative
> binomial regression: How should we model overdispersed count data?,
> Ecology, 88, 11, p. 2766-2772.
>
> Sincerely,
>
> Frédéric
>
> -----Message d'origine-----
> De : John Maindonald [mailto:john.maindonald at anu.edu.au]
> Envoyé : mercredi 10 novembre 2010 09:51
> À : Gosselin Frederic
> Cc : r-sig-mixed-models at r-project.org; tiflo at csli.stanford.edu
> Objet : Re: [R-sig-ME] quasi-binomial family in lme4
>
> I wonder if you have compared the results that you quote with the
> result you get with observation level random effects in a poisson
> model.
>
> As I see it, use of observation level random effects should, unless
> there is evidence that a multiplicative effect on the scale of the
> response is a better fit, replace use of the quasi- models in glm()
> as well as in generalised linear mixed models.
>
> John Maindonald email: john.maindonald at anu.edu.au
> phone : +61 2 (6125)3473 fax : +61 2(6125)5549
> Centre for Mathematics & Its Applications, Room 1194, John Dedman
> Mathematical Sciences Building (Building 27) Australian National
> University, Canberra ACT 0200.
> http://www.maths.anu.edu.au/~johnm
>
> On 10/11/2010, at 6:39 PM, Gosselin Frederic wrote:
>
>> Hi Florian,
>>
>> a different perpsective on the quasi-likelihood debate - that comes
>> out sporadically on this list:
>>
>> (i) I globally agree with the previous repliers that a fully
>> probabilistic solution looks better - at least aesthetically - than a
>> quasi-likelihood;
>>
>> (ii) however, as I have already mentioned on the list (cf. below),
>> earlier versions of lme4 give much more sensible results than the
>> latest versions:
>> http://markmail.org/message/s4abxhhdacqjkunm
>>
>> This is why in the following papers:
>> Elek Z., Dauffy-Richard & Gosselin F., 2010, Carabid species
>> responses to hybrid poplar plantation in floodplains in France,
>> Forest Ecology and Management, 260, 9, p. 1446-1455.
>>
>> and
>>
>> Vuidot A., Paillet Y., Archaux F. & Gosselin F. (In Press)
>> Influence of tree characteristics and forest management on tree
>> microhabitats in France, Biological Conservation.
>>
>> we used version the R version 2.5.1 and the associated lme4 version
>> (here with quasi-poisson, not quasi-bionomial).
>>
>> Hope this helps.
>>
>> Sincerely,
>>
>> Frédéric Gosselin
>> Engineer & Researcher (PhD) in Forest Ecology Cemagref Domaine des
>> Barres F-45290 Nogent sur Vernisson France
>>
>> http://www.cemagref.fr/les-contacts/les-pages-personnelles-professionn
>> elles/gosselin-frederic/english-short-scientific-cv
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list