[R] small bug in binom.test?

[Jerome Goudet]

At 22.01.2003  08:45 +0000, you wrote:
>Why do you think that?
>
>The problems binom.test(49,50,0.5) and binom.test(51,100,0.5) are
>symmetrical, so one would expect the same results for a two-sided test.
>
>The problem I guess is how a two-sided test is defined for a discrete
>distribution.  For one-sided tests one would use the probability of X >=11
>or X <= 9, and those are not equal.  For a two-sided test the code
>attempts to find a point in the opposite tail with at least as large a
>tail probability, and adds on that tail probability.  Thus for
>binom.test(11,100,p=0.1) it used P(X < 9 || X >= 11), and for
>binom.test(9,100,p=0.1) it used P(X <= 9 || X > 10), if I followed the
>code right.

I would have defined the two sided test as P(X<=9 || X>=11) (checking of 
course that if the two values are equal, the probability is not counted 
twice). Is this wrong? I've got the following code carrying out such a test:


binomial.test<-function(obs,size,p=0.5,alpha=0.05,alternative="two.sided",plotit=T){
#calcule la probabilité d'obtenir obs succès parmi size tirage sous 
l'hypothèse que la probabilité de succès est p (H0 prob=p)
#On peut spécifier si on souhaite un test bilatéral 
(alternative="two.sided")(H1 prob<>p)
#ou unilatéral (alternative="greater" pour H1 prob>p ou alternative="less" 
pour H1 prob<p)
#l'intervalle de confiance à (1-alpha)% n'est donné que si 
alternative="two.sided"
#exemple: binomial.test(23,257,prob=0.1) #voir poly p 11.2;
         p.hat=obs/size
         cent=floor(p*size) #l'espérance de la distribution
         ties=dbinom(cent,size,p)==dbinom(cent+1,size,p) #si 2 valeurs 
partagent la + haute probabilité
         xl=ifelse(obs<cent,obs,ifelse(ties,cent-(obs-(cent+1)),cent-(obs-cent)))
         xh=ifelse(obs<cent,ifelse(ties,(cent+1)+(cent-obs),cent+(cent-obs)),obs)
         if (plotit){
                 plot(0:size,dbinom(0:size,size,p),type="h")#type="h" 
signifie des barres verticales depuis l'axe des x
                      if (alternative=="greater")
                         points(obs:size,dbinom(obs:size,size,p),type="h",col="red",lwd=3)#col 
précise la couleur et lwd la largeur des traits
                      if (alternative=="less")
                         points(0:obs,dbinom(0:obs,size,p),type="h",col="red",lwd=3)
                      if (alternative=="two.sided")
                         points(c(0:xl,xh:size),dbinom(c(0:xl,xh:size),size,p),type="h",col="red",lwd=3)
         }
         if (alternative=="greater")
            return(list(p.hat=p.hat,pval=pbinom(obs-1,size,p,lower.tail=F)))
         if (alternative=="less")
            return(list(p.hat=p.hat,pval=pbinom(obs,size,p)))
         if (alternative=="two.sided"){
                 pval=ifelse(xl!=xh,
                             pbinom(xl,size,p)+pbinom(xh-1,size,p,lower.tail=F),
                             pbinom(xl,size,p)+pbinom(xh,size,p,lower.tail=F))
                 ic.li=qbinom(alpha/2,size,p.hat)/size
                 ic.ls=qbinom(1-alpha/2,size,p.hat)/size
                 return(list(p.hat=p.hat,pval=pval,ic.li=ic.li,ic.ls=ic.ls))
         }
}

>The great thing about R is that you can do
>
>debug(binom.test)
>
>and follow the calculations.
>
>
>On Wed, 22 Jan 2003, Jerome Goudet wrote:
>
> > Hi all,
> >
> > I am wondering whether there is a small bug in the binom.test function of
> > the ctest library (I'm using R 1.6.0 on windows 2000, but Splus 2000 seems
> > to have the same behaviour). Or perhaps I've misunderstood something.
> >
> > the command binom.test(11,100,p=0.1) and binom.test(9,100,p=0.1) give
> > different p-values (see below).  As 9 and 11 are equidistant from 10, the
> > mean of the distribution, I would have thought that the probabilities
> > should be the same. For instance, binom.test(49,50,0.5) and
> > binom.test(51,100,0.5) do give the same results. Any help wouldbe
> > appreciated. Jerome
> >
> >
> >
> >  > binom.test(11,100,0.1)
> >
> >          Exact binomial test
> >
> > data:  11 and 100
> > number of successes = 11, number of trials = 100, p-value = 0.7377
> > alternative hypothesis: true probability of success is not equal to 0.1
> > 95 percent confidence interval:
> >   0.05620702 0.18830113
> > sample estimates:
> > probability of success
> >                    0.11
> >
> >  > binom.test(9,100,0.1)
> >
> >          Exact binomial test
> >
> > data:  9 and 100
> > number of successes = 9, number of trials = 100, p-value = 0.8681
> > alternative hypothesis: true probability of success is not equal to 0.1
> > 95 percent confidence interval:
> >   0.04198360 0.16398226
> > sample estimates:
> > probability of success
> >                    0.09
> >
> >
> > Jerome GOUDET
> > Institut d'Ecologie, Bat. Biologie
> > Uni. Lausanne , CH-1015 Lausanne
> > Switzerland
> > Tel: +41 21 692 42 42    Fax: +41 21 692 42 65
> > Secr:+41 21 692 42 60
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > http://www.stat.math.ethz.ch/mailman/listinfo/r-help
> >
>
>--
>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