![[banner]](images/banner.png "Batsto Lake, Wharton State Forest, New Jersey")
Advertisement
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.
Examples in Summary and Analysis of Extension Program Evaluation
SAEEPER: Goodness-of-Fit Tests for Nominal Variables
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")}
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
# # #
The sign test is described in the Wilcoxon Signed-rank Test chapter.
See example below in the “Examples” section.
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
")
D1 = read.table(textConnection(Input),header=TRUE)
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
# # #
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
")
Data = read.table(textConnection(Input),header=TRUE)
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 is shown in the “Chi-square Goodness-of-Fit” section.
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
")
Gus = read.table(textConnection(Input),header=TRUE)
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
# # #