---
title: 'Technical Appendix C: Retention Analysis'
output:
  pdf_document: default
  html_document: default
---

# Part I: Exploratory Data Analysis & T-tests

## Data Wrangling 

```{r, message=FALSE, warning=FALSE}
# Import packages 
library(plyr)
library(dplyr)
library(lubridate)
library(tidyverse)
library(readxl)
library(ggplot2)
library(arm)
```

```{r}
# Import scholarship data 
scholarship <- read_excel("DR-826 Pittsburgh Promise scholarship use.xlsx", sheet = "Sheet1")

# Clean scholarship data
scholarship <- scholarship %>%
  mutate(QualifiedforCorePromise = ifelse(QualifiedforCorePromise == "yes",1,0),
         QualifiedforExtensionPromise = ifelse(QualifiedforExtensionPromise == "yes",1,0),
         EverReceivedPromiseAward = ifelse(EverReceivedPromiseAward == "yes",1,0),
         StillReceivingAward = ifelse(StillReceivingAward == "yes",1,0),
         StillEligible = ifelse(StillEligible == "yes",1,0),
         HighSchool = as.factor(HighSchool)) %>%
  rename(RandomID = "Random ID")
head(scholarship)
```

```{r}
# Import NSC data
nsc <- read_excel("DR-844 NSC data_not including 1920.xlsx", sheet = "NSC data")
nsc$RandomID <- as.character(nsc$RandomID)
head(nsc)
```

```{r}
# Import demographics data
demo1415 <- read_excel("DR-827 Demographics.xlsx", sheet = "Demographics for cohort 1415")
demo1516 <- read_excel("DR-827 Demographics.xlsx", sheet = "Demographics for cohort 1516")
demo1617 <- read_excel("DR-827 Demographics.xlsx", sheet = "Demographics for cohort 1617")

# Clean demographics data 
demographics <- rbind(demo1415,demo1516,demo1617)
demographics <- demographics %>% 
  dplyr::select(RandomID, Race) %>% 
  distinct()

# Delete students who have more than 1 races
demographics <- demographics %>% 
  filter(! RandomID %in% c(6008133,6040930,
                           6126484,6137723,
                           6207255,6213261,
                           6228923,6510341))
```


## Comparison I: Retention vs. Scholarship Receipt

### EDA

```{r, message=FALSE, warning=FALSE}
# Prepare data sets for the visualization of comparison plot
nsc_processed = nsc %>%
  filter(ENROLLMENT_BEGIN >0) %>%
  filter(ENROLLMENT_BEGIN %/% 10000 != 2016) %>%
  filter(COLLEGE_STATE == "PA") %>%
  mutate(ENROLLMENT_BEGIN = as.Date(paste(substr(ENROLLMENT_BEGIN, 1, 4),
                                          substr(ENROLLMENT_BEGIN, 5, 6),
                                          substr(ENROLLMENT_BEGIN, 7, 8),
                                          sep = "-")),
         ENROLLMENT_END = as.Date(paste(substr(ENROLLMENT_END, 1, 4),
                                        substr(ENROLLMENT_END, 5, 6),
                                        substr(ENROLLMENT_END, 7, 8),
                                        sep = "-")),
         elapsed = ENROLLMENT_END - ENROLLMENT_BEGIN)
```

```{r, message=FALSE, warning=FALSE}
# Issue: some students enroll in fall while others enroll in spring
# Create a new variable semester to indicate which semester students enroll in
nsc_retention = nsc_processed %>% 
  group_by(RandomID) %>% 
  summarise(retention = sum(elapsed),
            ENROLLMENT_BEGIN_year = year(min(ENROLLMENT_BEGIN)),
            semester = ifelse(month(min(ENROLLMENT_BEGIN))>6, "Fall","Spring"))
```

