Continuous

lm(crime_rate ~ pop_density, 
   data = communities) |> 
  coef() # extract coefficients

(Intercept) pop_density 
 -20.776945    3.093886 

ggplot(communities, 
       aes(x = pop_density, 
           y = crime_rate)) + 
  geom_point() + 
  geom_smooth(method = "lm")

Dummy variables

What is a regression model doing when you have a binary predictor?

lm(crime_rate ~ area, 
   data = communities) |>
  coef()

(Intercept)   areaUrban 
   10.43129    30.18226 

ggplot(communities, 
       aes(x = area, 
           y = crime_rate)) + 
  geom_point() + 
  geom_smooth(method = "lm")

A slight modification

lm(crime_rate ~ area, 
   data = communities) |>
  coef()

(Intercept)   areaUrban 
   10.43129    30.18226 

communities |> 
  mutate(urban = 
           as.numeric(area=="Urban")) |>
  ggplot(aes(x = urban, 
             y = crime_rate)) + 
  geom_point() + 
  geom_smooth(method = "lm")

When areaUrban is 0 (area is "Rural") the mean crime_rate is \(10.43\)

When areaUrban is 1 (area is "Urban") the mean crime_rate is \(10.43 + 30.18 = 40.61\)

Mean differences

Linear regression is just estimating means by group!

dummy_model <- 
  lm(crime_rate ~ area, 
     data = communities)
dummy_model |> coef()

(Intercept)   areaUrban 
   10.43129    30.18226 

communities |>
  group_by(area) |>
  summarize(crime_rate = 
              mean(crime_rate))

# A tibble: 2 × 2
  area  crime_rate
  <fct>      <dbl>
1 Rural       10.4
2 Urban       40.6

These are saying the exact same thing!

Hypothesis testing

Significance testing with dummy variables works like any others

summary(dummy_model)$coefficients

            Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 10.43129   1.270202  8.21231 6.660396e-15
areaUrban   30.18226   1.814574 16.63325 2.123665e-44

t.test(crime_rate ~ area, data = communities, var.equal=TRUE)

    Two Sample t-test

data:  crime_rate by area
t = -16.633, df = 298, p-value < 2.2e-16
alternative hypothesis: true difference in means between group Rural and group Urban is not equal to 0
95 percent confidence interval:
 -33.75326 -26.61125
sample estimates:
mean in group Rural mean in group Urban 
           10.43129            40.61355 

More categories

When a variable has more than two categories, all but one category is converted to a mutually exclusive binary variable

lm(crime_rate ~ incarceration, 
   data = communities) |> 
  tidy() |> 
  select(term, estimate)

# A tibble: 3 × 2
  term                estimate
  <chr>                  <dbl>
1 (Intercept)            17.7 
2 incarcerationMedium     3.22
3 incarcerationHigh      19.3 

communities |>
  group_by(incarceration) |>
  summarize(crime_rate = 
              mean(crime_rate))

# A tibble: 3 × 2
  incarceration crime_rate
  <fct>              <dbl>
1 Low                 17.7
2 Medium              20.9
3 High                37.0

Here our intercept represents the mean when incarceration is neither "Medium" nor "High". Our reference category is "Low"

ANOVA

Some of you may be familiar with ANOVA for analysing how 3+ category variables relate to a continuous outcome

aov(crime_rate ~ incarceration, data = communities) |> summary()

               Df Sum Sq Mean Sq F value  Pr(>F)    
incarceration   2  20991   10496   25.79 4.7e-11 ***
Residuals     297 120866     407                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aov(crime_rate ~ incarceration, data = communities) |> coef()

        (Intercept) incarcerationMedium   incarcerationHigh 
          17.726374            3.217159           19.265125 

As an lm()

summary(lm(crime_rate ~ incarceration, data = communities))

Call:
lm(formula = crime_rate ~ incarceration, data = communities)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.553 -14.044  -5.796  11.041  73.811 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           17.726      2.070   8.565 5.97e-16 ***
incarcerationMedium    3.217      2.850   1.129     0.26    
incarcerationHigh     19.265      2.897   6.649 1.41e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 20.17 on 297 degrees of freedom
Multiple R-squared:  0.148, Adjusted R-squared:  0.1422 
F-statistic: 25.79 on 2 and 297 DF,  p-value: 4.7e-11

Controlling for variables

Remember that we control for variables for the purpose of identifying something we’re interested in

Categorical controls

We include a categorical control \(Z\) if we believe some of \(X\) and \(Y\)’s apparent relationship is due to some \(Z\) groups having different average levels of X and Y

residualized_data <- communities |> group_by(area) |>
  mutate(pop_density_res =  pop_density - mean(pop_density), # Subtract mean!
         crime_rate_res  =  crime_rate  - mean(crime_rate))

residualized_data |> group_by(area) |>
  summarize(across(c(pop_density, crime_rate, pop_density_res, crime_rate_res),
                   mean)) 

