#Calculate t score for alpha=.05
t <- qt(.05,df=7)*-1

#SE
se <- 2/sqrt(8)

#Confidence interval
c(4-t*se,4+t*se)

set.seed(19102)
pa.poll <- rbinom(1000,1, .525)

mean(pa.poll)
var(pa.poll)
sd(pa.poll)

qt(.025, df=999)
se <- sd(pa.poll)/sqrt(1000)

#Confidence Interval
c(mean(pa.poll)-1.96*se ,mean(pa.poll)+1.96*se)


eval <- seq(.45,.55, .0001)
plot(eval, dnorm(eval, mean=.5, sd=sqrt(.249/1000)), type="l",
     main="Sampling Distribution of Fetterman Support: Tied Race")
abline(v=.5,lty=2)
abline(v=mean(pa.poll), lty=2, col="firebrick")

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

eval <- seq(.45,.55, .0001)
plot(eval, dnorm(eval, mean=.5, sd=sqrt(.249/1000)), type="l",
     main="Sampling Distribution of Fetterman Support: Tied Race")
abline(v=.5,lty=2)
abline(v=mean(pa.poll), lty=2, col="firebrick")
normal.shader(.45,.55,.5,sqrt(.249/1000), .53, greater = T, color = "firebrick")

pnorm(mean(pa.poll), mean=.5, sd=sqrt(.249/1000), lower.tail = F)


results <- rnorm(10000, mean = .5, sd = sqrt(.249/1000))
mean(results>.53)


#What value puts 5% in the right tail?
qnorm(.05, .5,sqrt(.249/1000), lower.tail=F)

eval <- seq(.40,.60, .0001)
plot(eval, dnorm(eval, mean=.5, sd=sqrt(.249/1000)), type="l",
     main="Sampling Distribution of Fetterman Support: Tied Race", xlab="Support",
     ylab="Density")
abline(v=.5,lty=2)
normal.shader(.4,.6, .5,sqrt(.249/1000), .526, greater=T, color="dodgerblue" )
abline(v=mean(pa.poll), lty=3, col="firebrick")


eval <- seq(.45,.55, .0001)
plot(eval, dnorm(eval, mean=.5, sd=sqrt(.249/1000)), type="l",
     main="Sampling Distribution of Fetterman Support: Tied Race")
abline(v=.5,lty=2)
normal.shader(.45,.55, .5,sqrt(.249/1000),qnorm(.025, mean=.5, sd=sqrt(.249/1000)), col="dodgerblue"  )
normal.shader(.45,.55, .5,sqrt(.249/1000),qnorm(.025, mean=.5, sd=sqrt(.249/1000), lower.tail=F), col="dodgerblue", greater=T)
abline(v=mean(pa.poll), lty=2, col="firebrick")

eval <- seq(.40,.60, .0001)
plot(eval, dnorm(eval, mean=.5, sd=sqrt(.249/1000)), type="l",
     main="Sampling Distribution of Fetterman Support: Tied Race", xlab="Support",
     ylab="Density")
abline(v=.5,lty=2)
normal.shader(.4,.6, .5, sqrt(.249/1000), .53, greater=T, col="firebrick")

pnorm(mean(pa.poll), mean=.5, sd=sqrt(.249/1000), lower.tail=F)

eval <- seq(.40,.60, .0001)
plot(eval, dnorm(eval, mean=.5, sd=sqrt(.249/1000)), type="l",
     main="Sampling Distribution of Fetterman Support: Tied Race", xlab="Support",
     ylab="Density")
abline(v=.5,lty=2)
normal.shader(.4,.6, .5, sqrt(.249/1000), .53, greater=T, col="firebrick")
#Need to look at area to the left of .47, which is our sample mean flipped to the
#other side of the null hypothesis:
normal.shader(.4,.6, .5, sqrt(.249/1000), .47, greater=F, col="firebrick")


pnorm(mean(pa.poll), mean=.5, sd=sqrt(.249/1000), lower.tail=F) + pnorm(.5 - (mean(pa.poll)-.5), mean=.5, sd=sqrt(.249/1000), lower.tail=T)

pnorm(mean(pa.poll), mean=.5, sd=sqrt(.249/1000), lower.tail=F)*2

#Construct a 95\% confidence interval
mean(pa.poll)
t.value <- qt(.025, df=999)*-1

ci <- c(mean(pa.poll)-t.value*sqrt(.249/1000), mean(pa.poll)+t.value*sqrt(.249/1000) )
library(dplyr)
between(.5,ci[1], ci[2] )
#Null is within the CI

#Hypothesis test:
p <- pnorm(mean(pa.poll), mean=.5, sd=sqrt(.249/1000), lower.tail=F)*2
p>.05

reject.null.ci <- rep(NA, 10000)
reject.null.hyp <- rep(NA, 10000)

for(i in 1:10000){
  samp <- rbinom(1000, 1, .525)
  xbar <- mean(samp)
  se <- sd(samp)/sqrt(1000)

#ci
  ci <- c(xbar-t.value*se, xbar+t.value*se)
  reject.null.ci[i] <- between(.5, ci[1], ci[2])
#hypothesis
#Use the absolute deviation from the null so the formula works whether
#xbar happens to fall above or below .5.
  p <- 2*pnorm(abs(xbar - .5)/se, lower.tail=F)
  reject.null.hyp[i] <- p>.05
}

