## Precept 3
library(dplyr)
library(ggplot2)

##--------------------------------------
## rnorm, pnorm, qnorm
##--------------------------------------
## r: random draws
## p: cumultive density function
## q: inverse cdf (quantiles)
norm <- rnorm(1000, mean = 0, sd = 1) # 1000 draws of standard normal
mean(norm)
sd(norm)

## plot pdf
ggplot() + geom_density(aes(x = norm))

pnorm(0.7, mean = 0, sd = 1) # CDF: probability of X<=0.7 

## create data to plot standard normal cdf
x <- seq(from = 0.01, to = 1, by = 0.01)
y <- rep(NA, length(x))
for (i in 1:length(x)) {
  y[i] <- qnorm(x[i])
}
cdf.data <- data.frame(xval = y, yval = x)

cdf.sn <- ggplot(data = cdf.data, aes(x = xval, y = yval))
cdf.sn + geom_line() +
  scale_x_continuous("x") +
  scale_y_continuous("F(x)") +
  labs(title = "CDF of Standard Normal") +
  theme_bw()

qnorm(0.025, mean = 0, sd = 1, lower.tail = TRUE) # value of z for 95% CI
qnorm(1 - 0.025, mean = 0, sd = 1, lower.tail = TRUE) # value of z for 95% CI
qnorm(.975, mean = 0, sd = 1)

qnorm(0.025, mean = 0, sd = 1, lower.tail = TRUE)
qnorm(1 - 0.05, mean = 0, sd = 1, lower.tail = TRUE)

##--------------------------------------
## Random Sample - Fulton Example
##--------------------------------------
set.seed(12345)

# Load data
load("fulton.RData")

# What does this data look like? 
head(precincts)[1:6]
head(precincts)[7:11]

names(precincts) # variable names
dim(precincts) # number of obs and number of variables

# GOAL: ESTIMATE THE PERCENT OF VOTERS WHO ARE BLACK IN FULTON COUNTY
# First, let's take a look at the density function for 'black'

library(ggplot2)
p.dens.bl <- ggplot()
p.dens.bl + geom_density(data = precincts, aes(x = black)) +
  xlim(-0.5, 1.5) + 
  labs(title = "Percent Black in Fulton County") +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

# Let's calculate the 'true' mean of percent black across precincts as a reference 
true.mean <- mean(precincts$black)
true.mean

# Now let's say we want to take a random sample of size 40 (precincts)
n.precincts <- nrow(precincts) # number of rows
N <- 40 # sample size

my.poll <- sample_n(precincts, size = N, replace = FALSE) # note sample_n is comes with dplyr

mu.b <- mean(my.poll$black) # sample mean
sd.b <- sd(my.poll$black) # sample standard deviation
mu.b
sd.b

# Let's calculate the two-sided 95% confidence interval
# How do we find z_{\alpha/2}?
# \alpha = 1-.95
# \alpha = .05

?qnorm
qnorm(0.05 / 2, mean = 0 , sd = 1, lower.tail = FALSE)

mu.b + 1.96*sd.b/sqrt(N)
mu.b - 1.96*sd.b/sqrt(N)

##--------------------------------------
## Many Samples! 
##--------------------------------------

## We're going to pretend we have 20 independent researchers out in
## the field polling precincts in Fulton. Each of the 20 will
## poll 40 precincts and then calculate the mean and sd and 95% CI.

set.seed(12345)		## again, remember seed setting!
N <- 40			## sample size (40 precincts)
pollsters <- 20		## number of pollsters (taking 20 random samples)

## containers
poll.means <- matrix(data = NA, ncol = 1, nrow = 20)
poll.sds   <- matrix(data = NA, ncol = 1, nrow = 20)

for (i in 1:pollsters) {
  poll.surveyed <- sample_n(precincts, size = N, replace = FALSE)
  poll.means[i] <- mean(poll.surveyed$black)
  poll.sds[i]   <-   sd(poll.surveyed$black)
}

## Note that we didn't save the 95% CI in each iteration. That's because
## we can create the CIs outside the loop much more quickly

ci.upper <- poll.means + 1.96*poll.sds/sqrt(N)
ci.lower <- poll.means - 1.96*poll.sds/sqrt(N)

## Let's plot the sampling distribution of the mean
ggplot() + geom_density(aes(x = poll.means)) 

## Now let's plot the confidence intervals
## We need to create a ruler that will denote the vertical locations of the CIs
ruler <- 1:pollsters

ggplot() + geom_point(aes(x = poll.means, y = ruler)) +
  geom_segment(aes(x = ci.lower, xend = ci.upper,
                   y = ruler, yend = ruler)) +
  geom_vline(xintercept = mean(precincts$black), 
             linetype = "dashed", color = "red")

df <- data.frame(poll.means, ruler)
ggplot(data = df, aes(x = poll.means, y = ruler)) + 
  geom_point()

# Check if confidence intervals contain the 'true' mean
ci_contains_mu <- ifelse(ci.upper > mean(precincts$black) &
                  ci.lower < mean(precincts$black), 
                TRUE, FALSE)

ggplot() + geom_point(aes(x = poll.means, y = ruler)) +
  geom_segment(aes(x = ci.lower, xend = ci.upper,
                   y = ruler, yend = ruler, color = ci_contains_mu)) +
  geom_vline(xintercept = mean(precincts$black), 
             linetype = "dashed", color = "red") 
  
## How many of the CIs contain the true mean?
mean(ci_contains_mu)

########################
## hypothesis testing ##
########################

## ok, let's say a fancy pollster thinks that there's
## no way that the turnout rate was more than 35%

## Let's use hypothesis testing to test his conclusion

## STEP 1: take a sample
## why? otherwise, you just know the truth
hypo.samp <- sample_n(precincts, 200, replace = FALSE)

## Step 2: take the sample mean (e.g., the estimator you are interested in)
##  and sample standard deviation

samp.mean <- mean(hypo.samp$turnout)
samp.sd <- sd(hypo.samp$turnout)

## Step 3: set up a test statistic
## it can be a lot of things, but generally people
## normalize the estimate to the standard normal and then
## check to see how unusual (or unusual) that test statistic would be

test.stat <- (samp.mean - .35)/(samp.sd/sqrt(N))

## step 4: calculate rejection regions

## so for alpha = 0.05 and 0.01:

qnorm(1 - .025)	## for alpha = 0.05 (why??)
qnorm(0 + .025)	## for alpha = 0.05 

qnorm(1 - .005)	## for alpha = 0.01 (why??)
qnorm(1 - .005)	## for alpha = 0.01 

## let's plot these against a standard normal to see what they look like

ruler <- seq(-6,6,.01)
plot(ruler, dnorm(ruler), type = "l", xlab = "standard normal")
abline(v = c(qnorm(1 - .025), qnorm(0 + .025)), col = "red")
abline(v = c(qnorm(1 - .005), qnorm(0 + .005)), col = "blue")
abline(v = test.stat, col = "green")
legend(x = "topleft", legend = c("alpha = .05", "alpha = .01", "test stat"), lwd = 2, col =
         c("red", "blue", "green"))

