Probability_and_Regression.Rmd

---
title: "Causal Inference: <br> *The Mixtape*"
subtitle: "<it>Probability and Regression Review</it>"
output: 
  learnr::tutorial:
    css: css/style.css
    highlight: "kate"
runtime: shiny_prerendered
---

## Welcome

This is material for the **Probability and Regression** chapter in Scott Cunningham's book, [Causal Inference: The Mixtape.](https://mixtape.scunning.com/)

### Packages needed

The first thing you need to do is install a few packages to make sure everything runs:

```{r, eval = FALSE}
install.packages("tidyverse")
install.packages("cli")
install.packages("haven")
install.packages("rmarkdown")
install.packages("learnr")
install.packages("estimatr")
install.packages("haven")

# This Chapter Only
install.packages("arm")
install.packages("mvtnorm")
install.packages("lme4")
install.packages("multiwayvcov")
install.packages("clusterSEs")
```

### Load

```{r load, warning=FALSE, message=FALSE}
library(learnr)
library(haven)
library(tidyverse)
library(estimatr)
library(stargazer)

# This Chapter Only
library(arm)
library(mvtnorm)
library(lme4)
library(multiwayvcov)
library(clusterSEs)

# 10 minute code time limit
options(tutorial.exercise.timelimit = 600)

# read_data function
read_data <- function(df) {
  full_path <- paste0("https://raw.github.com/scunning1975/mixtape/master/", df)
  return(haven::read_dta(full_path))
}
```


## OLS

```{r ols_1, exercise=TRUE, echo=FALSE}

cli::cli_alert_success("Generating Data")

set.seed(1)
tb <- tibble(
  x = rnorm(10000),
  u = rnorm(10000),
  y = 5.5*x + 12*u
) 

reg_tb <- tb %>% 
  lm(y ~ x, .) %>%
  print()

reg_tb$coefficients

tb <- tb %>% 
  mutate(
    yhat1 = predict(lm(y ~ x, .)),
    yhat2 = 0.0732608 + 5.685033*x, 
    uhat1 = residuals(lm(y ~ x, .)),
    uhat2 = y - yhat2
  )

summary(tb[-1:-3])


cli::cli_alert_success("Plotting Line of Best Fit")

ggplot(data= tb, aes(x=x, y=y)) + 
    ggtitle("OLS Regression Line") +
    geom_point(size = 0.05, color = "black", alpha = 0.5) +
    geom_smooth(aes(x=x, y=y), method = "lm", color = "black") +
    annotate("text", x = -1.5, y = 30, color = "red", 
             label = paste("Intercept = ", -0.0732608)) +
    annotate("text", x = 1.5, y = -30, color = "blue", 
             label = paste("Slope =", 5.685033))



```

#### Questions:

-   What is the predicted value of $y$ when $x = 0$?
-   How much do we estimate $y$ increases by when $x$ increases by one unit?
-   Assume we y was the natural log of some variable, and x was the natural log of some variable.  How do we interpret the coefficient on $x$ if it is a log-log regression?

```{r ols_2, exercise=TRUE, echo=FALSE}
tb <- tibble(
  x = 9*rnorm(10),
  u = 36*rnorm(10),
  y = 3 + 2*x + u,
  yhat = predict(lm(y ~ x)),
  uhat = residuals(lm(y ~ x))
)

cli::cli_alert_success("Results from OLS")
summary(tb)
colSums(tb)


```

#### Questions

-   What is the average of the residuals $\hat{u}$ from our regression?

```{r ols_3, exercise=TRUE, echo=FALSE}
cli::cli_alert_success("Running 1000 Simulations of OLS")

lm <- lapply(
  1:1000,
  function(x) tibble(
    x = 9*rnorm(10000),
    u = 36*rnorm(10000),
    y = 3 + 2*x + u
  ) %>% 
    lm(y ~ x, .)
)

as_tibble(t(sapply(lm, coef))) %>%
  summary(x)

as_tibble(t(sapply(lm, coef))) %>% 
  ggplot()+
  geom_histogram(aes(x), binwidth = 0.01)

```

#### Questions

-   Explain the concept of unbiasedness in the context of this simulation?  
-   On average, do we think the estimate is close to the true value of $\beta_1 = 2$?

