Birthdays as Instruments for Catholic School Effects

Overview

In this writeup, the goal is to walk through a third approach (in addition to adjusting for covariates and propensity score weighting) you can take towards tackling the issue of confounding. The following diagram illustrates what we’re hoping to achieve:

• We want to estimate the (average) causal effect of Catholic Schooling \(D\) on Post-Grad Earnings \(Y\).
• However, there is an unmeasured covariate, Work Ethic \(X\) (a fork relative to \(D \leftarrow X \rightarrow Y\)), which is confounding this effect \(D \rightarrow Y\).
• So, as Angrist and Krueger (1991) proposed, we can “work around” this confounding (literally, given how the diagram is set up) by using day of birth relative to school cohort as an instrument:
  • Since this day of birth has no direct causal relationship with work ethic, and
  • Has no direct causal relationship with Post-Grad Earnings (only an indirect effect on it through Catholic Schooling \(D\)),
  • This cohort-relative date of birth can play the same role that a random coin flip plays in the “gold standard” of randomized medical trials!
• Under these assumptions (which need to be argued for! And are argued for in Angrist and Krueger (1991)), we can recover the causal effect of \(D\) on \(Y\) via the “classical IV estimator”:

\[ \beta_{\text{IV}}^{D \rightarrow Y} = \frac{\text{Cov}[Z,Y]}{\text{Cov}[D,Y]} \]

Causal diagram illustrating the role of our proposed instrument \(Z\): a student’s day of birth within the cohort they enter school with. If \(Z\) has (a) no direct causal relationship with Work Ethic (\(X\)) and (b) no direct causal relationship with Post-Grad Earnings \(Y\) (that is, only an indirect effect on it through Catholic Schooling \(D\)), then we can recover the causal effect of Catholic Schooling on Post-Grad Earnings by computing \(\beta_{\text{IV}}^{D \rightarrow Y} = \frac{\text{Cov}[Z,Y]}{\text{Cov}[D,Y]}\)

In this writeup, we will generate simulated data on birthdates, school-dropout behavior, and eventual earnings, then use the above IV estimator formula to recover the causal effect of Catholic Schooling on Post-Grad Earnings. We can then verify the accuracy of this estimate, since we generated the data in the first place! The beauty of generative modeling 😉

Uniformly-Distributed Birthdays and Catholic Schooling

library(tidyverse) |> suppressPackageStartupMessages()
library(lubridate) |> suppressPackageStartupMessages()
library(extraDistr) |> suppressPackageStartupMessages()
library(scales) |> suppressPackageStartupMessages()
set.seed(5650)
n <- 10000
dob_vals <- lubridate::make_date(1999, 8, 27)
cath_vals <- rbern(n, prob=0.5)
offset_vals <- rdunif(n, min=0, max=364)
student_df <- tibble(day_of_yr=offset_vals, catholic=factor(cath_vals)) |>
  mutate(
    dob=dob_vals + duration(day_of_yr, unit="days") - duration(1, unit="days"),
    turns_18 = lubridate::as_date(dob + duration(18, unit="years"))
  ) |>
  arrange(dob)
#duration(sample(0:364, n()), unit = "days"))
student_df <- student_df |> mutate(
  year_abbr = as.character(lubridate::year(dob)),
  month_str = factor(
    paste0("",month.abb[lubridate::month(dob)],"",str_sub(year_abbr, 3, 4)),
    levels=c(
      "Aug99","Sep99","Oct99","Nov99","Dec99","Jan00","Feb00",
      "Mar00","Apr00","May00","Jun00","Jul00","Aug00"
    )
  )
)
student_df |> head()

day_of_yr	catholic	dob	turns_18	year_abbr	month_str
0	0	1999-08-26	2017-08-25	1999	Aug99
0	0	1999-08-26	2017-08-25	1999	Aug99
0	1	1999-08-26	2017-08-25	1999	Aug99
0	1	1999-08-26	2017-08-25	1999	Aug99
0	0	1999-08-26	2017-08-25	1999	Aug99
0	0	1999-08-26	2017-08-25	1999	Aug99

Number of Students Per Month

month_df <- student_df |> group_by(month_str, catholic) |> summarize(month_total=n())
month_df |> ggplot(aes(x=month_str,y=month_total, fill=catholic)) +
  geom_bar(stat='identity') +
  theme_classic()

Now the Dropping-Out Simulation

• First day of school year: Aug 26
• Last day of school: Jun 18
• \(\implies\) This cohort = all who turn 5 yrs old between Aug 27 2004 and Aug 26 2005
• \(\implies\) All born between Aug 27 1999 and Aug 26 2000