```{r}
# Prepare data sets for the visualization of comparison plot
nsc_promise = nsc_retention %>% 
  left_join(scholarship, by = "RandomID") %>%
  mutate(EverReceivedPromiseAward =
           ifelse(EverReceivedPromiseAward==1,"yes","no"),
         EverReceivedPromiseAward = replace_na(EverReceivedPromiseAward,"no"),
         EverReceivedPromiseAward = as.factor(EverReceivedPromiseAward),
         retention = as.numeric(retention)) %>%
  dplyr::select(-c(QualifiedforCorePromise,
            QualifiedforExtensionPromise,
            StillReceivingAward,
            StillEligible,HighSchool))
```

```{r}
# Visualization of retention comparison under receipt of scholarship
ggplot(nsc_promise, aes(x=as.factor(ENROLLMENT_BEGIN_year),
                        y=retention,
                        fill=EverReceivedPromiseAward)) +
  geom_boxplot() +
  labs(x = "Enrollment Begin",
       y = "Retention(days)",
       title = "Retention Comparison under Receipt of Scholarship") +
  theme_bw()
```

**Findings:**

* We see that students who received scholarship tend to have better retention except for year 2019. 
* However, notice that we only have 13 observations for year 2017, and 93 observations for year 2020.

```{r}
# Take a closer look by semester they enroll
nsc_promise %>% 
  group_by(ENROLLMENT_BEGIN_year, semester) %>% 
  count()

# Visualization of retention comparison under receipt of scholarship (by semester)
# Since we have limited observations for 2017 and all spring semesters, 
# the resulting plot is just for reference
ggplot(nsc_promise, aes(x=as.factor(ENROLLMENT_BEGIN_year),
                        y=retention,
                        fill=EverReceivedPromiseAward)) +
  geom_boxplot() +
  labs(x = "Enrollment Begin",
       y = "Retention(days)") +
  theme_bw() +
  facet_wrap(~semester)
```


### Check for Significance for 2018 & 2019 (T-tests)

```{r}
# Prepare data sets for 2018 & 2019
nsc_promise_2018 = nsc_promise %>% filter(ENROLLMENT_BEGIN_year == 2018)
nsc_promise_2019 = nsc_promise %>% filter(ENROLLMENT_BEGIN_year == 2019)

# Use Bartlett test to check for equal variance or not (2018)
bartlett.test(retention~EverReceivedPromiseAward, data = nsc_promise_2018) #unequal variance 
# T-test to check for significance of difference in retention (2018)
oneway.test(retention~EverReceivedPromiseAward, data = nsc_promise_2018, var.equal = FALSE)

# Use Bartlett test to check for equal variance or not (2019)
bartlett.test(retention~EverReceivedPromiseAward, data = nsc_promise_2019) #equal variance
# T-test to check for significance of difference in retention (2019)
oneway.test(retention~EverReceivedPromiseAward, data = nsc_promise_2019, var.equal = TRUE)
```

**Findings:**

* For 2018, the difference is significant; for 2019, the difference is not significant.
* Findings from significance tests align with what we observe in the box plot above.
* Thus, the difference in retention is quite obvious in 2018, whereas in 2019 is not very obvious. Maybe the retention difference would become more obvious for more senior
students.


## Comparison II: Retention vs.Race

### EDA

```{r}
# Check racial composition
demographics %>% group_by(Race) %>% count() %>% arrange(-n)
```

* We see that the majority of students are blacks and whites, so we will explore the difference in student retention between these two races and categorize the rest racial groups as "Other".
* Also, we restrict on students who went to college in PA.

```{r}
# Prepare data sets for the visualization of the comparison plot
demographics$RandomID <- as.character(demographics$RandomID)
nsc_demographics = nsc_promise %>%
  left_join(demographics, by = "RandomID") %>%
  mutate(Race = ifelse(Race != "White" & Race != "African American", 
                       "Others", 
                       ifelse(Race == "African American", 
                              "Black", 
                              "White"))) %>%
  na.omit(retention)
```

```{r}
# Visualization of retention comparison under race 
ggplot(nsc_demographics, aes(x=as.factor(ENROLLMENT_BEGIN_year),
                             y=retention, 
                             fill= Race)) +
  geom_boxplot() +
  labs(x = "Enrollment Begin",
       y = "Retention(days)",
       title = "Retention Comparison by Race") +
  theme_bw()
```


