library(tidyverse)

.22*.083


library(tidyverse)
set.seed(123)  # for reproducibility

n <- 10000

# Define correlation matrix
rho <- 0.5
Sigma <- matrix(c(1, rho, rho, 1), ncol = 2)

# Generate bivariate normal data
data <- MASS::mvrnorm(n = n, mu = c(0, 0), Sigma = Sigma)

# Convert to uniform(0,1) via the normal CDF
data_uniform <- pnorm(data)

# Extract vectors
p <- data_uniform[,1]
y <- data_uniform[,2]

sim.dat <- cbind.data.frame(p,y)
sim.dat <- sim.dat |> 
  mutate(y = if_else(y>=.5,1,0))

head(sim.dat)
cor(sim.dat$p, sim.dat$y)

result <- NA
for(i in 1:10000){
  result[i] <- mean(sample(sim.dat$y, 100, prob=sim.dat$p))
}
mean(result)
mean(sim.dat$y)

mean(result) - mean(sim.dat$y) 


cov(sim.dat$p, sim.dat$y)/mean(sim.dat$p)

national.targets <- rio::import("https://github.com/marctrussler/IIS-Data/raw/refs/heads/main/NationalWeightingTargets.Rds", trust=T)

sm <- rio::import("https://github.com/marctrussler/IIS-Data/raw/refs/heads/main/NBCAug25.Rds", trust=T)

sm |>
  mutate(approve.trump = if_else(trump_approval %in% c("Strongly approve","Somewhat approve"),T,F)) -> sm
mean(sm$approve.trump)


prop.table(table(sm$sex))

national.targets$sex

prop.table(table(sm$sex, sm$approve.trump),1)

sm$weight <- NA
sm$weight[sm$sex=="Male"] <- 0.4902613/0.473
sm$weight[sm$sex=="Female"] <- 0.5097387/0.527

head(sm)

weighted.mean(sm$sex=="Male", w = sm$weight)

mean(sm$approve.trump)
weighted.mean(sm$approve.trump, w=sm$weight)

prop.table(table(sm$educ))
national.targets$educ

sm |> 
  mutate(college = if_else(educ %in% c("College","Postgrad"),T,F)) -> sm
mean(sm$college)

national.targets$educ |> 
  filter(educ %in% c("Postgrad", "College")) |> 
  summarize(sum(Freq))


sm$weight <- NA
sm$weight[sm$college] <- .3304235/.545
sm$weight[!sm$college] <- (1-.3304235)/(1-.545)
head(sm)

weighted.mean(sm$college, w=sm$weight)

weighted.mean(sm$approve.trump, w=sm$weight)

sm |> 
  mutate(sex.college = paste(sex, college, sep="_")) -> sm

prop.table(table(sm$sex.college))

2*5*4*4*4

sm |> 
  group_by(sex, race, age4, educ,  vote.pres.2024) |> 
  summarise(n())

sm |> 
  mutate(college = if_else(college, "College","Non-College")) -> sm

sample.n <- table(sm$sex, sm$college)
sample.n

sm |> 
  group_by(sex, college) |> 
  summarise(n=n()) -> sample.n


#Also simplifying these weighting targets a bit to make the next stage clearer
national.targets$educ |> 
  mutate(educ= case_when(educ %in% c("Hs Or Less", "Some College") ~ "Non-College",
                         educ %in% c("College","Postgrad") ~ "College")) |> 
  group_by(educ) |> 
  summarize(Freq = sum(Freq)) |> 
  rename(college = educ,
         pop.prop = Freq) -> college
national.targets$sex |> 
  rename(pop.prop = Freq)-> sex

college
sex

#Gender Stage 1
sample.n |> 
  group_by(sex) |>
  #Calculate sample proportion for each group
  summarise(n = sum(n)) |> 
  mutate(sample.prop = n/sum(n)) |> 
  #Merge in the population proportion
  left_join(sex) |> 
  #Calculate cell weight as the population proportion divided by the sample proportion
  mutate(interim.weight = pop.prop/sample.prop) |>  
  select(sex, interim.weight) |> 
  #Merge that back to our original table of cell sizes
  right_join(sample.n) |> 
  #Calculate a weighted n, which is the original n multiplied by these weights
  mutate(weighted.n = n*interim.weight) -> weights

weights |> 
  group_by(sex) |> 
  summarise(n=sum(weighted.n)) |> 
  mutate(prop = n/sum(n))

weights |> 
  group_by(college) |> 
  summarise(n=sum(weighted.n)) |> 
  mutate(prop = n/sum(n))