# A tibble: 2 × 5
  area  pop_density crime_rate pop_density_res crime_rate_res
  <fct>       <dbl>      <dbl>           <dbl>          <dbl>
1 Rural        9.97       10.4       -5.63e-16       4.82e-16
2 Urban       20.0        40.6       -1.58e-15      -2.01e-15

Shifted distributions

Subtracting the mean from every value just shifts the distribution to be centred over zero.

This equalises the means of each group!

ggplot(residualized_data, 
       aes(x = pop_density,
           fill = area)) + 
  geom_histogram(alpha = 0.5, 
                 position = "identity")

ggplot(residualized_data, 
       aes(x = pop_density_res, 
           fill = area)) + 
  geom_histogram(alpha = 0.5, 
                 position = "identity")

Categorical controls

Controlling for a categorical variable is equivalent to centering the other variables by their group categories

lm(crime_rate ~ pop_density + area,
   data = communities) |> 
  tidy() |>
  select(term, estimate)

# A tibble: 3 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)   -23.0 
2 pop_density     3.36
3 areaUrban      -3.39

lm(crime_rate_res ~ pop_density_res,
   data = residualized_data)  |> 
  tidy() |> 
  select(term, estimate)

# A tibble: 2 × 2
  term            estimate
  <chr>              <dbl>
1 (Intercept)     4.03e-15
2 pop_density_res 3.36e+ 0

Note the identical pop_density slopes!

Visualising it

ggplot(communities, 
       aes(x = pop_density, 
           y = crime_rate, 
           color = area)) + 
  geom_point() + 
  geom_smooth(method = "lm", 
              color = "black") 

ggplot(residualized_data, 
       aes(x = pop_density_res,
           y = crime_rate_res,
           color = area))+ 
  geom_point() + 
  geom_smooth(method = "lm",
              color = "black")

Rural and urban data points now occupy a similar area due to subtracting their average \(X\) and \(Y\) values—this gives us their residuals

Residuals

Residuals are the difference between model predictions and the actual data

lm_res <- lm(crime_rate ~ pop_density, data = communities) |> augment()

augment() returns the original data with added residuals (.resid)

ggplot(lm_res, 
    aes(x = pop_density, 
        y = crime_rate)) + 
  geom_point() + 
  geom_smooth(method = "lm")

ggplot(lm_res,           # augment()
    aes(x = pop_density, # gives us
        y = .resid)) +   # residuals! 
  geom_point() + 
  geom_smooth(method = "lm")

car::avPlot()

This type of plot is sometimes called an added variable plot

ggplot(residualized_data, 
       aes(x = pop_density_res,
           y = crime_rate_res))+ 
  geom_point() + 
  geom_smooth(method = "lm")

lm(crime_rate ~ pop_density + 
   area, 
   data = communities) |> 
  car::avPlot(
    variable = "pop_density")

These show the regression line after controlling for other variables—and unlike ggplot2, avPlot() can accept more than one added variables at once

Continuous variables

We include continuous controls if some of \(X\) and \(Y\)’s relationship is because \(Z\) makes \(X\) and \(Y\) higher/lower—sometimes called spuriousness

residualized_data <- communities %>%
  mutate(area_num       = as.numeric(area == "Urban"),
         crime_rate_res = residuals(lm(crime_rate ~ pop_density, data = .)),
         area_res       = residuals(lm(area_num ~ pop_density, data = .)))

lm(crime_rate ~ pop_density + area, 
   data = communities) |> 
  tidy() |> select(term, estimate)

# A tibble: 3 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)   -23.0 
2 pop_density     3.36
3 areaUrban      -3.39

lm(crime_rate_res ~ area_res,
   data = residualized_data) |>
  tidy() |> select(term, estimate)

# A tibble: 2 × 2
  term         estimate
  <chr>           <dbl>
1 (Intercept)  7.77e-17
2 area_res    -3.39e+ 0

Visualising it

ggplot(residualized_data, 
       aes(x = area_num, 
           y = crime_rate, 
           color = area)) + 
  geom_point() + 
  geom_smooth(method = "lm", 
              color = "black") 

ggplot(residualized_data, 
       aes(x = area_res,
           y = crime_rate_res,
           color = area))+ 
  geom_point() + 
  geom_smooth(method = "lm",
              color = "black")

wat

Yep, that’s not a mistake—our model predicted how urban each point was and the residual is the unexpected level of urban on a 0 to 1 scale—weird!

ggplot(residualized_data, 
       aes(x = area_res,
           y = crime_rate_res,
           color = area))+ 
  geom_point() + 
  geom_smooth(method = "lm",
              color = "black")

lm(crime_rate ~ pop_density + 
       area, 
   data = communities) |> 
    car::avPlot(
        variable = "areaUrban")

Generating data

If you don’t like doing mathematical proofs to figure out if something works, you can often test things with simulation

sample_size <- 10000
sim_data <- 
  tibble(bd =
    rnorm(sample_size,
          mean = 3, sd = 1))