### Check for Significance for 2018 & 2019 (T-tests)

```{r}
# Prepare data set for the visualization of comparison plot
nsc_demographics2 = nsc_promise %>%
  left_join(demographics, by = "RandomID") %>%
  filter(Race %in% c("White","African American")) %>%
  mutate(Race = ifelse(Race == "White", "White","Black"))

# Prepare data sets for 2018 & 2019
nsc_demographics_2018 = nsc_demographics2 %>%
  filter(ENROLLMENT_BEGIN_year == 2018)
nsc_demographics_2019 = nsc_demographics2 %>%
  filter(ENROLLMENT_BEGIN_year == 2019)

# Use Bartlett test to check for equal variance or not (2018)
bartlett.test(retention~Race, data = nsc_demographics_2018) #unequal variance 
# T-test to check for significance of difference in retention (2018)
oneway.test(retention~Race, data = nsc_demographics_2018, var.equal = FALSE)

# Use Bartlett test to check for equal variance or not (2019)
bartlett.test(retention~Race, data = nsc_demographics_2019) #equal variance
# T-test to check for significance of difference in retention (2019)
oneway.test(retention~Race, data = nsc_demographics_2019, var.equal = TRUE)
```

**Findings:**

* The differences are significant for both 2018 and 2019.
* Findings from significance tests align with our what we observe in the box plot above.
  + Well aligns with what we observe for 2018 enrollment year
  + Not perfectly aligns with what we see for 2019 enrollment year (i.e. no significant difference in retention between racial groups)


## Comparison III: Retention vs. Interaction between Scholarship receipt and Race

### EDA 

```{r}
# Students who received the scholarship
# Visualization of retention comparison under race
nsc_demographics2 %>%
  filter(EverReceivedPromiseAward=="yes") %>%
  ggplot(aes(x=as.factor(ENROLLMENT_BEGIN_year), 
             y=retention, fill=Race)) +
  geom_boxplot() +
  labs(x = "Enrollment Begin",
       y = "Retention(days)",
       title = "Retention for Students Who Received Scholarship") +
  theme_bw()
```

```{r}
# Filter students who did not receive the scholarship
# Visualization of retention comparison under race
nsc_demographics2 %>%
  filter(EverReceivedPromiseAward=="no") %>%
  ggplot(aes(x=as.factor(ENROLLMENT_BEGIN_year), 
             y=retention, 
             fill=Race)) +
  geom_boxplot() +
  labs(x = "Enrollment Begin",
       y = "Retention(days)",
       title = "Retention for Students Who Not Received Scholarship") +
  theme_bw()
```


### Check for Significance for 2018 & 2019 

```{r}
# Prepare data sets for 2018 & 2019
nsc_demographics2_2018 = nsc_demographics2 %>%
  filter(ENROLLMENT_BEGIN_year == 2018,
         EverReceivedPromiseAward == "yes")
nsc_demographics2_2019 = nsc_demographics2 %>%
  filter(ENROLLMENT_BEGIN_year == 2019,
         EverReceivedPromiseAward == "yes")

# Use Bartlett test to check for equal variance or not (2018)
bartlett.test(retention~Race, data = nsc_demographics2_2018) #equal variance 
# T-test to check for significance of difference in retention (2018)
oneway.test(retention~Race, data = nsc_demographics2_2018, var.equal = TRUE)

# Use Bartlett test to check for equal variance or not (2019)
bartlett.test(retention~Race, data = nsc_demographics2_2019) #equal variance
# T-test to check for significance of difference in retention (2019)
oneway.test(retention~Race, data = nsc_demographics2_2019, var.equal = TRUE)
```

