Created
November 18, 2013 03:28
-
-
Save fditraglia/7521997 to your computer and use it in GitHub Desktop.
R Code to Accompany Solutions to Homework #10
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
group <- c('treatment', 'control', 'refused') | |
n.children <- 1000 * c(200, 200, 340) | |
n.polio <- c(57, 142, 157) | |
rate <- n.polio/n.children | |
polio.data <- data.frame(group, n.children, n.polio, rate) | |
polio.data | |
treatment <- subset(polio.data, group == 'treatment') | |
control <- subset(polio.data, group == 'control') | |
estimate <- control$rate - treatment$rate | |
SE <- sqrt( | |
control$rate * (1 - control$rate)/control$n.children | |
+ treatment$rate * (1 - treatment$rate)/control$n.children | |
) | |
ME <- qnorm(1 - 0.05/2) * SE | |
CI <- c(estimate - ME, estimate + ME) | |
estimate * 10^5 | |
ME * 10^5 | |
CI * 10^5 | |
x <- seq(0, 1, 0.01) | |
y <- dunif(x) | |
plot(x, y, main = "Uniform(0,1) Density", type = 'l', ylab = 'f(x)') | |
uniform.means <- replicate(10000, mean(runif(20))) | |
hist(uniform.means, main = "Uniform(0,1), n = 20") | |
x <- seq(0.01, 15, 0.01) | |
y <- dchisq(x, df = 5) | |
plot(x,y, type = 'l', main = "Chi-squared Density, df = 5", ylab = 'f(x)') | |
chisq.means <- replicate(10000, mean(rchisq(20, df = 5))) | |
hist(chisq.means, main = "Chi-squared(5), n = 20") | |
x <- c(0,1) | |
y <- dbinom(x, size = 1, p = 0.3) | |
plot(x, y, type = 'h', main = "Bernoulli(0.3) pmf", ylim = c(0,1), xlim = c(-1, 2), ylab = 'p(x)') | |
bernoulli.means <- replicate(10000, mean(rbinom(20, 1, 0.3))) | |
hist(bernoulli.means, main = "Bernoulli(0.3), n = 20") | |
p.hat <- 0.21 | |
n <- 1247 | |
SE <- sqrt(p.hat * (1 - p.hat)/n) | |
ME <- qnorm(1 - 0.01/2) * SE | |
LCL <- p.hat - ME | |
UCL <- p.hat + ME | |
c(LCL, UCL) | |
n <- 623 | |
p.M <- 0.24 | |
p.W <- 0.19 | |
SE.M <- sqrt(p.M * (1 - p.M)/n) | |
SE.W <- sqrt(p.W * (1 - p.W)/n) | |
SE <- sqrt(SE.M^2 + SE.W^2) | |
ME <- qnorm(1 - 0.01/2) * SE | |
diff <- p.M - p.W | |
LCL <- diff - ME | |
UCL <- diff + ME | |
c(LCL, UCL) | |
n.R <- 547 | |
p.R <- 0.27 | |
SE.R <- sqrt(p.R * (1 - p.R)/n.R) | |
n.O <- 623 | |
p.O <- 0.16 | |
SE.O <- sqrt(p.O * (1 - p.O)/n.O) | |
SE <- sqrt(SE.R^2 + SE.O^2) | |
ME <- qnorm(1 - 0.01/2) * SE | |
diff <- p.R - p.O | |
UCL <- diff + ME | |
LCL <- diff - ME | |
c(LCL, UCL) | |
data.url <- "http://www.ditraglia.com/econ103/survey_clean.csv" | |
survey <- read.csv(data.url) | |
anchoring <- survey[,c("rand.num", "africa.percent")] | |
lo <- subset(anchoring, rand.num == "10")$africa.percent | |
hi <- subset(anchoring, rand.num == "65")$africa.percent | |
lo <- lo[!is.na(lo)] | |
hi <- hi[!is.na(hi)] | |
SE.lo <- sd(lo)/sqrt(length(lo)) | |
SE.hi <- sd(hi)/sqrt(length(hi)) | |
SE <- sqrt(SE.hi^2 + SE.lo^2) | |
ME <- qnorm(0.975) * SE | |
diff <- mean(hi) - mean(lo) | |
LCL <- diff - ME | |
UCL <- diff + ME | |
c(LCL, UCL) | |
data.url <- "http://www.ditraglia.com/econ103/ex0222.csv" | |
test.scores <- read.csv(data.url) | |
head(test.scores) | |
test.men <- subset(test.scores, Gender == 'male')[,-1] | |
means.men <- apply(test.men, 2, mean) | |
var.men <- apply(test.men, 2, var) | |
test.men <- subset(test.scores, Gender == 'male')[,-1] | |
test.women <-subset(test.scores, Gender == 'female')[,-1] | |
means.men <- apply(test.men, 2, mean) | |
var.men <- apply(test.men, 2, var) | |
n.men <- nrow(test.men) | |
means.women <- apply(test.women, 2, mean) | |
var.women <- apply(test.women, 2, var) | |
n.women <- nrow(test.women) | |
diff.means <- means.men - means.women | |
SE <- sqrt(var.women/n.women + var.men/n.men) | |
ME <- qnorm(1 - 0.05/2) * SE | |
LCL <- diff.means - ME | |
UCL <- diff.means + ME | |
CI <- rbind(LCL, UCL) | |
round(diff.means, 2) | |
round(CI, 2) | |
data.url <- "http://www.ditraglia.com/econ103/case0202.csv" | |
twins <- read.csv(data.url) | |
head(twins) | |
mean.affected <- mean(twins$Affected) | |
var.affected <- var(twins$Affected) | |
n.affected <- length(twins$Affected) | |
mean.unaffected <- mean(twins$Unaffected) | |
var.unaffected <- var(twins$Unaffected) | |
n.unaffected <- length(twins$Unaffected) | |
diff.means <- mean.unaffected - mean.affected | |
SE.indep <- sqrt( | |
var.affected/n.affected | |
+ var.unaffected/n.unaffected) | |
ME.indep <- qnorm(1 - 0.05/2) * SE.indep | |
CI.indep <- c(diff.means - ME.indep, diff.means + ME.indep) | |
round(CI.indep, 3) | |
twin.diff <- twins$Unaffected - twins$Affected | |
n.twins <- length(twin.diff) | |
SE.paired <- sqrt(var(twin.diff)/n.twins) | |
ME.paired <- qnorm(1 - 0.05/2) * SE.paired | |
CI.paired <- c(diff.means - ME.paired, diff.means + ME.paired) | |
round(CI.paired, 3) | |
ME.t <- qt(1 - 0.05/2, df = n.twins - 1) * SE.paired | |
CI.paired.t <- c(diff.means - ME.t, diff.means + ME.t) | |
round(CI.paired.t, 3) | |
rbinom(1, size = 50, prob = 0.5)/50 | |
n <- 50 | |
p <- 0.5 | |
N.reps <- 100 | |
p.hat <- rbinom(N.reps, size = n, prob = p)/n | |
ME.hat <- qnorm(0.975) * sqrt(p.hat * (1 - p.hat) / n) | |
LCL.hat <- p.hat - ME.hat | |
UCL.hat <- p.hat + ME.hat | |
CI.hat <- cbind(LCL.hat, UCL.hat) | |
Coverage <- (LCL.hat <= p) & (p <= UCL.hat) | |
Coverage <- sum(Coverage)/N.reps | |
Coverage | |
n <- 10 | |
p <- 0.1 | |
N.reps <- 10000 | |
p.hat <- rbinom(N.reps, size = n, prob = p)/n | |
ME.hat <- qnorm(0.975) * sqrt(p.hat * (1 - p.hat) / n) | |
LCL.hat <- p.hat - ME.hat | |
UCL.hat <- p.hat + ME.hat | |
CI.hat <- cbind(LCL.hat, UCL.hat) | |
Coverage <- (LCL.hat <= p) & (p <= UCL.hat) | |
Coverage <- sum(Coverage)/N.reps | |
Coverage | |
p.tilde <- (n * p.hat + 2) / (n + 4) | |
ME.tilde <- qnorm(0.975) * sqrt(p.tilde * (1 - p.tilde) / (n + 4)) | |
LCL.tilde <- p.tilde - ME.tilde | |
UCL.tilde <- p.tilde + ME.tilde | |
CI.tilde <- cbind(LCL.tilde, UCL.tilde) | |
Cover.tilde <- (LCL.tilde <= p) & (p <= UCL.tilde) | |
Cover.tilde <- sum(Cover.tilde)/N.reps | |
Cover.tilde | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment