An R Companion for the Handbook of Biological Statistics

Salvatore S. Mangiafico

Exact Test of Goodness-of-Fit

The exact test goodness-of-fit can be performed with the binom.test function in the native stats package.  The arguments passed to the function are: the number of successes, the number of trials, and the hypothesized probability of success.  The probability can be entered as a decimal or a fraction.  Other options include the confidence level for the confidence interval about the proportion, and whether the function performs a one-sided or two-sided (two-tailed) test.  In most circumstances, the two-sided test is used.

Packages used in this chapter

The following commands will install these packages if they are not already installed:

if(!require(XNomial)){install.packages("XNomial")}
if(!require(pwr)){install.packages("pwr")}
if(!require(BSDA)){install.packages("BSDA")}

Introduction

When to use it

Null hypothesis

See the Handbook for information on these topics.

How the test works

Binomial test examples

### --------------------------------------------------------------
### Cat paw example, exact binomial test, pp. 30–31
### --------------------------------------------------------------
### In this example:
###   2 is the number of successes
###   10 is the number of trials
###   0.5 is the hypothesized probability of success

dbinom(2, 10, 0.5)            # Probability of single event only!
#   Not binomial test!

[1] 0.04394531

binom.test(2, 10, 0.5,
alternative="less",       # One-sided test
conf.level=0.95)

p-value = 0.05469

binom.test(2, 10, 0.5,
alternative="two.sided",  # Two-sided test
conf.level=0.95)

p-value = 0.1094

#     #     #

Probability density plot

### --------------------------------------------------------------
### Probability density plot, binomial distribution, p. 31
### --------------------------------------------------------------
# In this example:
#   You can change the values for trials and prob
#   You can change the values for xlab and ylab

trials = 10
prob = 0.5

x = seq(0, trials)                   # x is a sequence, 1 to trials
y = dbinom(x, size=trials, p=prob)   # y is the vector of heights

barplot (height=y,
names.arg=x,
xlab="Number of uses of right paw",
ylab="Probability under null hypothesis")

#     #     #

Comparing doubling a one-sided test and using a two-sided test

### --------------------------------------------------------------
### Cat hair example, exact binomial test, p. 31
32
###  Compares performing a one-sided test and doubling the
###    probability, and performing a two-sided test
### --------------------------------------------------------------

binom.test(7, 12, 3/4,
alternative="less",
conf.level=0.95)

p-value = 0.1576

Test = binom.test(7, 12, 3/4,             # Create an object called
alternative="less",     #  Test with the test
conf.level=0.95)        #  results.

2 * Test\$ p.value               # This extracts the p-value from the
#   test result, we called Test
#   and multiplies it by 2

[1] 0.3152874

binom.test(7, 12, 3/4, alternative="two.sided", conf.level=0.95)

p-value = 0.1893      # Equal to the "small p values" method in the Handbook

#     #     #

Sign test

The sign test is described in the Wilcoxon Signed-rank Test chapter.

See example below in the “Examples” section.

Post-hoc test

Post-hoc example with manual pairwise tests

A multinomial test can be conducted with the xmulti function in the package XNomial.  This can be followed with the individual binomial tests for each proportion, as post-hoc tests.

### --------------------------------------------------------------
### Post-hoc example, multinomial and binomial test, p. 33
### --------------------------------------------------------------

observed = c(72, 38, 20, 18)
expected = c(9, 3, 3, 1)

library(XNomial)
xmulti(observed,
expected,
detail = 2)         # 2: Reports three types of p-value

P value  (LLR)  =  0.003404  # log-likelihood ratio

P value (Prob)  =  0.002255  # exact probability

P value (Chisq) =  0.001608  # Chi-square probability

### Note last p-value below agrees with Handbook

successes   = 72
total       = 148
numerator   = 9
denominator = 16

binom.test(successes, total, numerator/denominator,
alternative="two.sided", conf.level=0.95)

p-value = 0.06822

successes   = 38
total       = 148
numerator   = 3
denominator = 16

binom.test(successes, total, numerator/denominator,
alternative="two.sided", conf.level=0.95)

p-value = 0.03504

successes   = 20
total       = 148
numerator   = 3
denominator = 16

binom.test(successes, total, numerator/denominator,
alternative="two.sided", conf.level=0.95)

p-value = 0.1139

successes   = 18
total       = 148
numerator   = 1
denominator = 16

binom.test(successes, total, numerator/denominator,
alternative="two.sided", conf.level=0.95)

p-value = 0.006057

#     #     #

Post-hoc test alternate method with custom function

When you need to do multiple similar tests, however, it is often possible to use the programming capabilities in R to do the tests more efficiently.  The following example may be somewhat difficult to follow for a beginner.  It creates a data frame and then adds a column called p.Value that contains the p-value from the binom.test performed on each row of the data frame.

### --------------------------------------------------------------
### Post-hoc example, multinomial and binomial test, p. 33
###    Alternate method for multiple tests
### --------------------------------------------------------------

