library(tidyverse)

set.seed(456)  # for reproducibility

n <- 10000

# Define correlation matrix
rho <- 0.2
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)

#Select people to respond based on probability of response

select <- sample(1:nrow(sim.dat), 1000, prob=sim.dat$p)
sim.dat |> 
  mutate(r = if_else(row_number() %in% select, T, F)) -> sim.dat


#Actual Error
mean(sim.dat$y[sim.dat$r]) - mean(sim.dat$y)

#Meng Equation
cor(sim.dat$r, sim.dat$y)*sqrt((10000-1000)/1000)*sd(sim.dat$y)

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

small.sample.correction <- sqrt((pop-1000)/(pop-1))
NINR.correction <- sqrt((pop-1000)/(1000))

par(mfrow=c(1,2))
plot(log10(pop), small.sample.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")

plot(log10(pop), NINR.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)

set.seed(124)  # for reproducibility

n <- 1000000

# Define correlation matrix
rho <- 0.2
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)

#####
#Nonprobability sample: randomly select 1000 peopel to contact
select <- sample(1:nrow(sim.dat), 1000, prob=sim.dat$p)
sim.dat |> 
  mutate(r = if_else(row_number() %in% select, T, F)) -> sim.dat
#Actual Error
mean(sim.dat$y[sim.dat$r]) - mean(sim.dat$y)
#Meng Equation for nonprob sample
cor(sim.dat$r, sim.dat$y)*sqrt((1000000-1000)/1000)*sd(sim.dat$y)

#####
#Probability sample: randomly select 5000 people to contact
select <- sample(1:nrow(sim.dat), 5000)
sim.dat |> 
  mutate(c = if_else(row_number() %in% select, T, F)) |> 
  filter(c) ->contact.dat
#Now simulate 1000 of these people responding based on p
select <- sample(1:nrow(contact.dat), 1000, prob=contact.dat$p)
contact.dat |> 
  mutate(r = if_else(row_number() %in% select, T, F)) -> contact.dat
#Actual Error
mean(contact.dat$y[contact.dat$r]) - mean(contact.dat$y)
#Meng Equation for nonprob sample
p <- 1000/5000
cor(contact.dat$r, contact.dat$y)*sqrt((1-p)/p)*sd(contact.dat$y)

denom <- ((.005^2)*(1e+06 - 1500))+1
1500/denom