#College Stage 1
#Start with the weights from previous stage
weights |> 
  group_by(college) |> 
  #Calculate the proportion of each group using weighted data from previous stage
  summarise(n = sum(weighted.n)) |> 
  mutate(sample.prop = n/sum(n)) |> 
  #Merge with population proportions
  left_join(college) |> 
  #Calculate cell weight by dividing population proportion by sample proportion
  mutate(stage.weight = pop.prop/sample.prop) |>  
  select(college, stage.weight) |> 
  #Merge back with starting data
  right_join(weights) |> 
  #Multiply the weight from this stage with the weight from the previous stage
  mutate(interim.weight = interim.weight*stage.weight,
         weighted.n = n*interim.weight) |> 
  select(-stage.weight)-> weights

weights

weights |> 
  group_by(sex) |> 
  summarise(n=sum(weighted.n)) |> 
  mutate(prop = n/sum(n))

weights |> 
  group_by(college) |> 
  summarise(n=sum(weighted.n)) |> 
  mutate(prop = n/sum(n))

#Sex Stage 2
weights |> 
  group_by(sex) |> 
  summarise(n = sum(weighted.n)) |> 
  mutate(sample.prop = n/sum(n)) |> 
  left_join(sex) |> 
  mutate(stage.weight = pop.prop/sample.prop) |>  
  select(sex, stage.weight) |> 
  right_join(weights) |> 
  mutate(interim.weight = interim.weight*stage.weight,
         weighted.n = n*interim.weight) |> 
  select(-stage.weight)-> weights


#College Stage 2

weights |> 
  group_by(college) |> 
  summarise(n = sum(weighted.n)) |> 
  mutate(sample.prop = n/sum(n)) |> 
  left_join(college) |> 
  mutate(stage.weight = pop.prop/sample.prop) |>  
  select(college, stage.weight) |> 
  right_join(weights) |> 
  mutate(interim.weight = interim.weight*stage.weight,
         weighted.n = n*interim.weight) |> 
  select(-stage.weight)-> weights

weights |> 
  group_by(sex) |> 
  summarise(n=sum(weighted.n)) |> 
  mutate(prop = n/sum(n))

weights |> 
  group_by(college) |> 
  summarise(n=sum(weighted.n)) |> 
  mutate(prop = n/sum(n))


sm |> 
  left_join(weights |>  rename(final.weight = interim.weight), join_by(sex, college)) -> sm

weighted.mean(sm$sex=="Female", w=sm$final.weight)

weighted.mean(sm$college=="College", w=sm$final.weight)


weighted.mean(sm$approve.trump, w=sm$final.weight)

#?pewmethods::rake_survey

sum(sm$final.weight)

weighting.targets <- national.targets[c("race","educ","age4","sex","vote.pres.2024")]

sm$rake.weights <- pewmethods::rake_survey(sm, weighting.targets)

#Saving a copy of this to use in subsequent lection
#sm |> 
#  select(-final.weight, -n, -weighted.n, -weight) |> 
#  rio::export("/Users/marctrussler/Documents/GitHub/IIS-Data/NBCAug25Weighted.Rds")


sm |> 
  mutate(educ = if_else(educ=="Hs Or Less", "HS Or Less",educ)) -> sm2
#pewmethods::rake_survey(sm2, weighting.targets)

sm2 <- sm
sm2$vote.pres.2024[2] <- NA
head(sm2)

sm2$weight <- pewmethods::rake_survey(sm2, weighting.targets)

sm2$weight[1]
sm2$weight[1000]

weighting.targets$race
#Correctly calculated weights
pewmethods::get_totals("race", sm, "rake.weights")
#Incorrectly calculated weights
pewmethods::get_totals("race",sm2,"weight")

summary(sm$rake.weights)

sm[sm$rake.weights==max(sm$rake.weights),]

sm[sm$rake.weights==min(sm$rake.weights),]

sm$trimmed.weights <- pewmethods::trim_weights(sm$rake.weights, lower_quantile = .01, upper_quantile=.99)
summary(sm$trimmed.weights)

pewmethods::get_totals("race", sm, "trimmed.weights")

sm$adj.trimmed.weights <- pewmethods::rake_survey(sm, weighting.targets, base_weight = "trimmed.weights")

summary(sm$rake.weights)
summary(sm$trimmed.weights)
summary(sm$adj.trimmed.weights)

pewmethods::get_totals("race", sm, "trimmed.weights")
pewmethods::get_totals("race", sm, "adj.trimmed.weights")

1 + (sd(sm$rake.weights)/mean(sm$rake.weights))^2
pewmethods::calculate_deff(sm$rake.weights)


1.96*sqrt(.25/468)