pop <- rbinom(1E6, 1, .5)
means <- NA

for(i in 1:10000){
  means[i] <- mean(sample(pop, 1000))
}
sd(means)

pop <- rbinom(1E8, 1, .5)
means <- NA

for(i in 1:10000){
  means[i] <- mean(sample(pop, 1000))
}
sd(means)

pop <- rbinom(1001, 1, .5)
means <- NA

for(i in 1:10000){
  means[i] <- mean(sample(pop, 1000))
}
sd(means)

pop <- 10^seq(3,7,0.1)

correction <- sqrt((pop-1000)/(pop-1))

plot(log10(pop), correction, axes=F, pch=16, 
     xlab="Population Size", 
     ylab="Finite Population Correction")
axis(side=1, at = 3:7, labels = 10^(3:7), )
axis(side=2, las=2)
abline(h=1, lty=2, col="firebrick")



set.seed(2105)
#Cluster sampling simulation
library(tidyverse)
# 40 classrooms of 25 people

classroom <- sort(rep(1:40,25))
student <- rep(1:25,40)

population <- cbind.data.frame(classroom, student)

#Randomly assign a latent level to each classroom
classroom <- 1:40
latent.level <- runif(40,-3.4,3.4)
latent <- cbind.data.frame(classroom, latent.level)

population |>  left_join(latent) -> population

#Generate student outcomes - Case 1, outcomes closely linked to latent levels
population$individual.prob <- population$latent.level + rnorm(nrow(population), mean=0, sd=2)
population$outcome <- if_else(pnorm(population$individual.prob)>.5,T,F)
population |> 
  group_by(classroom) |> 
  summarize(mean(outcome))


random.sample <- population[sample(1:nrow(population), 100),]
sort(unique(random.sample$classroom))

mean(random.sample$outcome)

#Calculate standard error
sqrt(.25/100)

#Estimate standard error through sampling: 

sample.mean <- NA

for(i in 1:1000){
  sample.mean[i] <- mean(sample(population$outcome, 100))
}
#Unbiased
mean(population$outcome)
mean(sample.mean)
#Standard Error (Approximately correct)
sd(sample.mean)

population |> 
  group_by(classroom) |> 
  summarize(mean(outcome))

#Do it once first:

sample <- population[population$classroom %in% sample(1:40, 4),]
table(sample$classroom)
mean(sample$outcome)

#Now simulate doing this many times: 

sample.means <- NA
for(i in 1:1000){
  sample <- population[population$classroom %in% sample(1:40, 4),]
  sample.means[i] <- mean(sample$outcome)
}



sd(sample.means)

plot(density(sample.means), main="Cluster Sampling Distribution")
abline(v=mean(sample.means), lty=2, col="darkblue", lwd=3)
abline(v=.49, lty=2, col="firebrick")

sd(sample.means)

plot(density(sample.means), main="Cluster Sampling Distribution")
abline(v=mean(sample.means), lty=2, col="darkblue", lwd=3)
abline(v=.49, lty=2, col="firebrick")
points(seq(0,1,.001), dnorm(seq(0,1,.001), mean=.49, sd=.05), col="orange", type="l")

population |> 
  group_by(classroom) |> 
  summarise(classroom.means = mean(outcome)) |> 
  summarise(var(classroom.means))

sqrt(.1/4)

set.seed(2105)
#Cluster sampling simulation
library(tidyverse)
# 40 classrooms of 25 people

classroom <- sort(rep(1:40,25))
student <- rep(1:25,40)

population2 <- cbind.data.frame(classroom, student)

#Randomly assign a latent level to each classroom
classroom <- 1:40
latent.level <- runif(40,-3.4,3.4)
latent <- cbind.data.frame(classroom, latent.level)

population2 |>  left_join(latent) -> population2

