source("https://raw.githubusercontent.com/marctrussler/IIS-Data/main/NormalShader.R")
library(MASS)

set.seed(19104)
false.negative <- rep(NA,10000)

n <- 2613
alpha <- .05

for(i in 1:10000){
samp <- rbinom(n,1, .524)
s <- sd(samp)
xbar <- mean(samp)
se <- s/sqrt(n)
t.val <- (xbar-.5)/se
#Two-tailed p-value via the absolute t so this works whether xbar happens
#to land above or below .5 (see the abs() explanation in the HypothesisTests
#chapter's CI/hypothesis equivalence loop).
p <- 2*pt(abs(t.val), df=n-1, lower.tail=F)
false.negative[i] <- p>=alpha
}
mean(false.negative)

x <- seq(.46,.57, .001)
n <- 2613
mu  <- .524
se <- sqrt((mu*(1-mu))/n)
alpha <-.05

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
legend("topleft", "True Sampling Distribution", lty=1, col="dodgerblue")

no <- .5

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ), greater=F, col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=F ), greater=T, col="firebrick")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))


plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")

#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )
#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)

n <- 2613
mu  <- .53
se <- sqrt((mu*(1-mu))/n)
no  <- .5
alpha <- .05

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")


#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )

#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)

n <- 2613
mu  <- .524
se <- sqrt((mu*(1-mu))/n)
no  <- .5
alpha <- .1

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")


#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )

#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)

n <- 5000
mu  <- .524
se <- sqrt((mu*(1-mu))/n)
no  <- .5
alpha <- .05

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")


#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )

#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)

n <- 10000
mu  <- .524
se <- sqrt((mu*(1-mu))/n)
no  <- .5
alpha <- .05

plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")


#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )

#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)

#Input
n <- 15
hA <- .55
hO <- .38
alpha <- .05

#SE of alternative
sd <- sqrt(hA*(1-hA))
seA <- sd/sqrt(n)

#SE of Null
sd <- sqrt(hO*(1-hO))
seO <- sd/sqrt(n)

#Calculate rejection zone for alpha
z.score <- qnorm(alpha/2)*-1
bounds <- c(hO+z.score*seO, hO-z.score*seO)

eval <- seq(.1,.7,.0001)
plot(eval, dnorm(eval, mean=hA, sd=seA), type="l", col="darkblue", main="What is Beta?")
legend("topleft","True Sampling Distribution", lty=1, col="darkblue")
points(eval, dnorm(eval, mean=hO, sd=seO), type="l", col="firebrick")
normal.shader(.1,.7, hO, seO, bounds[2], bounds[1], between=T, col="firebrick")
normal.shader(.1,.7, hA, seA, bounds[2], bounds[1], between=T, col="darkblue")
legend("topleft",c("True Sampling Distribution","Null Hyp. Sampling Distribution") , lty=c(1,1), col=c("darkblue", "firebrick"))

#Beta calculation:
pnorm(bounds[1], mean=hA, sd=seA) - pnorm(bounds[2], mean=hA, sd=seA) 

#Input
n <- 1000000
hA <- .385
hO <- .38
alpha <- .05

#SE of alternative
sd <- sqrt(hA*(1-hA))
seA <- sd/sqrt(n)

#SE of Null
sd <- sqrt(hO*(1-hO))
seO <- sd/sqrt(n)

#Calculate rejection zone for alpha
z.score <- qnorm(alpha/2)*-1
bounds <- c(hO+z.score*seO, hO-z.score*seO)

eval <- seq(.34,.44,.0001)
plot(eval, dnorm(eval, mean=hA, sd=seA), type="l", col="darkblue", main="What is Beta?")
legend("topleft","True Sampling Distribution", lty=1, col="darkblue")
points(eval, dnorm(eval, mean=hO, sd=seO), type="l", col="firebrick")
abline(v=bounds, lty=2, col="firebrick")
legend("topleft",c("True Sampling Distribution","Null Hyp. Sampling Distribution") , lty=c(1,1), col=c("darkblue", "firebrick"))

