# ~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~ CRP 241 Module 5 Day 1 ~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~
# Example 1
# Viability Study Data Example:
# Suppose researchers were interested in determining
# whether or not infants born at 24 weeks have a 50%
# chance of survival at 6 months.
# Single proportion test
prop.test(x=19,n=34,p=0.50,correct=FALSE)
##
## 1-sample proportions test without continuity correction
##
## data: 19 out of 34, null probability 0.5
## X-squared = 0.47059, df = 1, p-value = 0.4927
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.3945390 0.7111652
## sample estimates:
## p
## 0.5588235
# Chi-square goodness of fit test
chisq.test(x=c(19,15),p=c(0.50,0.50),correct=FALSE)
##
## Chi-squared test for given probabilities
##
## data: c(19, 15)
## X-squared = 0.47059, df = 1, p-value = 0.4927
# Example 2
# Firefighter Data Example:
# Suppose researchers were interested in determining
# whether or not death from CHD was equally likely
# across all duties.
# Data Dictionary:
# 1. FFID national firefighter ID
# 2. DUTY activity firefighter was performing when they
# experienced an on-duty death from CHD
# (1 = Fire suppression; 2 = Alarm response & return;
# 3 = Physical Training; 99 = Other)
# Download and load the firefighter dataset used in lecture:
download.file("http://www.duke.edu/~sgrambow/crp241data/firefighter.RData",
destfile = "firefighter.RData",quiet=TRUE,mode="wb",cacheOK=FALSE)
load("firefighter.RData")
# Compute counts for frequency table
freqs <- table(ff_data$DUTY)
freqs
##
## 1 2 3 99
## 144 138 56 111
# Compute proportions for frequency table
props <- prop.table(freqs)
props
##
## 1 2 3 99
## 0.3207127 0.3073497 0.1247216 0.2472160
round(props,2)
##
## 1 2 3 99
## 0.32 0.31 0.12 0.25
# Create bar plot for a descriptive figure
barplot(table(ff_data$DUTY),
ylim=c(0,150),
main='On-Duty Deaths from Heart Disease in Firefighters',
ylab='Frequency',
xlab='On-Duty Activity',
names=c('Fire Supression','Alarm Response','Physical Training','Other'))
box()
# Chi-square goodness of fit test
# - Option 1: Manually type in level counts
chisq.test(x=c(144,138,56,111),
p=c(0.25,0.25,0.25,0.25),
correct=FALSE)
##
## Chi-squared test for given probabilities
##
## data: c(144, 138, 56, 111)
## X-squared = 43.089, df = 3, p-value = 2.356e-09
# - Option 2: Use the ouptut of the table() function
chisq.test(table(ff_data$DUTY),
p=c(0.25,0.25,0.25,0.25),
correct=FALSE)
##
## Chi-squared test for given probabilities
##
## data: table(ff_data$DUTY)
## X-squared = 43.089, df = 3, p-value = 2.356e-09
# Understand the global Chi-square p-value
prop.test(x=144,n=449,p=0.25) # For 1 = Fire suppression
##
## 1-sample proportions test with continuity correction
##
## data: 144 out of 449, null probability 0.25
## X-squared = 11.6, df = 1, p-value = 0.0006596
## alternative hypothesis: true p is not equal to 0.25
## 95 percent confidence interval:
## 0.2781525 0.3663928
## sample estimates:
## p
## 0.3207127
prop.test(x=138,n=449,p=0.25) # For 2 = Alarm response & return
##
## 1-sample proportions test with continuity correction
##
## data: 138 out of 449, null probability 0.25
## X-squared = 7.5731, df = 1, p-value = 0.005924
## alternative hypothesis: true p is not equal to 0.25
## 95 percent confidence interval:
## 0.2653949 0.3526578
## sample estimates:
## p
## 0.3073497
prop.test(x= 56,n=449,p=0.25) # For 3 = Physical Training
##
## 1-sample proportions test with continuity correction
##
## data: 56 out of 449, null probability 0.25
## X-squared = 36.918, df = 1, p-value = 1.232e-09
## alternative hypothesis: true p is not equal to 0.25
## 95 percent confidence interval:
## 0.09631897 0.15972098
## sample estimates:
## p
## 0.1247216
prop.test(x=111,n=449,p=0.25) # For 99 = Other
##
## 1-sample proportions test with continuity correction
##
## data: 111 out of 449, null probability 0.25
## X-squared = 0.0066815, df = 1, p-value = 0.9349
## alternative hypothesis: true p is not equal to 0.25
## 95 percent confidence interval:
## 0.2085254 0.2903144
## sample estimates:
## p
## 0.247216
# End of Program