Input =("
Successes Total Numerator Denominator
72        148   9         16
38        148   3         16
20        148   3         16
18        148   1         16
")

Fun = function (x){
binom.test(x["Successes"],x["Total"],
x["Numerator"]/x["Denominator"])\$ p.value
}

D1\$ p.Value = apply(D1, 1, Fun)

D1

Successes Total Numerator Denominator     p.Value

1        72   148         9          16 0.068224131

2        38   148         3          16 0.035040215

3        20   148         3          16 0.113911643

4        18   148         1          16 0.006057012

#     #     #

Intrinsic hypothesis

Assumptions

See the Handbook for information on these topics.

Examples

Binomial test examples

### --------------------------------------------------------------
### Parasitoid examples, exact binomial test, p. 34
### --------------------------------------------------------------

binom.test(10, (17+10), 0.5,
alternative="two.sided",
conf.level=0.95)

p-value = 0.2478

binom.test(36, (7+36), 0.5,
alternative="two.sided",
conf.level=0.95)

p-value = 8.963e-06

#     #     #

### --------------------------------------------------------------
###  Drosophila example, exact binomial test, p. 34
### --------------------------------------------------------------

binom.test(140, (106+140), 0.5,
alternative="two.sided",
conf.level=0.95)

p-value = 0.03516

#     #     #

Sign test example

The following is an example of the two-sample dependent-samples sign test.  The data are arranged as a data frame in which each row contains the values for both measurements being compared for each experimental unit.  This is sometimes called “wide format” data.  The SIGN.test function in the BSDA package is used.  The option md=0 indicates that the expected difference in the medians is 0 (null hypothesis).  This function can also perform a one-sample sign test.

### --------------------------------------------------------------
###  Tree beetle example, two-sample sign test, p. 34
35
### --------------------------------------------------------------

Input =("
Row  Angiosperm.feeding  A.count  Gymonsperm.feeding   G.count
1    Corthylina           458     Pityophthorus           200
2    Scolytinae          5200     Hylastini_Tomacini      180
3    Acanthotomicus_P     123     Orhotomicus              11
4    Xyleborini_D        1500     Ipini                   195
5    Apion               1500     Antliarhininae           12
6    Belinae              150     Allocoryninae_Oxycorinae 30
7    H_Curculionidae    44002     Nemonychidae             85
8    H_Cerambycidae     25000     Aseminae_Spondylinae     78
9    Megalopodinae        400     Palophaginae              3
10   H_Chrysomelidae    33400     Aulocoscelinae_Orsod     26
")

library(BSDA)

SIGN.test(x = Data\$ A.count,
y = Data\$ B.count,
md = 0,
alternative = "two.sided",
conf.level = 0.95)

p-value = 0.001953

#     #     #

Binomial test examples

### --------------------------------------------------------------
###  First Mendel example, exact binomial test, p. 35
### --------------------------------------------------------------

binom.test(428, (428+152), 0.75, alternative="two.sided",
conf.level=0.95)

p-value = 0.5022          # Value is different than in the Handbook

#     See next example

#     #     #

### --------------------------------------------------------------
###  First Mendel example, exact binomial test, p. 35
###     Alternate method with XNomial package
### --------------------------------------------------------------

observed = c(428, 152)
expected = c(3, 1)

library(XNomial)

xmulti(observed,
expected,
detail = 2)             # 2: reports three types of p-value

P value  (LLR)  =  0.5331   # log-likelihood ratio

P value (Prob)  =  0.5022   # exact probability

P value (Chisq) =  0.5331   # Chi-square probability

### Note last p-value below agrees with Handbook

#     #     #

Multinomial test example

### --------------------------------------------------------------
###  Second Mendel example, multinomial exact test, p. 35
36
###     and SAS example, p. 38
### --------------------------------------------------------------

observed = c(315, 108, 101, 32)
expected = c(9, 3, 3, 1)

library(XNomial)

xmulti(observed,
expected,
detail = 2)              # reports three types of p-value

P value  (LLR)  =  0.9261    # log-likelihood ratio

P value (Prob)  =  0.9382    # exact probability

P value (Chisq) =  0.9272    # Chi-square probability

### Note last p-value below agrees with Handbook,

###   and agrees with SAS Exact Pr>=ChiSq

#     #     #

Graphing the results

Graphing is shown in the “Chi-square Goodness-of-Fit section.

Similar tests

The G–test goodness-of-fit and chi-square goodness-of-fit are presented elsewhere in this book.

How to do the test

Binomial test example where individual responses are counted

### --------------------------------------------------------------
### Cat paw example from SAS, exact binomial test, pp. 36–37
###      When responses need to be counted
### --------------------------------------------------------------

Input =("
Paw
right
left
right
right
right
right
left
right
right
right
")

Successes = sum(Gus\$ Paw == "left")      # Note the == operator
Failures  = sum(Gus\$ Paw == "right")

Total = Successes + Failures

Expected = 0.5

binom.test(Successes, Total, Expected,
alternative="less",           # One-sided test!
conf.level=0.95)

p-value = 0.05469

binom.test(Successes, Total, Expected,
alternative="two.sided",      # Two-sided test
conf.level=0.95)

p-value = 0.1094

#     #     #

Other SAS examples

R code for the other SAS example is shown in the examples in previous sections.

Power analysis

Power analysis for binomial test

### --------------------------------------------------------------
### Power analysis, binomial test, cat paw, p. 38
### --------------------------------------------------------------

P0 = 0.50
P1 = 0.40
H  = ES.h(P0,P1)               # This calculates effect size

library(pwr)

pwr.p.test(
h=H,
n=NULL,                  # NULL tells the function to
sig.level=0.05,          #     calculate this value
power=0.80,              # 1 minus Type II probability
alternative="two.sided")

n = 193.5839                # Slightly different than in Handbook

#     #     #