R Inference & Modeling: The Central Limit Theorem in Practice

The Central Limit Theorem in Practice

• Because X is the sum of random draws divided by a constant, the distribution of X is approximately normal.
• We can convert X to a standard normal random variable Z:
  Z= (x-E(x))/SE(x)
  
• The probability that X is within .01 of the actual value of p is:
  Pr(Z<=.01/sqrt(p(1-p)/N))-Pr(Z<=-.01/sqrt(p(1-p)/N))
• The Central Limit Theorem (CLT) still works if X is use in place of p. This is called a plug-in estimate. Hats over values denote estimates. Therefore:
  SE(X) = sqrt(X(1-X)/N)
• Using the CLT, the probability that X is within .01 of the actual value of p is:
  Pr(Z<=.01/sqrt(X(1-X)/N))-Pr(Z<=-.01/sqrt(X(1-X)/N))

Code: Computing the probability of X being within .01 of p

X_hat <- 0.48
se <- sqrt(X_hat*(1-X_hat)/25)
pnorm(0.01/se) - pnorm(-0.01/se)

Margin of Error:

• The margin of error is defined as 2 times the standard error of the estimate X.
• There is about a 95% chance that X will be within two standard errors of the actual parameter p.

Monte Carlo Simulation of CLT

• We can run Monte Carlo simulations to compare with theoretical results assuming a value of p.
• In practice, p, is unknown. We can corroborate theoretical results by running Monte Carlo simulations with one or several values of p.
• One practical choice for p when modeling is X, the observed value of X in a sample.

Code: Monte Carlo Simulation using a set value of p

p <- 0.45    # unknown p to estimate
N <- 1000

# simulate one poll of size N and determine x_hat
x <- sample(c(0,1), size = N, replace = TRUE, prob = c(1-p, p))
x_hat <- mean(x)

# simulate B polls of size N and determine average x_hat
B <- 10000    # number of replicates
N <- 1000    # sample size per replicate
x_hat <- replicate(B, {
    x <- sample(c(0,1), size = N, replace = TRUE, prob = c(1-p, p))
    mean(x)
})

Code: Histogram and QQ-Plot of Monte Carlo Simulation

library(tidyverse)
library(gridExtra)
p1 <- data.frame(x_hat = x_hat) %>%
    ggplot(aes(x_hat)) +
    geom_histogram(binwidth = 0.005, color = "black")
p2 <- data.frame(x_hat = x_hat) %>%
    ggplot(aes(sample = x_hat)) +
    stat_qq(dparams = list(mean = mean(x_hat), sd = sd(x_hat))) +
    geom_abline() +
    ylab("X_hat") +
    xlab("Theoretical normal")
grid.arrange(p1, p2, nrow=1)

CLT Spread

• The spread between two outcomes with probabilities p and 1-p is 2p-1.
• The expected value of the spread is 2X-1
• The standard error of the spread is 2SE(X)
• The margin of error of the spread is 2 times the margin of error of X.

CLT Bias: Example Large Poll

• An extremely large poll would theoretically be able to predict election results almost perfectly.
• These sample sizes are not practical. In addition to cost concerns, polling doesn't reach everyone in the population (eventual voters) with equal probability, and it also may include data from outside our population (people who will not end up voting).
• These systematic errors in polling are called bias. We will learn more about bias in the future.

Code: Plotting margin of error in an extremely large poll over a range of values of p

library(tidyverse)
N <- 100000
p <- seq(0.35, 0.65, length = 100)
SE <- sapply(p, function(x) 2*sqrt(x*(1-x)/N))
data.frame(p = p, SE = SE) %>%
    ggplot(aes(p, SE)) +
    geom_line()