#Beta calculation:
pnorm(bounds[1], mean=hA, sd=seA) - pnorm(bounds[2], mean=hA, sd=seA) 

set.seed(19103)
sigma<-rbind(c(1,.2), c(.2,1))
d <- as.data.frame(mvrnorm(n=500, mu=c(0,0), Sigma=sigma))
names(d) <- c("x","y")

plot(d$x, d$y )
correlation <- cor(d$x, d$y)
correlation

z <- d$x
z[5] <- NA
cor(z, d$y)

cor(z, d$y, use="pairwise.complete")

names(correlation)

cor.samp.dist <- rep(NA, 10000)

for(i in 1:10000){
 sigma<-rbind(c(1,.2), c(.2,1))
 samp <- as.data.frame(mvrnorm(n=500, mu=c(0,0), Sigma=sigma))
 names(samp) <- c("x","y")
 cor.samp.dist[i] <- cor(samp$x, samp$y)
}
plot(density(cor.samp.dist))
abline(v=0, lty=2)
sd(cor.samp.dist)

bs.samp <- d[sample(1:nrow(d), replace=T),]
cor(d$x, d$y)
cor(bs.samp$x, bs.samp$y)
#Not the same!


bs.samp <- d[sample(1:nrow(d), replace=F),]
cor(d$x, d$y)
cor(bs.samp$x, bs.samp$y)
#The same!

bs.samp.dist <- rep(NA, 10000)
for(i in 1:10000){
  bs.samp <- d[sample(1:nrow(d), replace=T),]
  bs.samp.dist[i] <- cor(bs.samp$x, bs.samp$y)
}

plot(density(cor.samp.dist), col="darkblue")
points(density(bs.samp.dist), col="firebrick", type="l")
abline(v=0, lty=2)
sd(cor.samp.dist)
sd(bs.samp.dist)

quantile(bs.samp.dist, .025)
quantile(bs.samp.dist, .975)

z.score <- cor(d$x, d$y)/sd(bs.samp.dist)
p <- pnorm(z.score, lower.tail=F)*2
p

set.seed(19104)
prime.sample <- rnorm(1000, mean=0, sd=5)
mean(prime.sample)
sd(prime.sample)

se.calc <- sd(prime.sample)/sqrt(1000)

eval <- seq(-1,1, .0001)
plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")

se.real <- 5/sqrt(1000)
plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")
points(eval, dnorm(eval, mean=0, sd=se.real), col="darkblue", type="l")
legend("topleft", c("Calculated SampDist", "Real SampDist"), lty=c(1,1), col=c("firebrick","darkblue"))



samp.dist <- rep(NA, 10000)

for(i in 1:10000){
  samp.dist[i] <- mean(rnorm(1000, mean=0, sd=5))
}

plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")
points(eval, dnorm(eval, mean=0, sd=se.real), col="darkblue", type="l")
points(density(samp.dist), col="forestgreen", type="l")
legend("topleft", c("Calculated SampDist", "Real SampDist", "Empirical SampDist"), lty=c(1,1,1), col=c("firebrick","darkblue", "forestgreen"))




par(mfrow=c(3,3))
eval.pop <- seq(-20,20,.01)

for(i in 1:9){
  samp <- rnorm(1000, mean=0, sd=5)
  plot(eval.pop, dnorm(eval.pop, mean=0, sd=5), lwd=2, lty=1, type="l",
       main=paste("Mean = ", round(mean(samp),2)))
  points(density(samp), col="forestgreen", type="l")
}


par(mfrow=c(3,3))
eval.pop <- seq(-20,20,.01)

for(i in 1:9){
  samp <- prime.sample[sample(1:length(prime.sample), length(prime.sample), replace=T)]
  plot(eval.pop, dnorm(eval.pop, mean=0, sd=5), lwd=2, lty=1, type="l",
       main=paste("Mean = ", round(mean(samp),2)))
  points(density(samp), col="magenta", type="l")
}


