# Libraries
pacman::p_load(haven, estimatr, texreg, janitor, tidyverse, skimr)

Models

M1

lm1 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_pct ~ female * fac_female, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_z ~ female * fac_female, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac", "female" = "female_stu", "fac_female" = "female_fac"))
  Course grade (pct) Course grade (z)
female_stu x female_fac -0.400 -0.018
  (1.225) (0.043)
female_stu 12.487*** 0.436***
  (0.796) (0.027)
female_fac 1.229 0.044
  (0.999) (0.034)
Num. obs. 28878 28878
***p < 0.01; **p < 0.05; *p < 0.1

M2

lm1 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_pct ~ female * fac_female + reservation_stu + age + b_i_jeemain_score + b_math_g1_score + b_physics_g1_score + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac + miss_reservation_stu + miss_age + miss_b_i_jeemain_score + miss_b_math_g1_score + miss_b_physics_g1_score +  miss_fac_associate_professor + miss_fac_professor + miss_fac_yearsinhighed + miss_fac_highest_degree_phd + miss_fac_highest_degree_phd_in_prog + miss_fac_degree_college_elite + miss_reservation_fac + b_ses_1 + b_ses_2 + b_ses_3 + b_ses_4 + b_ses_5 + b_ses_6 + b_i_ses_1 + b_i_ses_2 + b_i_ses_3 + b_ct_score + b_ql_score + miss_b_ct_score + miss_b_ql_score + b_i_enter_engineering_exam_2 + b_i_enter_engineering_exam_3 + b_i_enter_engineering_exam_NA + area_2 + area_3 + area_NA + hstype_2 + hstype_3 + hstype_4 + hstype_5 + hstype_6 + hstype_7 + hstype_NA + father_ed_ind_2 + father_ed_ind_3 + father_ed_ind_4 + father_ed_ind_5 + father_ed_ind_6 + father_ed_ind_7 + father_ed_ind_8 + father_ed_ind_NA + mother_ed_ind_2 + mother_ed_ind_3 + mother_ed_ind_4 + mother_ed_ind_5 + mother_ed_ind_6 + mother_ed_ind_7 + mother_ed_ind_8 + mother_ed_ind_NA + old_sib_factor_0 + old_sib_factor_2 + old_sib_factor_3 + old_sib_factor_4 + old_sib_factor_more + old_sib_factor_NA + young_sib_factor_0 + young_sib_factor_2 + young_sib_factor_3 + young_sib_factor_4 + young_sib_factor_more + young_sib_factor_NA + `school_years_english_factor_1-9` + `school_years_english_factor_10-12` + `school_years_english_factor_13-16` + school_years_english_factor_NA + district_clean, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_z ~ female * fac_female + reservation_stu + age + b_i_jeemain_score + b_math_g1_score + b_physics_g1_score + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac + miss_reservation_stu + miss_age + miss_b_i_jeemain_score + miss_b_math_g1_score + miss_b_physics_g1_score +  miss_fac_associate_professor + miss_fac_professor + miss_fac_yearsinhighed + miss_fac_highest_degree_phd + miss_fac_highest_degree_phd_in_prog + miss_fac_degree_college_elite + miss_reservation_fac + b_ses_1 + b_ses_2 + b_ses_3 + b_ses_4 + b_ses_5 + b_ses_6 + b_i_ses_1 + b_i_ses_2 + b_i_ses_3 + b_ct_score + b_ql_score + miss_b_ct_score + miss_b_ql_score + b_i_enter_engineering_exam_2 + b_i_enter_engineering_exam_3 + b_i_enter_engineering_exam_NA + area_2 + area_3 + area_NA + hstype_2 + hstype_3 + hstype_4 + hstype_5 + hstype_6 + hstype_7 + hstype_NA + father_ed_ind_2 + father_ed_ind_3 + father_ed_ind_4 + father_ed_ind_5 + father_ed_ind_6 + father_ed_ind_7 + father_ed_ind_8 + father_ed_ind_NA + mother_ed_ind_2 + mother_ed_ind_3 + mother_ed_ind_4 + mother_ed_ind_5 + mother_ed_ind_6 + mother_ed_ind_7 + mother_ed_ind_8 + mother_ed_ind_NA + old_sib_factor_0 + old_sib_factor_2 + old_sib_factor_3 + old_sib_factor_4 + old_sib_factor_more + old_sib_factor_NA + young_sib_factor_0 + young_sib_factor_2 + young_sib_factor_3 + young_sib_factor_4 + young_sib_factor_more + young_sib_factor_NA + `school_years_english_factor_1-9` + `school_years_english_factor_10-12` + `school_years_english_factor_13-16` + school_years_english_factor_NA + district_clean, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac", "female" = "female_stu", "fac_female" = "female_fac"))
  Course grade (pct) Course grade (z)
female_stu x female_fac 0.528 0.013
  (1.137) (0.039)
female_stu 12.747*** 0.447***
  (0.805) (0.028)
female_fac 1.195 0.042
  (0.826) (0.028)
Num. obs. 28469 28469
***p < 0.01; **p < 0.05; *p < 0.1

M3

lm1 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_pct ~ female * fac_female + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_z ~ female * fac_female + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac", "female" = "female_stu", "fac_female" = "female_fac"))
  Course grade (pct) Course grade (z)
female_stu x female_fac -0.424 -0.019
  (1.232) (0.043)
female_stu 12.512*** 0.436***
  (0.797) (0.028)
female_fac 1.279 0.044
  (0.908) (0.031)
Num. obs. 28878 28878
***p < 0.01; **p < 0.05; *p < 0.1

M4

lm1 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_pct ~ female * fac_female + reservation_stu + age + b_i_jeemain_score + b_math_g1_score + b_physics_g1_score + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac + miss_reservation_stu + miss_age + miss_b_i_jeemain_score + miss_b_math_g1_score + miss_b_physics_g1_score +  miss_fac_associate_professor + miss_fac_professor + miss_fac_yearsinhighed + miss_fac_highest_degree_phd + miss_fac_highest_degree_phd_in_prog + miss_fac_degree_college_elite + miss_reservation_fac + b_ses_1 + b_ses_2 + b_ses_3 + b_ses_4 + b_ses_5 + b_ses_6 + b_i_ses_1 + b_i_ses_2 + b_i_ses_3 + b_ct_score + b_ql_score + miss_b_ct_score + miss_b_ql_score + b_i_enter_engineering_exam_2 + b_i_enter_engineering_exam_3 + b_i_enter_engineering_exam_NA + area_2 + area_3 + area_NA + hstype_2 + hstype_3 + hstype_4 + hstype_5 + hstype_6 + hstype_7 + hstype_NA + father_ed_ind_2 + father_ed_ind_3 + father_ed_ind_4 + father_ed_ind_5 + father_ed_ind_6 + father_ed_ind_7 + father_ed_ind_8 + father_ed_ind_NA + mother_ed_ind_2 + mother_ed_ind_3 + mother_ed_ind_4 + mother_ed_ind_5 + mother_ed_ind_6 + mother_ed_ind_7 + mother_ed_ind_8 + mother_ed_ind_NA + old_sib_factor_0 + old_sib_factor_2 + old_sib_factor_3 + old_sib_factor_4 + old_sib_factor_more + old_sib_factor_NA + young_sib_factor_0 + young_sib_factor_2 + young_sib_factor_3 + young_sib_factor_4 + young_sib_factor_more + young_sib_factor_NA + `school_years_english_factor_1-9` + `school_years_english_factor_10-12` + `school_years_english_factor_13-16` + school_years_english_factor_NA + district_clean + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_z ~ female * fac_female + reservation_stu + age + b_i_jeemain_score + b_math_g1_score + b_physics_g1_score + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac + miss_reservation_stu + miss_age + miss_b_i_jeemain_score + miss_b_math_g1_score + miss_b_physics_g1_score +  miss_fac_associate_professor + miss_fac_professor + miss_fac_yearsinhighed + miss_fac_highest_degree_phd + miss_fac_highest_degree_phd_in_prog + miss_fac_degree_college_elite + miss_reservation_fac + b_ses_1 + b_ses_2 + b_ses_3 + b_ses_4 + b_ses_5 + b_ses_6 + b_i_ses_1 + b_i_ses_2 + b_i_ses_3 + b_ct_score + b_ql_score + miss_b_ct_score + miss_b_ql_score + b_i_enter_engineering_exam_2 + b_i_enter_engineering_exam_3 + b_i_enter_engineering_exam_NA + area_2 + area_3 + area_NA + hstype_2 + hstype_3 + hstype_4 + hstype_5 + hstype_6 + hstype_7 + hstype_NA + father_ed_ind_2 + father_ed_ind_3 + father_ed_ind_4 + father_ed_ind_5 + father_ed_ind_6 + father_ed_ind_7 + father_ed_ind_8 + father_ed_ind_NA + mother_ed_ind_2 + mother_ed_ind_3 + mother_ed_ind_4 + mother_ed_ind_5 + mother_ed_ind_6 + mother_ed_ind_7 + mother_ed_ind_8 + mother_ed_ind_NA + old_sib_factor_0 + old_sib_factor_2 + old_sib_factor_3 + old_sib_factor_4 + old_sib_factor_more + old_sib_factor_NA + young_sib_factor_0 + young_sib_factor_2 + young_sib_factor_3 + young_sib_factor_4 + young_sib_factor_more + young_sib_factor_NA + `school_years_english_factor_1-9` + `school_years_english_factor_10-12` + `school_years_english_factor_13-16` + school_years_english_factor_NA + district_clean + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac, data = ., fixed_effects = ~ course_fe, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac", "female" = "female_stu", "fac_female" = "female_fac"))
  Course grade (pct) Course grade (z)
female_stu x female_fac 0.528 0.013
  (1.137) (0.039)
female_stu 12.747*** 0.447***
  (0.805) (0.028)
female_fac 1.195 0.042
  (0.826) (0.028)
Num. obs. 28469 28469
***p < 0.01; **p < 0.05; *p < 0.1

M5

lm1 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_pct ~ female * fac_female + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac, data = ., fixed_effects = ~ course_fe + stdid, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_z ~ female * fac_female + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac, data = ., fixed_effects = ~ course_fe + stdid, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac", "fac_female" = "female_fac"))
  Course grade (pct) Course grade (z)
female_stu x female_fac 2.227** 0.075**
  (0.890) (0.030)
female_fac 0.088 0.006
  (0.747) (0.025)
Num. obs. 28878 28878
***p < 0.01; **p < 0.05; *p < 0.1

M6

lm1 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_pct ~ female * fac_female + reservation_stu + age + b_i_jeemain_score + b_math_g1_score + b_physics_g1_score + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac + miss_reservation_stu + miss_age + miss_b_i_jeemain_score + miss_b_math_g1_score + miss_b_physics_g1_score +  miss_fac_associate_professor + miss_fac_professor + miss_fac_yearsinhighed + miss_fac_highest_degree_phd + miss_fac_highest_degree_phd_in_prog + miss_fac_degree_college_elite + miss_reservation_fac + b_ses_1 + b_ses_2 + b_ses_3 + b_ses_4 + b_ses_5 + b_ses_6 + b_i_ses_1 + b_i_ses_2 + b_i_ses_3 + b_ct_score + b_ql_score + miss_b_ct_score + miss_b_ql_score + b_i_enter_engineering_exam_2 + b_i_enter_engineering_exam_3 + b_i_enter_engineering_exam_NA + area_2 + area_3 + area_NA + hstype_2 + hstype_3 + hstype_4 + hstype_5 + hstype_6 + hstype_7 + hstype_NA + father_ed_ind_2 + father_ed_ind_3 + father_ed_ind_4 + father_ed_ind_5 + father_ed_ind_6 + father_ed_ind_7 + father_ed_ind_8 + father_ed_ind_NA + mother_ed_ind_2 + mother_ed_ind_3 + mother_ed_ind_4 + mother_ed_ind_5 + mother_ed_ind_6 + mother_ed_ind_7 + mother_ed_ind_8 + mother_ed_ind_NA + old_sib_factor_0 + old_sib_factor_2 + old_sib_factor_3 + old_sib_factor_4 + old_sib_factor_more + old_sib_factor_NA + young_sib_factor_0 + young_sib_factor_2 + young_sib_factor_3 + young_sib_factor_4 + young_sib_factor_more + young_sib_factor_NA + `school_years_english_factor_1-9` + `school_years_english_factor_10-12` + `school_years_english_factor_13-16` + school_years_english_factor_NA + district_clean, data = ., fixed_effects = ~ course_fe + facid, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_z ~ female * fac_female + reservation_stu + age + b_i_jeemain_score + b_math_g1_score + b_physics_g1_score + fac_associate_professor + fac_professor + fac_yearsinhighed + fac_highest_degree_phd + fac_highest_degree_phd_in_prog + fac_degree_college_elite + reservation_fac + miss_reservation_stu + miss_age + miss_b_i_jeemain_score + miss_b_math_g1_score + miss_b_physics_g1_score +  miss_fac_associate_professor + miss_fac_professor + miss_fac_yearsinhighed + miss_fac_highest_degree_phd + miss_fac_highest_degree_phd_in_prog + miss_fac_degree_college_elite + miss_reservation_fac + b_ses_1 + b_ses_2 + b_ses_3 + b_ses_4 + b_ses_5 + b_ses_6 + b_i_ses_1 + b_i_ses_2 + b_i_ses_3 + b_ct_score + b_ql_score + miss_b_ct_score + miss_b_ql_score + b_i_enter_engineering_exam_2 + b_i_enter_engineering_exam_3 + b_i_enter_engineering_exam_NA + area_2 + area_3 + area_NA + hstype_2 + hstype_3 + hstype_4 + hstype_5 + hstype_6 + hstype_7 + hstype_NA + father_ed_ind_2 + father_ed_ind_3 + father_ed_ind_4 + father_ed_ind_5 + father_ed_ind_6 + father_ed_ind_7 + father_ed_ind_8 + father_ed_ind_NA + mother_ed_ind_2 + mother_ed_ind_3 + mother_ed_ind_4 + mother_ed_ind_5 + mother_ed_ind_6 + mother_ed_ind_7 + mother_ed_ind_8 + mother_ed_ind_NA + old_sib_factor_0 + old_sib_factor_2 + old_sib_factor_3 + old_sib_factor_4 + old_sib_factor_more + old_sib_factor_NA + young_sib_factor_0 + young_sib_factor_2 + young_sib_factor_3 + young_sib_factor_4 + young_sib_factor_more + young_sib_factor_NA + `school_years_english_factor_1-9` + `school_years_english_factor_10-12` + `school_years_english_factor_13-16` + school_years_english_factor_NA + district_clean, data = ., fixed_effects = ~ course_fe + facid, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac", "female" = "female_stu"))
  Course grade (pct) Course grade (z)
female_stu x female_fac 0.925 0.026
  (1.171) (0.040)
female_stu 12.618*** 0.442***
  (0.820) (0.028)
Num. obs. 28469 28469
***p < 0.01; **p < 0.05; *p < 0.1

M7

lm1 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_pct ~ female * fac_female, data = ., fixed_effects = ~ course_fe + stdid + facid, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(course_fe = str_c(department_id, course_name, sep = "_")) %>%
  lm_robust(course_rank_z ~ female * fac_female, data = ., fixed_effects = ~ course_fe + stdid + facid, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac"))
  Course grade (pct) Course grade (z)
female_stu x female_fac 2.601*** 0.087***
  (0.926) (0.031)
Num. obs. 28878 28878
***p < 0.01; **p < 0.05; *p < 0.1

M8

lm1 <-
  df_grades_sfix %>%
  mutate(
    course_fe = str_c(department_id, course_name, sep = "_"),
    classroom_fe = str_c(course_fe, facid, sep = "_")
  ) %>%
  lm_robust(course_rank_pct ~ female * fac_female, data = ., fixed_effects = ~ course_fe + stdid + classroom_fe, se_type = "stata", clusters = facid)

lm2 <-
  df_grades_sfix %>%
  mutate(
    course_fe = str_c(department_id, course_name, sep = "_"),
    classroom_fe = str_c(course_fe, facid, sep = "_")
  ) %>%
  lm_robust(course_rank_z ~ female * fac_female, data = ., fixed_effects = ~ course_fe + stdid + classroom_fe, se_type = "stata", clusters = facid)

knitreg(list(lm1, lm2), custom.note = "%stars", stars = c(0.01, 0.05, 0.1), include.ci = F, dcolumn = T, booktabs = T, float.pos = "H", caption = "", digits = 3, include.nclust = F, include.rsquared = F, include.adjrs = F, include.rmse = F, custom.model.names = c("Course grade (pct)", "Course grade (z)"), custom.coef.map = list("female:fac_female" = "female_stu x female_fac"))
  Course grade (pct) Course grade (z)
female_stu x female_fac 2.638*** 0.087***
  (0.953) (0.032)
Num. obs. 28878 28878
***p < 0.01; **p < 0.05; *p < 0.1