```{r}
# Prepare data sets for 2018 & 2019
nsc_demographics3_2018 = nsc_demographics2 %>%
  filter(ENROLLMENT_BEGIN_year == 2018,
         EverReceivedPromiseAward == "no")
nsc_demographics3_2019 = nsc_demographics2 %>%
  filter(ENROLLMENT_BEGIN_year == 2019,
         EverReceivedPromiseAward == "no")

# Use Bartlett test to check for equal variance or not (2018)
bartlett.test(retention~Race, data = nsc_demographics3_2018) #equal variance 
# T-test to check for significance of difference in retention (2018)
oneway.test(retention~Race, data = nsc_demographics3_2018, var.equal = TRUE)

# Use Bartlett test to check for equal variance or not (2019)
bartlett.test(retention~Race, data = nsc_demographics3_2019) #equal variance
# T-test to check for significance of difference in retention (2019)
oneway.test(retention~Race, data = nsc_demographics3_2019, var.equal = TRUE)
```

**Findings:**

* Among students who received the scholarship, we see that for 2018, there is a significant racial difference, but for 2019 the difference is not significant.
* Among students who did not receive the scholarship, we see reverse happens, the difference is not significant for 2018 but significant for 2019.

\newpage

# Part II: Multivariate Regression Analysis for Whole Range

## Data Wrangling

```{r, message=FALSE, warning=FALSE}
# Import data sets
# data_qualified is the final data set used in logistic regression modeling from scholarship analysis
# retention_demographics is the joined data set between NSC and demographics 
data_qualified <- read_csv("data_qualified.csv")
retention_demographics <- read_csv("retention_demographics.csv") %>% dplyr::select(-c(1))
retention_demographics$EverReceivedPromiseAward <- 
  ifelse(retention_demographics$EverReceivedPromiseAward == "yes", 1, 0)
```

```{r, message=FALSE, warning=FALSE}
# Inner join data_qualified and retention_demographics to obtain data set used for regression
retention_reg <- inner_join(data_qualified, retention_demographics, by = c("RandomID", 
                                                                           "Race", 
                                                                           "GradYear", 
                                                                           "EverReceivedPromiseAward"))

# Categorize race as white, black, and other
retention_reg$Race <- ifelse(retention_reg$Race != "White" & retention_reg$Race != "African American", 
                             "Other", 
                             retention_reg$Race)
```


## First Trial on Linear Regression

### Linear Regression for Students Starting College in 2018

```{r}
# Filter students who started college in 2018
retention2018 <- retention_reg %>%
  filter(ENROLLMENT_BEGIN_year == 2018) #485 obs
```

```{r}
# Construct null and full model for linear regression 2018
md.null2018 = lm(retention ~ Race + EverReceivedPromiseAward + Race*EverReceivedPromiseAward, 
                 data = retention2018)
md.full2018 = lm(retention~AttendanceRate+Num_AP+Num_CTE+KeystoneMean+Race+Gender+ELLStatus+
                   IEPGroup+EconDisad+SAT_Total+CumulativeGPA+MagnetInd+
                   EverReceivedPromiseAward+QualifiedforCorePromise+semester+
                   Race*EverReceivedPromiseAward, data = retention2018)
```

```{r}
# Stepwise variable selection (both directions)
md.2018 = step(md.null2018, 
               list(lower = formula(md.null2018), 
                    upper = formula(md.full2018)), 
               direction = "both", trace = 0)
summary(md.2018)
```

```{r}
# Model diagnostics 
plot(md.2018)
```


### Linear Regression for Students Starting College in 2019

```{r}
retention2019 <- retention_reg %>%
  filter(ENROLLMENT_BEGIN_year == 2019) #663 obs
```

```{r}
# Construct null and full model for linear regression 2019
md.null2019 = lm(retention ~ Race + EverReceivedPromiseAward + Race*EverReceivedPromiseAward, 
                 data = retention2019)
md.full2019 = lm(retention~AttendanceRate+Num_AP+Num_CTE+KeystoneMean+Race+Gender+ELLStatus+
                   IEPGroup+EconDisad+SAT_Total+CumulativeGPA+MagnetInd+
                   EverReceivedPromiseAward+QualifiedforCorePromise+semester+
                   Race*EverReceivedPromiseAward, data = retention2019)
```

