---
title: 'Technical Appendix B: Scholarship Analysis'
output:
  html_document:
    df_print: paged
---

# Data Wrangling

```{r, message=FALSE, warning=FALSE}
# Library packages
library(tidyverse)
library(dplyr)
library(visdat)
library(stats)
library(ggplot2)
library(DHARMa)
library(arm)
library(readxl)
```

```{r}
# Import data
scholarship <- read_excel("DR-826 Pittsburgh Promise scholarship use.xlsx", sheet = "Sheet1")
NSC <- read_excel("DR-844 NSC data_not including 1920.xlsx", sheet = "NSC data")
sat <- read_excel("DR-834 SAT scores.xlsx", sheet = "SAT Scores")
ap <- read_excel("DR-834 AP scores.xlsx", sheet = "AP Exam Scores")
keystone <- read_excel("DR-836 Keystone assessment results.xlsx", sheet = "Keystone data")
enrollment <- read_excel("DR-828 Enrollment.xlsx", sheet = "Enrollment records")
```

```{r}
# Clean and reformat the data above 
# Create new data sets for data compilation later

# Create magnet_model data
# Serve as indicator of whether students attend magnet school
enrollment_new <- enrollment %>%
  mutate_at(c(2:14), as.factor) %>%
  mutate(EntryDate=as.Date(as.character(EntryDate), tryFormats = c("%Y%m%d"))) %>%
  mutate(WithdrawDate=as.Date(as.character(WithdrawDate), tryFormats = c("%Y%m%d")))

magnet_model <- enrollment_new %>%
  filter(FullMagnetInd == "1") %>%
  distinct(RandomID, EnrolledSchoolID, SchoolReportName, FullMagnetInd)
magnet_model$RandomID <- as.character(magnet_model$RandomID)

# Create keystone_model data
# Provide information about students' keystone scores 
keystone_new <- keystone %>%
  mutate_at(c(1,3:7), as.factor) %>%
  mutate(AdminScaleScore=as.numeric(AdminScaleScore))

keystone_model <- keystone_new %>%
  na.omit() %>%
  group_by(RandomID)%>%
  dplyr::summarise(ScoreMean = mean(AdminScaleScore))
keystone_model$RandomID <- as.character(keystone_model$RandomID)

# Create sat_model data 
# Provide information about each student's highest SAT score
sat_new <- sat %>%
  mutate_at(2, as.factor) %>%
  mutate(LATEST_ASSESSMENT_DATE=as.Date(LATEST_ASSESSMENT_DATE)) %>%
  mutate_at(c(4:6), as.numeric)

sat_rep <- sat_new$RandomID[duplicated(sat_new$RandomID)]

sat_model <- sat_new %>% 
  distinct(RandomID, Latest_SAT_Total, .keep_all = TRUE)
sat_model$RandomID <- as.character(sat_model$RandomID)

# Create nsc_model data
# Joined data set between NSC and scholarship
nsc2 <- NSC %>%
  dplyr::select(-c("Cohort")) %>%
  mutate(GradYear = substr(as.character(HIGH_SCHOOL_GRAD_DATE),1,4)) %>%
  mutate(EnrollYear = substr(as.character(ENROLLMENT_BEGIN),1,4)) %>%
  filter(EnrollYear == GradYear) 

nsc2$RandomID <- as.character(nsc2$RandomID)
scholarship$GradYear <- as.character(scholarship$GradYear)
names(scholarship)[1] <- "RandomID"

nsc_scholarship <- inner_join(nsc2, scholarship, by = c("RandomID", "GradYear"))

nsc_model <- nsc_scholarship %>%
  dplyr::select(-c("HIGH_SCHOOL_GRAD_DATE",
            "ENROLLMENT_BEGIN",
            "ENROLLMENT_END")) %>%
  distinct()
nsc_model$RandomID <- as.character(nsc_model$RandomID)
```