table(reject.null.ci, reject.null.hyp)
#Same!

library(rio)
sm.az <- import("https://github.com/marctrussler/IIS-Data/raw/main/AZFinalWeeks.csv")
head(sm.az)

table(sm.az$senate.topline)

sm.az$DemSenateVote <- NA
sm.az$DemSenateVote[sm.az$senate.topline=="Democrat"]<- 1
sm.az$DemSenateVote[sm.az$senate.topline=="Republican"]<- 0
mean(sm.az$DemSenateVote,na.rm=T)
table(sm.az$DemSenateVote)
sum(table(sm.az$DemSenateVote))

mean(sm.az$DemSenateVote,na.rm=T)

s <- sd(sm.az$DemSenateVote,na.rm=T)

eval <- seq(.40,.60, .0001)
plot(eval, dnorm(eval, mean=.5, sd=s/sqrt(2613)), type="l",
     main="Sampling Distribution of Kelly Support: Tied Race", xlab="Support",
     ylab="Density")
abline(v=.5,lty=2)
normal.shader(.4,.6, .5, s/sqrt(2613), .5 - (.524-.5), greater=F, color="firebrick")
normal.shader(.4,.6, .5, s/sqrt(2613), .524, greater=T, color="firebrick")

2*pnorm(abs(.524 - .5)/(s/sqrt(2613)), lower.tail=F)

t <- (.524-.5)/(s/sqrt(2613))

2*pt(abs(t), df=2612, lower.tail=F)

set.seed(19104)
samp <- rbinom(1000,1, prob=.5)

mean(samp)

se <- sd(samp)/sqrt(1000)

eval <- seq(.40,.60, .0001)
rejection.region <- qnorm(.025, mean=.5, sd=se, lower.tail=F)

plot(eval, dnorm(eval, mean=.5, sd=se), type="l",
     main="Sampling Distribution under Null Hypothesis", xlab="Support",
     ylab="Density")
abline(v=.5,lty=2)
normal.shader(.4,.6, .5, se , rejection.region, greater=T, color="firebrick")
normal.shader(.4,.6, .5, se,.5 - (rejection.region-.5) , greater=F, color="firebrick")



#Use the absolute deviation from the null so we get the right two-tailed
#p-value regardless of which side of .5 our sample mean fell on.
2*pnorm(abs(mean(samp) - .5)/se, lower.tail=F)

p <- NA
for(i in 1:10000){
samp <- rbinom(1000,1, prob=.5)
#Shortcut with the t.test function:
p[i] <- t.test(samp, mu=.5)$p.value
}

mean(p<.05)


p <- NA
for(i in 1:10000){
samp <- rbinom(1000,1, prob=.5)
#Shortcut with the t.test function:
p[i] <- t.test(samp, mu=.5)$p.value
}

mean(p<.1)


#We want to know the probability of getting at least one head
#Which is 1 minus the probability of getting no heads
1-dbinom(0, 1, .05)
1-dbinom(0, 2, .05)


plot(1:100, 1-dbinom(0, 1:100, .05), type="l", pch=16, 
     main="Multiple Comparisons",
     xlab="Number of Hypothesis Tests",
     ylab = "Probability of at least one false positive")

m <- 10
coin <- c(0,1)
any.sig <- NA

for(j in 1:1000){
  sig <- NA
  for(i in 1:m){
    flips <- sample(coin, 100, replace=T)
    sig[i] <- t.test(flips, mu=.5)$p.value <.05
  }
  any.sig[j] <- any(sig)
}
mean(any.sig)

m <- 10
coin <- c(0,1)
any.sig <- NA

for(j in 1:1000){
  sig <- NA
  for(i in 1:m){
    flips <- sample(coin, 100, replace=T)
    sig[i] <- t.test(flips, mu=.5)$p.value <(.05/m)
  }
  any.sig[j] <- any(sig)
}
mean(any.sig)

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 using the absolute t so this works whether xbar landed
#above or below .5.
p <- 2*pt(abs(t.val), df=n-1, lower.tail=F)
false.negative[i] <- p>=alpha
}
mean(false.negative)

library(rio)
sm.az <- import("https://github.com/marctrussler/IIS-Data/raw/main/AZFinalWeeks.csv")
head(sm.az)

table(sm.az$biden.approval)
sm.az$approve.biden[sm.az$biden.approval %in% c("Strongly approve","Somewhat approve")] <- 1
sm.az$approve.biden[sm.az$biden.approval %in% c("Strongly disapprove","Somewhat disapprove")] <- 0
table(sm.az$biden.approval, sm.az$approve.biden)
#
mean(sm.az$approve.biden, na.rm=T)
sd(sm.az$approve.biden,na.rm=T)
sum(table(sm.az$approve.biden))

mean(sm.az$approve.biden,na.rm=T)

se <- sd(sm.az$approve.biden,na.rm=T)/sqrt(3024)

t.score <- (mean(sm.az$approve.biden,na.rm=T)- .38)/se

2*pt(abs(t.score), df=3023, lower.tail=F)

#Default is null=0
t.test(sm.az$approve.biden)

t.test(sm.az$approve.biden, mu=.38)

#One-tailed: Biden's approval is *greater than* .38
t.test(sm.az$approve.biden, mu=.38, alternative="greater")