```{r}
# Stepwise variable selection (both directions)
md.2019 = step(md.null2019, 
               list(lower = formula(md.null2019), 
                    upper = formula(md.full2019)), 
               direction = "both", trace = 0)
summary(md.2019)
```

```{r}
# Model diagnostics 
plot(md.2019)
```


## Second Trial on Poisson Regression

### Poisson Regression for Students Starting College in 2018

```{r}
# Construct null possion model and full poisson model
poisson2018.null = glm(retention~Race+EverReceivedPromiseAward+Race*EverReceivedPromiseAward, 
                       family = "poisson", data = retention2018)
poisson2018.full = glm(retention~AttendanceRate+Num_AP+Num_CTE+KeystoneMean+Race+Gender+
                         ELLStatus+IEPGroup+EconDisad+SAT_Total+CumulativeGPA+MagnetInd+
                         EverReceivedPromiseAward+QualifiedforCorePromise+semester+
                         Race*EverReceivedPromiseAward, family = "poisson", data = retention2018)
```

```{r}
# Stepwise variable selection (both directions)
poisson2018 = step(poisson2018.full, 
                   list(lower = formula(poisson2018.null), 
                        upper = formula(poisson2018.full)), 
                   direction = "both", 
                   trace = 0)
summary(poisson2018)
```

```{r, message=FALSE}
# Model diagnostics 
# Residual deviance 
pchisq(poisson2018$deviance, 
       poisson2018$df.residual, 
       lower.tail = FALSE)

# Diagnostics plot
plot(poisson2018)

# Binned plot
x.2018 = predict(poisson2018)
y.2018 = residuals(poisson2018)
binnedplot(x.2018,y.2018)
```


### Poisson Regression for Students Starting College in 2019

```{r}
# Construct null possion model and full poisson model
poisson2019.null = glm(retention~Race+EverReceivedPromiseAward+Race*EverReceivedPromiseAward, 
                       family = "poisson", data = retention2019)
poisson2019.full = glm(retention~AttendanceRate+Num_AP+Num_CTE+KeystoneMean+Race+Gender+
                         ELLStatus+IEPGroup+EconDisad+SAT_Total+CumulativeGPA+MagnetInd+
                         EverReceivedPromiseAward+QualifiedforCorePromise+semester+
                         Race*EverReceivedPromiseAward, family = "poisson", data = retention2019)
```

```{r}
# Stepwise variable selection (both directions)
poisson2019 = step(poisson2019.full, list(lower = formula(poisson2019.null), 
                                          upper = formula(poisson2019.full)), 
                   direction = "both", trace = 0)
summary(poisson2019)
```

```{r, message=FALSE}
# Model diagnostics 
# Residual deviance 
pchisq(poisson2019$deviance, poisson2019$df.residual, lower.tail = FALSE)

# Diagnostics plot
plot(poisson2019)

# Binned plot
x.2019 = predict(poisson2019)
y.2019 = residuals(poisson2019)
binnedplot(x.2019,y.2019)
```


## Comparison between Poisson Regression and Linear Regression

### Retention for Students Starting College in 2018

```{r}
# Compare fitness of linear regression and poisson regression
# md.2018 for linear regression of 2018 retention
# poisson2018 for poisson regression of 2018 retention
anova(md.2018, poisson2018)
```


### Retention for Students Starting College in 2019

```{r}
# Compare fitness of linear regression and poisson regression
# md.2019 for linear regression of 2019 retention
# poisson2019 for poisson regression of 2019 retention
anova(md.2019, poisson2019)
```

**Conclusions:**

* Since the p-value is very small ($< 2.2\times 10^{-16}$), we say that poisson model fits better than linear regression model for 2018 retention.
* Since the p-value is also very small ($< 2.2\times 10^{-16}$), we say that poisson model fits better than linear regression model for 2019 retention.
* Thus, we will employ poisson regression model for students' retention in college.

\newpage

# Part III: Multivariate Regression Analysis for Box Range

## Data Wrangling