sim_data |>
  summarize(bd_mean = mean(bd), 
            bd_sd   = sd(bd))

# A tibble: 1 × 2
  bd_mean bd_sd
    <dbl> <dbl>
1    3.00 0.997

ggplot(sim_data, 
       aes(x = bd)) +
  geom_density()

Simulation

Now we can add other variables to our data to generate a dataset that obeys the DAG on the right

sim_data <- sim_data |>
  mutate(
    x = rnorm(n(), 0.3*bd, 1),
    y = rnorm(n(), 3 + 0.3*x + 0.3*bd, 1),
    coll = rnorm(n(), -0.2*x + -0.2*y, 1))

Estimates with back doors

Because we know the data generating process—we created it—we know what our estimates should be

lm(y ~ x + bd,
   data = sim_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 3 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    3.00 
2 x              0.301
3 bd             0.299

lm(y ~ x , # Leaving out bd 
   data = sim_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 2 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    3.82 
2 x              0.384

Excluding bd makes the effect of x larger than it should be—it is magnified

Estimates with colliders

Again, we know what effects should be identified and how to identify them. What if we add the collider?

lm(y ~ x + bd,
   data = sim_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 3 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    3.00 
2 x              0.301
3 bd             0.299

lm(y ~ x + bd + coll,
   data = sim_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 4 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    2.88 
2 x              0.249
3 bd             0.285
4 coll          -0.204

Selection (collider) bias

Let’s create the data from CRM yesterday.

coll_data <- tibble(
  Skill     = runif(500, 0, 10), 
  Frequency = runif(500, 0, 10), 
  Arrested  = Frequency > Skill) 
arrested <- coll_data |> filter(Arrested) |> mutate(Who = "Arrested")
everyone <- coll_data |> mutate(Who = "Everyone")

What’s going on here?

• No relationship between Skill at crime and Frequency of crime
  • Just two uniformly distributed variables
• Anyone who commits more crime than their skill gets arrested
  • Only these people are in arrested

Collider model

lm(Frequency ~ Skill,
   data = everyone) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 2 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)  4.86   
2 Skill       -0.00325

lm(Frequency ~ Skill,
   data = arrested) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 2 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    4.44 
2 Skill          0.635

What about controlling?

lm(Frequency ~ Skill + Arrested,
   data = everyone) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 3 × 2
  term         estimate
  <chr>           <dbl>
1 (Intercept)    -0.513
2 Skill           0.548
3 ArrestedTRUE    5.26 

lm(Frequency ~ Skill,
   data = arrested) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 2 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    4.44 
2 Skill          0.635

Visualizing collider bias

bind_rows(arrested, everyone) |>
  ggplot(aes(x = Skill, y = Frequency, color = Who)) + 
    facet_wrap(~Who) + 
    geom_point() +
    geom_smooth(method = "lm", color = "black") +
    labs(x = "Skill at Crime", y = "Frequency of Crime") +
    theme_minimal(base_size = 16) + theme(legend.position = "none")

An advanced example

You don’t need to know this…

… but it shows what’s really going on

Non-randomized treatment

Let’s create example data with actual potential outcomes

n <- 10000
po_data <- tibble(
  bd = runif( n, 0,        1), # Random uniform variable
  x  = rbinom(n, 1,        bd), # Random binary variable
  y0 = rnorm( n, 2*bd,     1),  # untreated outcome
  y1 = rnorm( n, 2*bd + 1, 1)) |> # treated outcome
  mutate(y = ifelse(x==1, y1, y0)) # Treatment just selects outcome

• x is a binary treatment with effect size of 1
• y0 is the outcome if the unit is untreated
• y1 is the outcome if the unit is treated
• bd is a backdoor predicting treatment and y
  • bd improves both potential outcomes the same amount
  • bd also makes treatment more likely

Estimation

Like usual, if we can see the back door, we can identify the effect of \(X\) on \(Y\)

lm(y ~ x + bd, data = po_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 3 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)   0.0209
2 x             0.987 
3 bd            1.95  

lm(y ~ x, data = po_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 2 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    0.663
2 x              1.63 

sample(0:1, 20, replace = TRUE) # randomly select 0 or 1 20 times

 [1] 0 0 0 0 0 1 1 1 0 1 0 0 1 1 0 0 0 1 0 1

Random assignment

Let’s select random units to receive the treatment

po_data <- po_data |>
  mutate(treat = sample(0:1, n(), replace=TRUE),
         yt = ifelse(treat==1, y1, y0))

All the treatment does is determine which potential outcome we see

lm(y ~ x + bd, data = po_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 3 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)   0.0209
2 x             0.987 
3 bd            1.95  

lm(yt ~ treat, data = po_data) |>
  tidy() |> 
  select(term, estimate)

# A tibble: 2 × 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)    1.01 
2 treat          0.950

Randomization DAGs

Randomization breaks back doors—even when we can’t observe the responsible variables!