#Generate student outcomes - Case 1, outcomes closely linked to latent levels
population2$individual.prob <- population2$latent.level + rnorm(nrow(population2), mean=0, sd=10)
population2$outcome <- if_else(pnorm(population2$individual.prob)>.5,T,F)
population2 |> 
  group_by(classroom) |> 
  summarize(classroom.means = mean(outcome)) |> 
  summarize(var(classroom.means))


sqrt(.018/4)


sample.means <- NA
for(i in 1:1000){
  sample <- population2[population2$classroom %in% sample(1:40, 4),]
  sample.means[i] <- mean(sample$outcome)
}
sd(sample.means)


sqrt(.1/10)

sample.means <- NA
for(i in 1:1000){
sample <- population[population$classroom %in% sample(1:40, 10) & 
                        population$student %in% sample(1:25,10),]
  sample.means[i] <- mean(sample$outcome)
}
sd(sample.means)

#Creating fake data. Ignore this. 
library(tidyverse)
library(MASS)
set.seed(123)
n <- 100000
cor_matrix <- matrix(c(
  1.00,  0.20,  0.20, -0.10,  0.40,  # male
  0.20,  1.00,  0.30,  0.25,  0.50,  # white
  0.20,  0.30,  1.00,  0.35,  0.45,  # rich
 -0.10,  0.25,  0.35,  1.00, -0.30,  # college
  0.40,  0.50,  0.45, -0.30,  1.00   # vote.republican
), nrow = 5, byrow = TRUE)

# Simulate latent traits
latent <- mvrnorm(n = n, mu = rep(0, 5), Sigma = cor_matrix)

# Thresholds to get realistic marginal distributions:
# ~50% male, ~60% white, ~25% rich, ~55% college, ~40% vote.republican
male             <- latent[, 1] > 0
white            <- latent[, 2] > -0.25
rich             <- latent[, 3] > 0.67
college          <- latent[, 4] > -0.13
vote.republican  <- latent[, 5] > 0

detach(package:MASS)

# Combine into a data frame
pop <- data.frame(
  male = male,
  white = white,
  rich = rich,
  college = college,
  vote.republican = vote.republican
)

#Create strata

pop |> 
  mutate(race.education = case_when(white==F & college==F ~ "Non-white non-college",
                                    white==F & college==T ~ "Non-white college",
                                    white==T & college==F ~ "White non-college",
                                    white==T & college==T ~ "White college"),
         race.education = factor(race.education, levels = c("Non-white non-college",
                                                            "Non-white college",
                                                            "White non-college",
                                                            "White college"))) ->pop

pop |>  select(vote.republican, race.education) -> pop


head(pop)
mean(pop$vote.republican)

pop |> 
  group_by(race.education) |> 
  summarise(size = (n()/100000)*100,
            rep.percent = mean(vote.republican)*100)

sample.means.srs <- NA

for(i in 1:1000){
  sample.means.srs[i] <- mean(pop$vote.republican[sample(1:nrow(pop), 400)])
}

sd(sample.means.srs)

plot(density(sample.means.srs), main="Sampling Distribution")
abline(v=.501, lty=2, col="firebrick", lwd=2)


sample.means.strat1 <- NA