```{r}
# Import data sets generated before 
# Please refer to Appendix A for more details in code

# attendance_rate is the joined data set between attendance and scholarship
attendance_model <- read.csv("attendance_rate.csv")
attendance_model$RandomID <- as.character(attendance_model$RandomID)

# cte_scholarship is the joined data set between CTE and scholarship
cte_model <- read.csv("cte_scholarship.csv")
cte_model$RandomID <- as.character(cte_model$RandomID)

# senior_gpa is the cumulative GPA data generated from the original GPA data
senior_gpa <- read.csv("senior_gpa.csv")
senior_gpa$RandomID <- as.character(senior_gpa$RandomID)

# demo_scholarship_final is the joined data set between demographics and scholarship
demographics_model <- read.csv("demo_scholarship_final.csv")
demographics_model$RandomID <- as.character(demographics_model$RandomID)

# ap_scholarship is the joined data set between ap and scholarship
ap_model <- read.csv("ap_scholarship.csv")
ap_model$RandomID <- as.character(ap_model$RandomID)
```

```{r}
# Data cleaning and reformatting 
# Left join all data sets
data_join <- attendance_model %>%
  left_join(ap_model, by = c("RandomID"="RandomID")) %>%
  left_join(cte_model, by = c("RandomID"="RandomID")) %>%
  left_join(keystone_model, by = c("RandomID"="RandomID")) %>%
  left_join(demographics_model, by = c("RandomID"="RandomID")) %>%
  left_join(sat_model, by = c("RandomID"="RandomID")) %>%
  left_join(senior_gpa, by = c("RandomID"="RandomID")) %>%
  left_join(magnet_model, by = c("RandomID"="RandomID")) 

# Select variables by column names
# Delete duplicated rows
data_variables <- data_join %>%
  dplyr::select(RandomID, GradYear.x, excused_pct, excused_pct_cate, Num_AP, num_cte, 
         ScoreMean, Race, Gender, ELLStatus, IEPGroup, EconDisad, Latest_SAT_Total, 
         CumulativeGPA, FullMagnetInd, EverReceivedPromiseAward.x) %>%
  distinct()

# Deal with NAs
# Replace missing num_AP, num_cte, FullMagnetInd with 0
data_na <- data_variables %>%
  replace_na(list(Num_AP=0, num_cte=0, FullMagnetInd=0))

# Visualize NAs
vis_miss(data_na) 

# Omit all NAs
data_omit <- data_na %>%
  na.omit() 

# Change response variable to 0,1
data_omit$EverReceivedPromiseAward.x <- as.numeric(as.factor(data_omit$EverReceivedPromiseAward.x))-1

# Change variables into appropriate format
data_everReceived <- data_omit %>%
  mutate_at(c(2, 4,8:12,15:16), as.factor)

# Rename columns
colnames(data_everReceived)<- c("RandomID","GradYear","AttendanceRate", "AttendaceRateCate", 
                                "Num_AP", "Num_CTE", "KeystoneMean", "Race", "Gender", 
                                "ELLStatus", "IEPGroup", "EconDisad", "SAT_Total", 
                                "CumulativeGPA","MagnetInd", "EverReceivedPromiseAward")
data_everReceived$Race_new <- ifelse(data_everReceived$Race != "White" & data_everReceived$Race != "African American", "Others", as.character(data_everReceived$Race))
data_everReceived$Race_new <- as.factor(data_everReceived$Race_new)
```

\newpage 

# Part I: Logistic Regression Analysis for All Students 

```{r}
# Construct null model and full model 
model_null <- glm(EverReceivedPromiseAward ~ Race_new+Gender, data = data_everReceived, 
                  family = "binomial")
model_full <- glm(EverReceivedPromiseAward ~ AttendanceRate+Num_AP+Num_CTE+KeystoneMean+Race_new
              +Gender+ELLStatus+IEPGroup+EconDisad+SAT_Total+CumulativeGPA+MagnetInd, 
              data = data_everReceived, family = "binomial")
```

```{r}
# Stepwise selection on AIC
# Backwards selection
backwards <- step(model_full, scope=list(lower=formula(model_null),upper=formula(model_full)), 
                  direction="backward", trace = 0)

# Forwards selection 
forwards <- step(model_null, scope=list(lower=formula(model_null), upper=formula(model_full)), 
                 direction="forward", trace = 0)

# Selection on both Directions
bothways <- step(model_null, list(lower=formula(model_null), upper=formula(model_full)), 
                 direction="both", trace=0)

# Summary of the selected models
# Choose the selected model on both directions as our final model
summary(bothways)
formula(backwards)
formula(forwards)
formula(bothways)
```

**Interpretations for all significant variables:**