student_df <- student_df |> mutate(
  can_drop_out = as.numeric(turns_18 <= lubridate::make_date(2018, 6, 18))
)
student_df |> group_by(dob) |> summarize(sum(can_drop_out)) |> head()

dob	sum(can_drop_out)
1999-08-26	25
1999-08-27	30
1999-08-28	30
1999-08-29	29
1999-08-30	23
1999-08-31	24

What percentage of students will actually have a chance to drop out? i.e., What proportion of students are born so that they turn 18 before the last day of school?

(elig_prop <- student_df |> summarize(elig_prop=mean(can_drop_out)) |> pull())

[1] 0.8174

# National mean 8.6%
raw_dropout_rate <- 0.086
(dropout_rate <- (1/elig_prop) * raw_dropout_rate)

[1] 0.1052116

• DCPS School Year Calendar

# Simulate desire to drop out
set.seed(5651)
wtd_vals <- rbern(nrow(student_df), prob=dropout_rate)
student_df$wants_to_drop <- wtd_vals
student_df <- student_df |> mutate(
  dropout = as.numeric(wants_to_drop==1 & can_drop_out==1)
)
mean(student_df$dropout)

[1] 0.0846

student_df |> group_by(turns_18) |> summarize(
  dob_dropouts=sum(dropout),
  wants_dropout=sum(wants_to_drop),
  total=n(), 
) |> head()

turns_18	dob_dropouts	wants_dropout	total
2017-08-25	1	1	25
2017-08-26	5	5	30
2017-08-27	7	7	30
2017-08-28	3	3	29
2017-08-29	1	1	23
2017-08-30	3	3	24

Compute Cumulative Number of Dropouts Over Time