```{r ols_4, exercise=TRUE, echo=FALSE}
auto <- read_data("auto.dta") %>% 
  mutate(length = length - mean(length))

lm1 <- lm(price ~ length, auto)
lm2 <- lm(price ~ length + weight + headroom + mpg, auto)


coef_lm1 <- lm1$coefficients
coef_lm2 <- lm2$coefficients
resid_lm2 <- lm2$residuals 

y_single <- tibble(price = coef_lm1[1] + coef_lm1[2]*auto$length, 
                   length = auto$length)

y_multi <- tibble(price = coef_lm1[1] + coef_lm2[2]*auto$length, 
                  length = auto$length)


ggplot(auto) + 
  geom_point(aes(x = length, y = price)) +
  geom_smooth(aes(x = length, y = price), data = y_multi, color = "blue") +
  geom_smooth(aes(x = length, y = price), data = y_single, color="red")
```

#### Questions

-   What happened to the coefficient on length after controlling for weight, headroom, and mpg in the regression?

## Clustering Standard Errors

### Cluster robust standard errors

People will try to scare you by challenging how you constructed your standard errors. Heteroskedastic errors, though, aren't the only thing you should be worried about when it comes to inference. Some phenomena do not affect observations individually, but they do affect groups of observations that involve individuals. And then they affect those individuals within the group in a common way. Say you want to estimate the effect of class size on student achievement, but you know that there exist unobservable things (like the teacher) that affect all the students equally. If we can commit to independence of these unobservables across classes, but individual student unobservables are correlated within a class, then we have a situation in which we need to cluster the standard errors. Before we dive into an example, I'd like to start with a simulation to illustrate the problem.

As a baseline for this simulation, let's begin by simulating nonclustered data and analyze least squares estimates of that nonclustered data. This will help firm up our understanding of the problems that occur with least squares when data is clustered.

First, I will create a function to generate our Monte Carlo simulation.


```{r}
#- Analysis of Clustered Data
#- Courtesy of Dr. Yuki Yanai, 
#- http://yukiyanai.github.io/teaching/rm1/contents/R/clustered-data-analysis.html

cli::cli_alert_success("Simulation Functions Created")

gen_cluster <- function(param = c(.1, .5), n = 1000, n_cluster = 50, rho = .5) {
  # Function to generate clustered data
  # Required package: mvtnorm
  
  # individual level
  Sigma_i <- matrix(c(1, 0, 0, 1 - rho), ncol = 2)
  values_i <- rmvnorm(n = n, sigma = Sigma_i)
  
  # cluster level
  cluster_name <- rep(1:n_cluster, each = n / n_cluster)
  Sigma_cl <- matrix(c(1, 0, 0, rho), ncol = 2)
  values_cl <- rmvnorm(n = n_cluster, sigma = Sigma_cl)
  
  # predictor var consists of individual- and cluster-level components
  x <- values_i[ , 1] + rep(values_cl[ , 1], each = n / n_cluster)
  
  # error consists of individual- and cluster-level components
  error <- values_i[ , 2] + rep(values_cl[ , 2], each = n / n_cluster)
  
  # data generating process
  y <- param[1] + param[2]*x + error
  
  df <- data.frame(x, y, cluster = cluster_name)
  return(df)
}

# Simulate a dataset with clusters and fit OLS
# Calculate cluster-robust SE when cluster_robust = TRUE
cluster_sim <- function(param = c(.1, .5), n = 1000, n_cluster = 50,
                        rho = .5, cluster_robust = FALSE) {
  # Required packages: mvtnorm, multiwayvcov
  df <- gen_cluster(param = param, n = n , n_cluster = n_cluster, rho = rho)
  fit <- lm(y ~ x, data = df)
  b1 <- coef(fit)[2]
  if (!cluster_robust) {
    Sigma <- vcov(fit)
    se <- sqrt(diag(Sigma)[2])
    b1_ci95 <- confint(fit)[2, ]
  } else { # cluster-robust SE
    Sigma <- cluster.vcov(fit, ~ cluster)
    se <- sqrt(diag(Sigma)[2])
    t_critical <- qt(.025, df = n - 2, lower.tail = FALSE)
    lower <- b1 - t_critical*se
    upper <- b1 + t_critical*se
    b1_ci95 <- c(lower, upper)
  }
  return(c(b1, se, b1_ci95))
}

# Function to iterate the simulation. A data frame is returned.
run_cluster_sim <- function(n_sims = 1000, param = c(.1, .5), n = 1000,
                            n_cluster = 50, rho = .5, cluster_robust = FALSE) {
  # Required packages: mvtnorm, multiwayvcov, dplyr
  df <- replicate(n_sims, cluster_sim(param = param, n = n, rho = rho,
                                      n_cluster = n_cluster,
                                      cluster_robust = cluster_robust))
  df <- as.data.frame(t(df))
  names(df) <- c('b1', 'se_b1', 'ci95_lower', 'ci95_upper')
  df <- df %>% 
    mutate(id = 1:n(),
           param_caught = ci95_lower <= param[2] & ci95_upper >= param[2])
  return(df)
}


```