* When holding everything else fixed:      
  + the odds to receive Promise for white students are 30.6% lower than black students;    
  + the odds to receive Promise for student who are in magnet school are 38.3% higher than those in non-magnet school;      
  + the odds for student who are not in ELL group are 200.7% higher than those in ELL group;      
  + for 0.1 increase in GPA, we expect 9.4% increase in the odds of receiving Promise;      
  + for 0.01 increase in attendance rate, we expect 8.4% increase in the odds of receiving Promise;       
  + for 100 increase in Keystone mean, we expect 55.5% decrease in the odds of receiving Promise.      
  
```{r}
# Model diagnostics 
plot(bothways)

y <- resid(bothways)
x <- predict(bothways)
binnedplot(x,y)

sim <- simulateResiduals(bothways)
plot(sim)
```

**Conclusion of model diagnostics:**

* From the binned plot, we know that the selected model is a good fit. 
  + gray lines refer to the 95% confidence interval;
  + most points fall within the confidence interval;
  + residuals are randomly scattered with no obvious pattern. 
* From QQ residuals plot, the residuals follow normal distribution.
* From residuals vs predicted plot, variance assumption of residuals is not violated by the selected model. 

\newpage
 
# Part II: Logistic Regression Analysis for Qualified Students 

```{r}
# Prepare data sets for the analysis on only qualified students 
data_qualified <- merge(data_everReceived, 
                        attendance_model[, c("RandomID", "QualifiedforCorePromise")], by="RandomID")

data_qualifiedYes <- data_qualified %>%
  filter(QualifiedforCorePromise == "yes")

data_qualifiedYes$Race_new <- ifelse(data_qualifiedYes$Race != "White" & data_qualifiedYes$Race != "African American", "Other", as.character(data_qualifiedYes$Race))
data_qualifiedYes$Race_new <- as.factor(data_qualifiedYes$Race_new)
data_qualifiedYes$QualifiedforCorePromise <- as.factor(data_qualifiedYes$QualifiedforCorePromise)
```

```{r}
# Construct null model and full model
model_null <- glm(EverReceivedPromiseAward ~ Race_new+Gender, data = data_qualifiedYes, family = "binomial")
model_full <- glm(EverReceivedPromiseAward ~ AttendanceRate+Num_AP+Num_CTE+KeystoneMean+Race_new
              +Gender+ELLStatus+IEPGroup+EconDisad+SAT_Total+CumulativeGPA+MagnetInd, 
              data = data_qualifiedYes, family = "binomial")
```

```{r}
# Stepwise selection on AIC
# Backwards selection 
backwards <- step(model_full, scope=list(lower=formula(model_null), upper=formula(model_full)), 
                  direction="backward", trace = 0)
 
# Forward selection 
forwards <- step(model_null, scope=list(lower=formula(model_null),upper=formula(model_full)), 
                 direction="forward", trace = 0)

# Selection on both directions
bothways <- step(model_null, list(lower=formula(model_null),upper=formula(model_full)), 
                 direction="both", trace=0)

# Summary of selected model
# Choose the selected model on both directions as our final model
summary(bothways)
formula(backwards)
formula(forwards)
formula(bothways)
```

**Interpretation for all significant variables:**

* When holding everything else fixed:      
  + the odds to receive Promise for student who are in magnet school are 25.3% higher than those in non-magnet school;      
  + the odds for student who are not in ELL group are 119% higher than those in ELL group;       
  + for 0.1 increase in GPA, we expect 4.9% increase in the odds of receiving Promise;      
  + for 0.01 increase in attendance rate, we expect 11% increase in the odds of receiving Promise;      
  + for 100 increase in Keystone mean, we expect 57.9% decrease in the odds of receiving Promise.      

```{r}
# Model diagnostics 
plot(bothways)

y <- resid(bothways)
x <- predict(bothways)
binnedplot(x,y)

sim <- simulateResiduals(bothways)
plot(sim)
```

**Conclusion of model diagnostics:**

* From the binned plot, we know that the selected model is a good fit. 
  + gray lines refer to the 95% confidence interval;
  + most points fall within the confidence interval;
  + residuals are randomly scattered with no obvious pattern. 
* From QQ residuals plot, the residuals follow normal distribution.
* From residuals vs predicted plot, variance assumption of residuals is not violated by the selected model. 












