download.file("http://www.openintro.org/stat/data/atheism.RData", destfile = "atheism.RData")
load("atheism.RData")
us12 <- subset(atheism, nationality == "United States" & year == "2012")
us12$nationality <- as.factor(as.character(us12$nationality))
( us12prop <- prop.table(table(us12$nationality, us12$response)))
##
## atheist non-atheist
## United States 0.0499002 0.9500998
inference(us12$response, est = "proportion", type = "ci", method = "theoretical",
success = "atheist")
## Warning: package 'BHH2' was built under R version 4.1.3
## Single proportion -- success: atheist
## Summary statistics:
## p_hat = 0.0499 ; n = 1002
## Check conditions: number of successes = 50 ; number of failures = 952
## Standard error = 0.0069
## 95 % Confidence interval = ( 0.0364 , 0.0634 )
per12 <- subset(atheism, nationality == "Peru" & year == "2012")
per12$nationality <- as.factor(as.character(per12$nationality))
( per12prop <- prop.table(table(per12$nationality, per12$response)) )
##
## atheist non-atheist
## Peru 0.02982601 0.97017399
inference(per12$response, est = "proportion", type = "ci", method = "theoretical", success = "atheist")
## Single proportion -- success: atheist
## Summary statistics:
## p_hat = 0.0298 ; n = 1207
## Check conditions: number of successes = 36 ; number of failures = 1171
## Standard error = 0.0049
## 95 % Confidence interval = ( 0.0202 , 0.0394 )
jap12 <- subset(atheism, nationality == "Japan" & year == "2012")
jap12$nationality <- as.factor(as.character(jap12$nationality))
( jap12prop <- prop.table(table(jap12$nationality, jap12$response)) )
##
## atheist non-atheist
## Japan 0.3069307 0.6930693
inference(jap12$response, est = "proportion", type = "ci", method = "theoretical", success = "atheist")
## Single proportion -- success: atheist
## Summary statistics:
## p_hat = 0.3069 ; n = 1212
## Check conditions: number of successes = 372 ; number of failures = 840
## Standard error = 0.0132
## 95 % Confidence interval = ( 0.281 , 0.3329 )
n <- 1000
p <- seq(0, 1, 0.01)
me <- 2 * sqrt(p * (1 - p)/n)
plot(me ~ p, ylab = "Margin of Error", xlab = "Population Proportion")
p <- 0.1
n <- 1040
p_hats <- rep(0, 5000)
for(i in 1:5000){
samp <- sample(c("atheist", "non_atheist"), n, replace = TRUE, prob = c(p, 1-p))
p_hats[i] <- sum(samp == "atheist")/n
}
hist(p_hats, main = "p = 0.1, n = 1040", xlim = c(0, 0.18))
p_too <- c(0.1, 0.1, 0.02, 0.02)
n_too <- c(1040, 400, 1040, 400)
p_hats_too <- data.frame(c(rep(0, 5000)), c(rep(0, 5000)), c(rep(0, 5000)), c(rep(0, 5000)))
for(i in 1:5000){
samp <- sample(c("atheist", "non_atheist"), n_too[1], replace = TRUE, prob = c(p_too[1], 1-p_too[1]))
p_hats_too[i, 1] <- sum(samp == "atheist")/n_too[1]
samp <- sample(c("atheist", "non_atheist"), n_too[2], replace = TRUE, prob = c(p_too[2], 1-p_too[2]))
p_hats_too[i, 2] <- sum(samp == "atheist")/n_too[2]
samp <- sample(c("atheist", "non_atheist"), n_too[3], replace = TRUE, prob = c(p_too[3], 1-p_too[3]))
p_hats_too[i, 3] <- sum(samp == "atheist")/n_too[3]
samp <- sample(c("atheist", "non_atheist"), n_too[4], replace = TRUE, prob = c(p_too[4], 1-p_too[4]))
p_hats_too[i, 4] <- sum(samp == "atheist")/n_too[4]
}
par(mfrow = c(2, 2))
hist(p_hats_too[, 1], xlab = "", main = "p = 0.1, n = 1040", xlim = c(0, 0.18))
hist(p_hats_too[, 2], xlab = "", main = "p = 0.1, n = 400", xlim = c(0, 0.18))
hist(p_hats_too[, 3], xlab = "", main = "p = 0.02, n = 1040", xlim = c(0, 0.18))
hist(p_hats_too[, 4], xlab = "", main = "p = 0.02, n = 400", xlim = c(0, 0.18))