for(i in 1:1000){
  sample<- bind_rows(pop |> filter(race.education=="Non-white non-college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="Non-white college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="White non-college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="White college") |>  slice_sample(n=100))
  sample.means.strat1[i] <- mean(sample$vote.republican)
}
plot(density(sample.means.srs), main="Sampling Distribution")
points(density(sample.means.strat1), type="l", col="blue")
abline(v=.501, lty=2, col="firebrick", lwd=2)

sd(sample.means.strat1)


mean(sample.means.strat1)

pop |> 
  filter(race.education=="Non-white college") |> 
  summarise(n())

pop |> 
  filter(race.education=="White college") |> 
  summarise(n())

  sample<- bind_rows(pop |> filter(race.education=="Non-white non-college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="Non-white college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="White non-college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="White college") |>  slice_sample(n=100)) 

  #Assign weights
  sample |> 
    mutate(weight = case_when(race.education=="Non-white non-college" ~ 1/(100/21750),
                              race.education=="Non-white college" ~ 1/(100/18392),
                              race.education=="White non-college" ~1/(100/23133),
                              race.education=="White college" ~ 1/(100/36725))) -> sample

weighted.mean(sample$vote.republican, w=sample$weight)

sample.means.strat2 <- NA

for(i in 1:1000){
    sample<- bind_rows(pop |> filter(race.education=="Non-white non-college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="Non-white college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="White non-college") |>  slice_sample(n=100),
                                          pop |> filter(race.education=="White college") |>  slice_sample(n=100)) 

  #Assign weights
  sample |> 
    mutate(weight = case_when(race.education=="Non-white non-college" ~ 1/(100/21750),
                              race.education=="Non-white college" ~ 1/(100/18392),
                              race.education=="White non-college" ~1/(100/23133),
                              race.education=="White college" ~ 1/(100/36725))) -> sample

sample.means.strat2[i] <- weighted.mean(sample$vote.republican, w=sample$weight)
}

plot(density(sample.means.srs), main="Sampling Distribution")
points(density(sample.means.strat1), type="l", col="blue")
points(density(sample.means.strat2), type="l", col="orange")
abline(v=.501, lty=2, col="firebrick", lwd=2)

sd(sample.means.strat2)




#Sample sizes:
round(400*.21750)
round(400*.18392)
round(400*.23133)
round(400*.36725)

sample.means.strat3 <- NA

for(i in 1:1000){
  sample<- bind_rows(pop |> filter(race.education=="Non-white non-college") |>  slice_sample(n=87),
                                          pop |> filter(race.education=="Non-white college") |>  slice_sample(n=74),
                                          pop |> filter(race.education=="White non-college") |>  slice_sample(n=93),
                                          pop |> filter(race.education=="White college") |>  slice_sample(n=147))
  sample.means.strat3[i] <- mean(sample$vote.republican)
}
plot(density(sample.means.srs), main="Sampling Distribution")
points(density(sample.means.strat1), type="l", col="blue")
points(density(sample.means.strat2), type="l", col="orange")
points(density(sample.means.strat3), type="l", col="magenta")

abline(v=.501, lty=2, col="firebrick", lwd=2)

sd(sample.means.strat3)


#Creating fake data. Ignore this. 
library(MASS)
set.seed(123)
n <- 100000
cor_matrix <- matrix(c(
  1.00,  0.20,  0.20, -0.10,  0.40,  # male
  0.20,  1.00,  0.30,  0.25,  0.50,  # white
  0.20,  0.30,  1.00,  0.35,  0.45,  # rich
  -0.10,  0.25,  0.35,  1.00, -0.30,  # college
  0.40,  0.50,  0.45, -0.30,  1.00   # vote.republican
), nrow = 5, byrow = TRUE)

# Simulate latent traits
latent <- mvrnorm(n = n, mu = rep(0, 5), Sigma = cor_matrix)

# Thresholds to get realistic marginal distributions:
# ~50% male, ~60% white, ~25% rich, ~55% college, ~40% vote.republican
male             <- latent[, 1] > 0
white            <- latent[, 2] > -0.25
rich             <- latent[, 3] > 0.67
college          <- latent[, 4] > -0.13
vote.republican  <- latent[, 5] > 0

detach(package:MASS)

# Combine into a data frame
pop2 <- data.frame(
  male = male,
  white = white,
  rich = rich,
  college = college,
  vote.republican = vote.republican,
  latent.republican = latent[,5]
)
pop2 |>  select(vote.republican, latent.republican) -> pop2



head(pop2)

#Create strata

pop2 |> 
  mutate(republican.strata =ntile(latent.republican, 10)) -> pop2

head(pop2)

sample.means.strat4 <- NA
for(i in 1:1000){
samp <- list()
for(j in 1:10){
samp[[j]] <- pop2 |> filter(republican.strata==j) |> slice_sample(n=40)  
}

samp <- bind_rows(samp)


sample.means.strat4[i] <- mean(samp$vote.republican)
}

mean(sample.means.strat4)
sd(sample.means.strat4)

#Neyman allocation

#1. For each strata calculate the proportion and variance:

pop |> 
  group_by(race.education) |> 
  summarize(n = n(), 
            var = var(vote.republican)) |> 
  mutate(prop = n/sum(n)) -> ney


#2. Calculate the neyman allocation percentage by multiplying the proportion and variance, and
#then normalizing to 1. 

ney |> 
  mutate(product = var*prop,
         allocation = product/sum(product)) -> ney

#3. Use the allocations to get an $n$ for each group
ney |> 
  mutate(sampling.n = round(allocation*400)) -> ney

#4. Calculate inverse probability weights

ney |> 
  mutate(weight = 1/(sampling.n/n)) -> ney

#4. Merge this back in to population

pop |> 
  left_join(ney, join_by(race.education)) -> pop

#5. Sample according to these allocations
   sample<- bind_rows(pop |> filter(race.education=="Non-white non-college") |>  slice_sample(n=ney$sampling.n[1]),
                                          pop |> filter(race.education=="Non-white college") |>  slice_sample(n=ney$sampling.n[2]),
                                          pop |> filter(race.education=="White non-college") |>  slice_sample(n=ney$sampling.n[3]),
                                          pop |> filter(race.education=="White college") |>  slice_sample(n=ney$sampling.n[4])) 

#.6 Calculate weighted mean
 weighted.mean(sample$vote.republican, w=sample$weight)

#7. Iterate this many times to see sampling distribution
 
neyman.means <- NA
 for(i in 1:1000){
 sample<- bind_rows(pop |> filter(race.education=="Non-white non-college") |>  slice_sample(n=ney$sampling.n[1]),
                                          pop |> filter(race.education=="Non-white college") |>  slice_sample(n=ney$sampling.n[2]),
                                          pop |> filter(race.education=="White non-college") |>  slice_sample(n=ney$sampling.n[3]),
                                          pop |> filter(race.education=="White college") |>  slice_sample(n=ney$sampling.n[4])) 

#.6 Calculate weighted mean
 neyman.means[i] <- weighted.mean(sample$vote.republican, w=sample$weight)
}
 
mean(neyman.means)
sd(neyman.means)

#I guess it's a bit smaller? I don't know this seems like a pretty big hassle!


#States and population:
state_pops <- data.frame(
  state = c("AL","AK","AZ","AR","CA","CO","CT","DE","FL","GA",
            "HI","ID","IL","IN","IA","KS","KY","LA","ME","MD",
            "MA","MI","MN","MS","MO","MT","NE","NV","NH","NJ",
            "NM","NY","NC","ND","OH","OK","OR","PA","RI","SC",
            "SD","TN","TX","UT","VT","VA","WA","WV","WI","WY"),
  pop = c(
    5157699, 740133, 7582384, 3088354, 39431263,
    5939456, 3605944, 1013876, 22244823, 10912876,
    1455271, 1979278, 12518071, 6833037, 3197307,
    2937880, 4536153, 4625500, 1356458, 6201046,
    6941634, 10034113, 5709752, 2930210, 6209670,
    1122667, 1992876, 3358392, 1408345, 9267130,
    2121027, 20058984, 10809634, 779094, 11777593,
    4009574, 4275583, 12932753, 1067582, 5339274,
    895376, 7064187, 30381457, 3497197, 647230, 8819590,
    7785786, 1793716, 5893718, 576851
  )
)

state_pops$prob <- state_pops$pop / sum(state_pops$pop)

cbind(state_pops$state, round(state_pops$prob*1000))