bs.samp.dist <- rep(NA,10000)

for(i in 1:10000){
  bs.samp.dist[i] <- mean(prime.sample[sample(1:length(prime.sample), length(prime.sample), replace=T)])
}

bs.samp.dist <- bs.samp.dist - mean(bs.samp.dist)

plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")
points(eval, dnorm(eval, mean=0, sd=se.real), col="darkblue", type="l")
points(density(samp.dist), col="forestgreen", type="l")
points(density(bs.samp.dist), col="magenta", type="l")
legend("topleft", c("Calculated SampDist", "Real SampDist", "Empirical SampDist", "BS SampDist"), lty=c(1,1,1,1), col=c("firebrick", "darkblue","forestgreen","magenta"))

library(rio)
anes <- import("https://github.com/marctrussler/IIS-Data/raw/main/ANESFinalProjectData.csv")


#5 and 6 are Sanders and Warren
anes$vote.centrist.dem[anes$V201021 %in% c(1,2,3,4)] <- 1
anes$vote.centrist.dem[anes$V201021 %in% c(5,6)] <- 0
table(anes$vote.centrist.dem)

#1 is very liberal and 7 is very conservative.
anes$ideology <- anes$V201200
anes$ideology[anes$ideology %in% c(-9,-8,99)] <- NA
table(anes$ideology)

table(anes$vote.centrist.dem, anes$ideology)

cor(anes$ideology, anes$vote.centrist.dem, use="pairwise.complete")

anes$white.v.black[anes$V201549x==1] <- 1
anes$white.v.black[anes$V201549x==2] <- 0
table(anes$white.v.black)

#Among White Americans
cor(anes$ideology[anes$white.v.black==1], anes$vote.centrist.dem[anes$white.v.black==1], use="pairwise.complete")
#Among Black Americans
cor(anes$ideology[anes$white.v.black==0], anes$vote.centrist.dem[anes$white.v.black==0], use="pairwise.complete")
#Difference 
cor.delta <- cor(anes$ideology[anes$white.v.black==1], anes$vote.centrist.dem[anes$white.v.black==1], use="pairwise.complete") - cor(anes$ideology[anes$white.v.black==0], anes$vote.centrist.dem[anes$white.v.black==0], use="pairwise.complete")
cor.delta


boot.data <- anes[c("ideology","vote.centrist.dem","white.v.black")]
boot.data <- boot.data[complete.cases(boot.data),]

bs.cor.delta <- rep(NA,10000)
for(i in 1:10000){
  #Generate the bootstrap dataset
   bs.data <- boot.data[sample(1:nrow(boot.data), nrow(boot.data), replace=T),]
  #Calculat the same estimate as I calculated in the prime dataset, saving the result each time.
   bs.cor.delta[i] <- cor(bs.data$ideology[bs.data$white.v.black==1],bs.data$vote.centrist.dem[bs.data$white.v.black==1],use="pairwise.complete") - cor(bs.data$ideology[bs.data$white.v.black==0],bs.data$vote.centrist.dem[bs.data$white.v.black==0], use="pairwise.complete")
}


plot(density(bs.cor.delta), main="Bootstrap Sampling Distribution")
abline(v=0, lty=2)

#Our original estimate
cor.delta
#CI
quantile(bs.cor.delta, .025)
quantile(bs.cor.delta, .975)

se <- sd(bs.cor.delta)

#For my own sanity, make sure that the sample distribution is normally distributed:
eval <- seq(-.1, .5, .0001)
plot(density(bs.cor.delta), main="Bootstrap Sampling Distribution")
points(eval, dnorm(eval, mean=cor.delta, sd=se), col="firebrick", type="l")
abline(v=0, lty=2)

#For null of 0, how many SEs is our test statistic away?

z.score <- (cor.delta-0)/se

#Evalualte under the standard normal

pnorm(z.score, lower.tail=F)*2