```{r}
# Current GPA and attendance Cutoff: GPA 2.5, attendance 90%
# Create different cutoffs for GPA and attendance 
gpa_upper = seq(2.9, 3.2, by = 0.1)
gpa_lower = seq(2.1, 1.8, by = -0.1)
attn_upper = seq(0.94, 1, by = 0.02)
attn_lower = seq(0.86,0.8, by = -0.02)
range_df = data.frame(attn_lower,attn_upper, gpa_lower,gpa_upper)
range_df
```

```{r}
# Clean data 
# retention_reg data is the final data set used for multivariate regression analysis 
# for whole range of GPA and attendance
retention_joint = read.csv("retention_reg.csv", stringsAsFactors = T) 
retention_joint$RandomID <- as.character(retention_joint$RandomID)
retention_joint <- retention_joint %>% 
  mutate(Gender = ifelse(Gender == "Female",1,0),
         ELLStatus = ifelse(ELLStatus == "ELL",1,0),
         EconDisad = ifelse(EconDisad == "Free Lunch",1,0),
         QualifiedforCorePromise = ifelse(QualifiedforCorePromise == "yes",1,0)) %>%
  dplyr::select(-c(X, AttendaceRateCate, GradYear, retention)) %>% 
  inner_join(nsc_retention, by = c("RandomID","ENROLLMENT_BEGIN_year", "semester")) %>%
  mutate(retention = as.numeric(retention))
head(retention_joint)
```

```{r}
# Prepare data sets for box range analysis
# gpa_scholarship_attendance data is the joined data set of gpa, scholarship, and attendance
joint_df = read.csv("gpa_sholarship_attendance.csv")
names(joint_df)[11] <- "present_pct"
names(joint_df)[12] <- "present_pct_cate"
joint_df$RandomID <- as.character(joint_df$RandomID)

box_df = joint_df %>% 
  dplyr::select(RandomID, CumulativeGPA, excused_pct) %>% 
  inner_join(nsc_promise,  by = "RandomID") %>% 
  filter(CumulativeGPA >=2, CumulativeGPA<=3) %>%
  filter(excused_pct >=0.85, excused_pct<=0.95)
```


## EDA

```{r}
# Visualize the box ranges 
ggplot(joint_df, aes(x =CumulativeGPA, y= excused_pct)) +
  geom_point() +
  labs(y = "Attendance") +
  geom_rect(aes(xmin = 2.1, xmax = 2.9, ymin = 0.86, ymax = 0.94), color = "red", alpha = 0)+
  geom_rect(aes(xmin = 2, xmax = 3, ymin = 0.84, ymax = 0.96), color = "red", alpha = 0)+
  geom_rect(aes(xmin = 1.9, xmax = 3.1, ymin = 0.82, ymax = 0.98), color = "red", alpha = 0)+
  geom_rect(aes(xmin = 1.8, xmax = 3.2, ymin = 0.8, ymax = 1), color = "red", alpha = 0) + 
  theme_bw() +
  labs(title = "Box Range for Treatment Analysis(2017-2020)")
```

```{r}
# Exclude semester and IEPGroup because of insufficient factor levels
get_box_df <- function(df, attn_lo, attn_high, gpa_lo, gpa_high){
  res = df %>%
    filter(AttendanceRate >=attn_lo, AttendanceRate<=attn_high) %>%
    filter(CumulativeGPA >=gpa_lo, CumulativeGPA<=gpa_high) 
  return(res)
}
```


## T-test for the treatment (EverReceivedPromiseAward)

```{r}
# T-test for 2018
t.test.summary.2018 <- data.frame(estimate=numeric(),
                 lower=numeric(), 
                 upper = numeric(),
                 pval=numeric(),
                 significant = logical()) 

# Conduct T-test for each box range 
year = 2018
for(ind in 1:4){
  row = unlist(range_df[ind,])
  box_df = get_box_df(retention_joint,row[1], row[2],row[3],row[4]) 
  box_df_year = box_df %>% filter(ENROLLMENT_BEGIN_year == year)
  a = t.test(retention~EverReceivedPromiseAward, data = box_df_year)
  res = data.frame(meanDiff=round(a$estimate[2]-a$estimate[1],3),
                 lower=round(a$conf.int[1],3), 
                 upper = round(a$conf.int[2],3),
                 pval=round(a$p.value,3), 
                 significant = a$p.value<=0.05) 
  t.test.summary.2018 = rbind(t.test.summary.2018, res)
  
}
rownames(t.test.summary.2018) <- paste0("Box", 1:4)
t.test.summary.2018
```