### Case 1: No Correlation

For the first simulation, the data will be generated with $\beta_1 = 0$ and no correlation between the error terms.

```{r cluster_1, exercise = TRUE, eval = FALSE}

# Distribution of the estimator and confidence intervals
cli::cli_alert_success("Running Simulation")
sim_params <- c(.4, 0)   # beta1 = 0: no effect of x on y
sim_nocluster <- run_cluster_sim(n_sims = 1000, param = sim_params, rho = 0)

cli::cli_alert_success("Plot: Distribution of Betas")
(ggplot(sim_nocluster, aes(b1)) +
  geom_histogram(color = 'black') +
  geom_vline(xintercept = sim_params[2], color = 'red'))

cli::cli_alert_success("Plot: Confidence Intervals")
(ggplot(sample_n(sim_nocluster, 100),
                         aes(x = reorder(id, b1), y = b1, 
                             ymin = ci95_lower, ymax = ci95_upper,
                             color = param_caught)) +
  geom_hline(yintercept = sim_params[2], linetype = 'dashed') +
  geom_pointrange() +
  labs(x = 'sim ID', y = 'b1', title = 'Randomly Chosen 100 95% CIs') +
  scale_color_discrete(name = 'True param value', labels = c('missed', 'hit')) +
  coord_flip())



1 - sum(sim_nocluster$param_caught)/nrow(sim_nocluster)

```

#### Questions:

-   What point does the least squares estimate appear to be centered on?
-   Setting the significance level at 5%, we should incorrectly reject the null that $\beta_1=0$ about 5% of the time in our simulations. About what percent of the time does the 95% confidence intervals contain the true value of $\beta_1 = 0$?

### Case 2: Clustered Data

Now let's resimulate our data with observations that are no longer independent draws in a given cluster of observations, but the true value of $\beta_1$ still is 0.

```{r cluster_2, exercise=TRUE, echo=FALSE}

sim_params <- c(.4, 0)   # beta1 = 0: no effect of x on y
sim_cluster_ols <- run_cluster_sim(n_sims = 1000, param = sim_params)

cli::cli_alert_success("Plot: Distribution of Betas")
(ggplot(sim_cluster_ols, aes(b1)) +
  geom_histogram(color = 'black') +
  geom_vline(xintercept = sim_params[2], color = 'red'))


cli::cli_alert_success("Plot: Confidence Intervals")
(ggplot(sample_n(sim_cluster_ols, 100),
                         aes(x = reorder(id, b1), y = b1, 
                             ymin = ci95_lower, ymax = ci95_upper,
                             color = param_caught)) +
  geom_hline(yintercept = sim_params[2], linetype = 'dashed') +
  geom_pointrange() +
  labs(x = 'sim ID', y = 'b1', title = 'Randomly Chosen 100 95% CIs') +
  scale_color_discrete(name = 'True param value', labels = c('missed', 'hit')) +
  coord_flip())

1 - sum(sim_cluster_ols$param_caught)/nrow(sim_cluster_ols)
```

#### Questions:

-   When the errors are clustered, does the distribution of $\hat{\beta}_1$ estimates get wider or narrower?
-   When the errors are clustered, do we incorrectly reject the null more or less frequently?