dob_df <- student_df |> group_by(dob, catholic) |> summarize(
  dob_dropout=sum(dropout)
)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by dob and catholic.
ℹ Output is grouped by dob.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(dob, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

#dob_df$drop_todate <- cumsum(dob_df$dob_dropout)
#dob_df |> group_by(catholic) |> summarize(drop_todate=cumsum(dob_dropout)) |> ungroup()
dob_df <- dob_df |> group_by(catholic) |> arrange(dob) |>
  mutate(drop_todate=cumsum(dob_dropout)) |>
  ungroup()

Plot Cumulative Dropouts by Day

dob_df |> ggplot(aes(x=dob, y=drop_todate, color=catholic)) +
  geom_line(linewidth=0.5) +
  geom_vline(
    xintercept=lubridate::make_date(2000, 6, 18),
    linetype='dashed',
  ) +
  theme_classic()

Plot Cumulative Dropouts by Month

dob_df <- dob_df |> mutate(
  year_abbr = as.character(lubridate::year(dob)),
  month_str = factor(
    paste0("",month.abb[lubridate::month(dob)],"",str_sub(year_abbr, 3, 4)),
    levels=c(
      "Aug99","Sep99","Oct99","Nov99","Dec99","Jan00","Feb00",
      "Mar00","Apr00","May00","Jun00","Jul00","Aug00"
    )
  )
)
month_drop_df <- dob_df |> group_by(month_str, catholic) |>
  summarize(month_dropout=sum(dob_dropout))

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by month_str and catholic.
ℹ Output is grouped by month_str.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(month_str, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

month_drop_df$total_dropouts <- cumsum(month_drop_df$month_dropout)
month_drop_df <- month_drop_df |> mutate(
  next_month = lead(month_str, default="Aug00"),
  prev_month_dropouts = lag(total_dropouts, default=0),
  change = total_dropouts - prev_month_dropouts
)
# month_drop_df
month_drop_df |> ggplot(aes(x=month_str, y=total_dropouts, fill=catholic)) +
  geom_bar(stat='identity') +
  # geom_segment(aes(xend=next_month), linetype="dashed") +
  geom_text(
    aes(label=month_dropout),
    position=position_dodge(width=0.9),
    vjust=-0.25, size=3.5, show.legend=FALSE
  ) +
  theme_classic()

(Sanity check: wtd_df = students who want to drop, per month)

wtd_df <- student_df |> group_by(month_str, catholic) |>
  summarize(total_wtd=sum(wants_to_drop))

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by month_str and catholic.
ℹ Output is grouped by month_str.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(month_str, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

wtd_df |> ggplot(aes(x=month_str, y=total_wtd, fill=catholic)) +
  geom_bar(stat='identity') +
  theme_classic()

Compute Total Years (with Decimals) of Schooling

grad_day <- lubridate::make_date(2018, 6, 18)
student_df <- student_df |> mutate(
  schooling_start = lubridate::make_date(2005, 8, 26),
  schooling_end = as_date(ifelse(dropout==1, as_date(turns_18), as_date(grad_day))),
  school_days = difftime(schooling_end, schooling_start, units="days"),
  school_yrs = interval(schooling_start, schooling_end) / years(1)
)
student_df |> head()

day_of_yr	catholic	dob	turns_18	year_abbr	month_str	can_drop_out	wants_to_drop	dropout	schooling_start	schooling_end	school_days	school_yrs
0	0	1999-08-26	2017-08-25	1999	Aug99	1	0	0	2005-08-26	2018-06-18	4679 days	12.81096
0	0	1999-08-26	2017-08-25	1999	Aug99	1	0	0	2005-08-26	2018-06-18	4679 days	12.81096
0	1	1999-08-26	2017-08-25	1999	Aug99	1	0	0	2005-08-26	2018-06-18	4679 days	12.81096
0	1	1999-08-26	2017-08-25	1999	Aug99	1	0	0	2005-08-26	2018-06-18	4679 days	12.81096
0	0	1999-08-26	2017-08-25	1999	Aug99	1	0	0	2005-08-26	2018-06-18	4679 days	12.81096
0	0	1999-08-26	2017-08-25	1999	Aug99	1	0	0	2005-08-26	2018-06-18	4679 days	12.81096

Simulated Earnings \(n\) Years Later

student_df$earnings_noise <- rnorm(nrow(student_df), 0, 500)
student_df <- student_df |> mutate(
  earnings = 10000 * school_yrs + 5000 * as.numeric(catholic) + earnings_noise
)
student_df |> ggplot(aes(x=earnings, fill=catholic)) + 
  geom_density() +
  theme_classic()

student_df |> ggplot(aes(x=school_yrs, fill=catholic)) +
  geom_density(alpha=0.5) +
  theme_classic()

Mean Years of Schooling by Birth Month

schooling_df <- student_df |> group_by(month_str, catholic) |> summarize(
  mean_yrs=mean(school_yrs),
  n=n()
)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by month_str and catholic.
ℹ Output is grouped by month_str.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(month_str, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# Plot mean school days by month
schooling_df |> ggplot(aes(x=month_str, y=mean_yrs, color=catholic)) +
  geom_vline(aes(xintercept=month_str), linetype='dashed', linewidth=0.5, alpha=0.18) +
  geom_point(size=2) +
  theme_classic()

(max_schooling <- max(student_df$school_yrs))

[1] 12.81096

q1 <- c("Aug99","Sep99","Oct99","Nov99")
q2 <- c("Dec99","Jan00","Feb00")
q3 <- c("Mar00","Apr00","May00")
q4 <- c("Jun00","Jul00","Aug00")
schooling_df <- schooling_df |> mutate(
  qtr = case_match(
    month_str,
    q1 ~ 1,
    q2 ~ 2,
    q3 ~ 3,
    q4 ~ 4
  )
)

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `qtr = case_match(month_str, q1 ~ 1, q2 ~ 2, q3 ~ 3, q4 ~ 4)`.
ℹ In group 1: `month_str = Aug99`.
Caused by warning:
! `case_match()` was deprecated in dplyr 1.2.0.
ℹ Please use `recode_values()` instead.

qtr_df <- schooling_df |> group_by(qtr, catholic) |> summarize(
  mean_yrs_qtr=mean(mean_yrs),
  n_qtr=sum(n)
)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by qtr and catholic.
ℹ Output is grouped by qtr.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(qtr, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

qtr_df |> ggplot(aes(x=qtr, y=mean_yrs_qtr, color=catholic)) +
  geom_point(
    stat='identity', size=2, position=position_dodge2(width=0.4)
  ) +
  geom_segment(
    aes(yend=max_schooling),
    # position=position_jitterdodge(
    #   #dodge.width=0.5
    # ),
    position=position_dodge2(width=0.4),
    linetype='dotted', linewidth=0.5
  ) +
  # ylim(12.759, 12.82) +
  geom_hline(
    aes(
      yintercept=max_schooling,
      #linetype='Max Possible'
    ),
    linewidth=0.5
  ) +
  # scale_linetype_manual("", values=c("dashed")) +
  # geom_text(x=2.5, y=12.81, label='Maximum Possible (Non-Dropout Amount)', vjust=-1) +
  labs(
    title="Schooling Gap By Quarter of Birth",
    x="DOB Quarter",
    y="Mean Years of Schooling"
  ) +
  theme_classic(base_size=14)

schooling_df <- schooling_df |> mutate(
  half = factor(
    ifelse(qtr==1 | qtr==2, 0, 1),
    levels=c(0,1)
  )
)
half_df <- schooling_df |> group_by(half, catholic) |> summarize(
  mean_yrs_half=mean(mean_yrs),
  n_half=sum(n)
)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by half and catholic.
ℹ Output is grouped by half.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(half, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

half_df |> ggplot(aes(x=half, y=mean_yrs_half, color=catholic)) +
  geom_point(
    stat='identity', size=2.5,
    position=position_dodge2(width=0.4)
  ) +
  geom_segment(
    aes(yend=max_schooling),
    position=position_dodge2(width=0.4),
    linetype='dotted', linewidth=0.6
  ) +
  # ylim(12.759, 12.82) +
  geom_hline(
    aes(
      yintercept=max_schooling,
      #linetype='Max Possible'
    ),
    linewidth=0.5
  ) +
  # scale_linetype_manual("", values=c("dashed")) +
  # geom_text(x=1.5, y=12.81, label='Maximum Possible (Non-Dropout Amount)', vjust=-1) +
  labs(
    title="Schooling Gap By DOB Half",
    x="DOB Half",
    y="Mean Years of Schooling"
  ) +
  #xlim(0.5, 2.5) +
  # ylim(12.76, 12.82) +
  theme_classic(base_size=14)

h0 <- c("Aug99","Sep99","Oct99","Nov99","Dec99","Jan00","Feb00")
h1 <- c("Mar00","Apr00","May00","Jun00","Jul00","Aug00")
student_df <- student_df |> mutate(
  half = factor(
    ifelse(month_str %in% h0, 0, 1),
    levels=c(0,1)
  )
)
earn_df <- student_df |> group_by(half, catholic) |>
  summarize(mean_earnings=mean(earnings))

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by half and catholic.
ℹ Output is grouped by half.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(half, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

earn_df

half	catholic	mean_earnings
0	0	132518.2
0	1	137524.1
1	0	132994.3
1	1	138005.8

# earn_df |> pivot_wider(-half)
earn_df |> ggplot(aes(x=half, y=log(mean_earnings), color=catholic, group=catholic)) +
  geom_point(stat='identity', size=2.5) +
  geom_line() +
  # ylim(12.759, 12.82) +
  # scale_linetype_manual("", values=c("dashed")) +
  labs(
    title="Earnings Gap By DOB Half",
    x="DOB Half",
    y="Mean Earnings"
  ) +
  #xlim(0.5, 2.5) +
  # ylim(127400, 128100) +
  theme_classic(base_size=14)

student_df <- student_df |> mutate(
  qtr = case_match(
    month_str,
    q1 ~ 1,
    q2 ~ 2,
    q3 ~ 3,
    q4 ~ 4
  )
)
earn_df <- student_df |> group_by(qtr, catholic) |>
  summarize(mean_earnings=mean(earnings))

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by qtr and catholic.
ℹ Output is grouped by qtr.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(qtr, catholic))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

earn_df

qtr	catholic	mean_earnings
1	0	132374.0
1	1	137431.5
2	0	132665.2
2	1	137625.5
3	0	132902.5
3	1	137922.4
4	0	133094.5
4	1	138098.6

# earn_df |> pivot_wider(-half)
earn_df |> ggplot(aes(x=qtr, y=log(mean_earnings), color=catholic, group=catholic)) +
  geom_point(stat='identity', size=2.5) +
  geom_line() +
  # ylim(12.759, 12.82) +
  # scale_linetype_manual("", values=c("dashed")) +
  labs(
    title="Earnings Gap By DOB Half",
    x="DOB Quarter",
    y="Mean Earnings"
  ) +
  #xlim(0.5, 2.5) +
  # ylim(127400, 128100) +
  theme_classic(base_size=14)

lm(mean_earnings ~ as.numeric(qtr), data=earn_df)

Call:
lm(formula = mean_earnings ~ as.numeric(qtr), data = earn_df)

Coefficients:
    (Intercept)  as.numeric(qtr)  
       134677.1            234.9  

lm(as.numeric(catholic) ~ as.numeric(qtr), data=earn_df)

Call:
lm(formula = as.numeric(catholic) ~ as.numeric(qtr), data = earn_df)

Coefficients:
    (Intercept)  as.numeric(qtr)  
      1.500e+00        2.402e-16  

lm(mean_earnings ~ as.numeric(catholic), data=earn_df)

Call:
lm(formula = mean_earnings ~ as.numeric(catholic), data = earn_df)

Coefficients:
         (Intercept)  as.numeric(catholic)  
              127749                  5010  

References

Angrist, Joshua D., and Alan B. Krueger. 1991. “Does Compulsory School Attendance Affect Schooling and Earnings?” The Quarterly Journal of Economics 106 (4): 979–1014. https://doi.org/10.2307/2937954.