```{r}
# T-test for 2019
t.test.summary.2019 <- data.frame(estimate=numeric(),
                 lower=numeric(), 
                 upper = numeric(),
                 pval=numeric(),
                 significant = logical()) 

# Conduct T-test for each box range 
year = 2019
for(ind in 1:4){
  row = unlist(range_df[ind,])
  box_df = get_box_df(retention_joint,row[1], row[2],row[3],row[4]) 
  box_df_year = box_df %>% filter(ENROLLMENT_BEGIN_year == year)
  a = t.test(retention~EverReceivedPromiseAward, data = box_df_year)
  res = data.frame(meanDiff=round(a$estimate[2]-a$estimate[1],3),
                 lower=round(a$conf.int[1],3), 
                 upper = round(a$conf.int[2],3),
                 pval=round(a$p.value,3), 
                 significant = a$p.value<=0.05) 
  t.test.summary.2019 = rbind(t.test.summary.2019, res)
  
}
rownames(t.test.summary.2019) <- paste0("Box", 1:4)
t.test.summary.2019
```

**Findings:**

* We see that for 2018, the treatment effect is significant for all boxes except for the first one. This is probably because of limited observations in the smallest box (24 observations).
* We see that for 2019, none of the boxes has significant treatment effect. This aligns with our box plot for the whole data set above.


## Multivariate Linear Regression

* We perform linear regression on retention with all other predicting variables. 
* The linear regression serves as an extra step to verify the treatment effect of scholarship on retention.

```{r,warning=FALSE}
# Define null model and full model for multivariate linear regression
form.0 = as.formula("retention ~ Race + EverReceivedPromiseAward + Race*EverReceivedPromiseAward")
form.full = as.formula("retention ~ AttendanceRate + Num_AP + Num_CTE + KeystoneMean +
                       Race + Gender + ELLStatus + EconDisad + SAT_Total + CumulativeGPA +
                       MagnetInd + EverReceivedPromiseAward + QualifiedforCorePromise + 
                       Race*EverReceivedPromiseAward")

# Create a function box_yearly to perform multivariate linear regression for each box range 
box_yearly <- function(range_ind, year = 2018){
  row = unlist(range_df[range_ind,])
  box_df = get_box_df(retention_joint,row[1], row[2],row[3],row[4]) 
  box_df_year = box_df %>% filter(ENROLLMENT_BEGIN_year == year)
  print(paste("Number of students:", nrow(box_df_year)))
  lm.full.year = lm(form.full, data = box_df_year)
  res = step(lm.full.year, list(lower = form.0, upper = form.full), direction = "both", trace = 0)
  return(res)
} 
```

```{r}
# Construct linear regression model for each box range of 2018 retention
res.2018 = lapply(1:4, box_yearly)

# Model diagnostics for the last box
plot(res.2018[[4]])

# Coefficient output for the last box
(round((summary(res.2018[[4]]))$coef[,"Pr(>|t|)"],3))
```

```{r,warning = FALSE}
# Construct linear regression model for each box range of 2019 retention
res.2019 = lapply(1:4, box_yearly, year = 2019)

# Model diagnostics for the last box
plot(res.2019[[4]])

# Coefficient output for the last box
(round((summary(res.2019[[4]]))$coef[,"Pr(>|t|)"],3))
```

**Notice:**

* For simplicity, we only demonstrate the model diagnosis for the last box. 
* Replacing the index "4" with 1 through 3 to see model performance and coefficients for other 3 smaller boxes.
* More results of the linear regression will be discussed in the paper.



