## Sourcing ci methods
source("http://myweb.uiowa.edu/pbreheny/571/f14/labs/binom.bayes.R")
prop.wald <- function(x, n, level=0.95) {
  pi.hat <- x/n
  SE <- sqrt(pi.hat*(1-pi.hat)/n)
  z <- qnorm(c((1-level)/2, 1-(1-level)/2))
  pi.hat + z*SE
}
prop.test(31,39) ## Just to show how prop.test works

## Coverage simulation loop
source("http://myweb.uiowa.edu/pbreheny/571/f14/labs/displayProgressBar.R")
N <- 10000
n <- 20
pi <- c(0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.5, 0.68, 0.84, 0.92, 0.96, 0.98, 0.99)
coverage <- matrix(NA, length(pi), 4, dimnames=list(pi, c("CP", "Bayes", "Score", "Wald")))
for (j in 1:length(pi)) {
  x <- rbinom(N, prob=pi[j], size=n)
  covered <- matrix(NA, N, 4)
  for (i in 1:N) {
    ci <- binom.test(x[i], n)$conf
    covered[i,1] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
    ci <- binom.bayes(x[i], n)$ci.hpd
    covered[i,2] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
    ci <- prop.test(x[i], n, correct=FALSE)$conf
    covered[i,3] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
    ci <- prop.wald(x[i], n)
    covered[i,4] <- (ci[1] < pi[j]) & (pi[j] < ci[2])
  }
  coverage[j,] <- apply(covered, 2, mean)
  displayProgressBar(j, length(pi))
}
col <- hcl(seq(15, 375, len = 5), 150, 60)[1:4] ## Setting up some colors
matplot(pi, coverage, type="l", lwd=3, lty=1, col=col)
legend("bottom", legend=colnames(coverage), col=col, lwd=3)
matplot(pi, coverage, type="l", lwd=3, lty=1, col=col, ylim=c(0.8,1))
legend("bottom", legend=colnames(coverage), col=col, lwd=3)
abline(h=0.95, col="gray", lwd=3)

## Cystic fibrosis study
cf <- read.delim("http://myweb.uiowa.edu/pbreheny/data/cystic-fibrosis.txt")
Diff <- cf$Placebo - cf$Drug
hist(Diff, col="gray", border="white")
t.test(Diff)

## NHANES lipid data
lipids <- read.delim("http://myweb.uiowa.edu/pbreheny/data/lipids.txt")
pop <- lipids$TRG

## NHANES sim mean
N <- 10000          ## Number of simulations to run
n <- 25             ## Sample size
covered <- numeric(N)
for (i in 1:N) {
  sam <- sample(pop, n)
  covered[i] <- t.test(sam, mu=mean(pop))$p.value > 0.05
}
mean(covered)

## NHANES sim variance
N <- 10000          ## Number of simulations to run
n <- 25             ## Sample size
covered <- numeric(N)
for (i in 1:N) {
  sam <- sample(pop, n)
  ci <- var(sam)*(n-1)/qchisq(c(0.975,0.025), n-1)
  covered[i] <- (ci[1] < var(pop)) & (var(pop) < ci[2])
}
mean(covered)