Re: [R] unexpected result in glm (family=poisson) for data with an only zero response in one factor

From: vito muggeo <vmuggeo_at_dssm.unipa.it>
Date: Wed 13 Sep 2006 - 10:49:59 GMT

Dear Antonin,
It is a statistical problem: the well-known monotone likelihood.

In this case ML estimate does not exist (or equals infinity) and Wald approximations (ob which SE are based) does not hold.

However LRT is valid and provides reliable results.

As far as I know, the only software dealing with monotone likelihood problems in loglinear models is LogXact by cytel corporation.

best,
vito

Antonin Ferry wrote:
> Dear members,
> here is my trouble: My data consists of counts of trapped insects in different attractive traps. I usually use GLMs with a poisson error distribution to find out the differences between my traitments (and to look at other factor effects). But for some dataset where one traitment contains only zeros, GLM with poisson family fail to find any difference between this particular traitment and anyother one (even with traitment that have trapped a lot of insects). GLMs with gaussian family does not seem to have the same problem but GLMs with binomial family does.
> I'm not sure if it is a statistical problem or if it comes from R... in the latter case I think some solution exists (perhaps in the options of the glm() function ?).
> Thank you for your help.
>
>
> Here I figure out an exemple to past in the console:
>
> ## START ##############################################################################
> # Take a data set of counts for two traitments, one containing only zeros
> A=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
> B=c(1,0,0,0,2,1,0,0,1,2,0,0,0,1,2,2,0,1,1,0,1,0,2,1,1,0,1,2,0,1,0,1,1,1,0,1,1,1,0,1)
> traitment=c(rep("A",40),rep("B",40))
> response=c(A,B)
> mydata=data.frame(traitment ,response)
>
>
> # Make a GLM on this dataset , with "family=poisson"
>
> g=glm(response~traitment, data=mydata, family=poisson)
> anova.glm(g,test="Chisq")
> # There is an effect of the traitment ...
>
> summary(g)
> # But traitment A does not differ from traitment B ! ! ! (the pvalue is always close from 1 in such cases)
>
> # Now if you replace only one zero of the A reponse to 1, the GLM works properly:
> mydata[1,2]=1
> g=glm(response~traitment, data=mydata, family=poisson)
> anova.glm(g,test="Chisq")
> summary(g)
> ##################################################################################### END ##
>
>
>
> Antonin Ferry (PhD)
>
> "Laboratoire d'Ecobiologie des Insectes Parasitoides"
> http://www.parasitoides.univ-rennes1.fr
> Université de Renes1, FRANCE
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
====================================
Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 6626240
fax: 091 485726/485612

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Wed Sep 13 21:54:05 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Wed 13 Sep 2006 - 12:30:05 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.