[R] glm binomial with no successes
Gavin Simpson
gavin.simpson at ucl.ac.uk
Wed Feb 27 18:00:33 CET 2008
[Following up for my own personal education so a bit OT!]
Naively, I would have thought that package multcomp would be of use
here. So I tried, for my own comprehension and education, to answer the
OP's question using multcomp. Here's what I got:
## make this reproducible (I hope)
set.seed(1234)
s <- c(rpois(8, 4), rep(0, 4))
f <- rpois(12, 30)
tr <- gl(3, 4)
sf <- cbind(s,f)
## fit the glm
mod <- glm(sf ~ tr, family=binomial)
summary(mod) ## tr2 and tr3 not different from reference level tr1
anova(mod, test = "Chisq") ## tr is signif
## multiple comparison of levels of tr
require(multcomp)
mod.glht <- glht(mod, linfct = mcp(tr = "Tukey"))
mod.glht
summary(mod.glht)
If I interpret this correctly, both summary(mod) and summary(mod.glht)
suggest that there are no significant differences between the 3 levels
of tr, but that tr, as a whole, is better than the Null model (as shown
by anova(mod) )?
Is my interpretation correct, for this specific example, or am I abusing
multcomp and statistics in this case?
Thanks for your time and indulgence of a more statistically-related than
R-related question
All the best,
G
On Wed, 2008-02-27 at 12:51 +0100, juli pausas wrote:
> Thank you very much for your reply.
> Then I understand that would not be correct to perform the test in
> summary for testing the significance of the different levels of a
> factor in relation to the first level, including when there are more
> than 2 levels, as in my real case; at least for binomial regressions.
> So here a more close-to-real example, with a 3-level factor
>
> s <- c(rpois(8, 4), rep(0, 4))
> f <- rpois(12, 30)
> tr <- gl(3, 4)
> sf <- cbind(s,f)
> drop1(glm(sf ~ tr, family="binomial"), test="Chisq") # significant
> summary(glm(sf ~ tr, family="binomial")) # the 3rd level
> is not significant from the 1st
>
> So I understand that I need to explite the data and perform the two
> tests separately:
>
> drop1(glm(sf ~ tr, family="binomial", subset=(tr %in% c("1", "2"))),
> test="Chisq") # ns as expected
>
> drop1(glm(sf ~ tr, family="binomial", subset=(tr %in% c("1", "3"))),
> test="Chisq") # significant, as expected
>
> Is this the correct approach?
> Many thanks
>
> Juli
>
> On Wed, Feb 27, 2008 at 12:13 PM, Prof Brian Ripley
> <ripley at stats.ox.ac.uk> wrote:
> > On Wed, 27 Feb 2008, juli pausas wrote:
> >
> > > Dear all,
> > > I have a question on glm, family binomial. I do not see significant
> > > differences between the levels of a factor (treatment) if all data for
> > > a level is 0; and replacing a 0 for a 1 (in fact reducing the
> > > difference), then I detect the significant difference that I expected.
> >
> > This is because you are using the wrong test, one with negligible power.
> > See MASS4 pp.197-8 -- you need to use the LRT, as in
> >
> > > drop1(glm(sf ~ tr, family="binomial"), test="Chisq")
> > Single term deletions
> >
> > Model:
> > sf ~ tr
> > Df Deviance AIC LRT Pr(Chi)
> > <none> 1.595 17.730
> > tr 1 24.244 38.379 22.649 1.944e-06
> >
> > (and in your example you can replace 'drop1' by 'anova').
> >
> >
> > > Is there a way to overcome this problem? or this is an expected
> > > behaviour ? Here is an example:
> > >
> > > s <- c(2,4,4,5,0,0,0,0)
> > > f <- c(31,28,28,28,32,37,34,35)
> > > tr <- gl(2, 4)
> > > sf <- cbind(s,f) # numbers of successes and failures
> > > summary(glm(sf ~ tr, family="binomial")) # tr ns
> > >
> > > sf[8,1] <- 1
> > > summary(glm(sf ~ tr, family="binomial")) # tr significative **
> > >
> > > Thanks for any suggestion
> > >
> > > Juli
> > >
> > > --
> > > http://www.ceam.es/pausas
> > >
> > > ______________________________________________
> > > R-help at r-project.org 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.
> > >
> >
> > --
> > Brian D. Ripley, ripley at stats.ox.ac.uk
> > Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
> > University of Oxford, Tel: +44 1865 272861 (self)
> > 1 South Parks Road, +44 1865 272866 (PA)
> > Oxford OX1 3TG, UK Fax: +44 1865 272595
> >
>
>
>
--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Dr. Gavin Simpson [t] +44 (0)20 7679 0522
ECRC, UCL Geography, [f] +44 (0)20 7679 0565
Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
More information about the R-help
mailing list