final.Rmd

---
title: "Ravenswood Peer School Working File"
author: "Qingyang Zhang & Dorna Abdi"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Import data

For more data field references: https://caaspp-elpac.ets.org/caaspp/ResearchFileFormatSB?ps=true&lstTestYear=2022&lstTestType=B

```{r importing_packages, include=FALSE}
library(tidyverse)
library(tidymodels)
library(readxl)
library(dplyr)
library(ggplot2)
library(Hmisc) #https://CRAN.R-project.org/package=Hmisc
library(stargazer)
library(RANN) #https://cran.r-project.org/web/packages/RANN/RANN.pdf

library(stats) #hclust
library(FactoMineR)
library(factoextra)
library(cluster)
library(ggpubr)
library(corrplot)

library(caret) #KNN regression
library(glmnet) #LASSO regression
library(plyr)
library(reshape2)

set.seed(42)
```

```{r loading_data, include=FALSE}
caaspp22 <- read.table('data/caaspp.txt',sep='^', header=TRUE)
student_groups <- read_excel('data/lookup.xlsx', sheet = "student_groups")
entities <- read_csv('data/entities.csv')
#ravenswood <- read_csv('data/rvsw.csv') 
#bh <- read_csv('data/bh.csv') 
```

```{r exporting, include=FALSE}
#write.csv(entities,file='data/entities.csv')
#write.csv(caaspp22,file='data/caaspp22.csv')
#ravenswood <- caaspp22 %>% filter(District.Code=='68999')
#write.csv(ravenswood,file='data/rvsw.csv')
#bh <- caaspp22 %>% filter(School.Code=='6044309')
#write.csv(bh,file='data/bh.csv')
```

# Wrangling

<div class="alert alert-danger">
TODO:
<form>
  <input type="checkbox">
  <label> generate an equivalent dataset for Test.ID==2, i.e. math results </label><br>
  <input type="checkbox">
  <label> decide on grade list > type of schools </label><br>
  <input type="checkbox">
  <label> there are still 57 missing values in score when merging with demo </label><br>
</form>

Decisions:

* keep charter schools for now
* drop "declined to state" and keep only one variable for parental education
* use "Currently EL/or not" as the one variable for ELLs
</div>

```{r full-data}
select_cols <- c('School.Code',
                  'Student.Group.ID',
                  'Students.Tested')

# All Students	001
# Disability Status	Reported disabilities	128
# Economic Status	Economically disadvantaged	031
# English-Language Fluency EL (English learner)	160
# Race and Ethnicity	
## American Indian or Alaska Native	075, Asian 076, Black or African American 074
## Filipino	077, Hispanic or Latino	078, Native Hawaiian or Pacific Islander 079
## White 080, Two or more races	144
# Gender Female	004
# Parent Education
## Not a high school graduate	090
select_groups <- c(1, 128, 31, 160, 75, 76, 74, 77, 78, 79, 80, 144, 4, 90)

caaspp22_ela_pivoted <- caaspp22 %>% 
  # Grade 13 means all grades, Test.ID = 1 (ELA)
  filter(Grade==13 & Test.ID==1) %>%
  # public schools (7), Direct(9)/Locally(10) Funded Charter School
  filter(Type.ID==7 | Type.ID==9 | Type.ID==10) %>%
  select(all_of(select_cols)) %>%
  filter(Student.Group.ID %in% select_groups) %>%
  # merging on student group name
  left_join(y=student_groups, by = c("Student.Group.ID" = "Demographic ID")) %>% 
  rename("Demographic.ID.Num" = "Demographic ID Num",
         "Demographic.Name" = "Demographic Name") %>%
  select(all_of(c('School.Code',
                  'Students.Tested',
                  'Demographic.Name'))) %>%
  # pivot from long from to wide
  # each student group for a school is pivoted from a row to a val in column
  pivot_wider(names_from = Demographic.Name, values_from = 'Students.Tested')
```

After we pivot the original dataset, we have 10257 rows (public schools, direct funded charter school and locally funded) and 16 variables (school code, school size, and other 14 student group size). for the school size as well as the size of each student group, the value in the cell is simply the cardinality/count.

Abnormality:
* School Code `161242` (New Haven Unified District Level Program - Alameda New Haven Unified) is in english but not in math 
* School Code `1030873` (Miles W. Culwell Community Day - Fresno Coalinga-Huron Unified) is in math but not in english
* So in total 10258 schools, with only 10256 schools with both subject tested. 

## Select Sample

### Grade Spans
```{r school_type}
caaspp22_grade <- caaspp22 |> 
  group_by(School.Code) |> # 10,258 Unique School Codes
  filter(School.Code != 0 & Grade != 13) |>
  mutate(Grade.Count = n_distinct(Grade)) |>
  mutate(Grade.List = paste0(unique(Grade), collapse = ','))

grade_count <- caaspp22 |> 
  group_by(School.Code) |> # 10,258 Unique School Codes
  filter(School.Code != 0 & Grade != 13) |>
  summarise(Grade.Count = n_distinct(Grade))

grade_list <- caaspp22 |> 
  group_by(School.Code) |> # 10,258 Unique School Codes
  filter(School.Code != 0 & Grade != 13) |>
  summarise(Grade.List = paste0(unique(Grade), collapse = ','))

grades <- grade_count |>
  left_join(y=grade_list, by=c("School.Code" = "School.Code"))

school_type_table <- as.data.frame(
  table(grades$Grade.List, grades$Grade.Count)) |>
  setNames(c('grade_list', 'grade_count', 'freq')) |>
  filter(freq != 0)

# write.csv(school_type_table,file='data/school_type.csv')

typical_span <- c(
  '3,4,5', '3,4,5,6', '3,4,5,6,7,8', '3,4,5,6,7,8,11', '6,7,8', '6,7,8,11', '7,8', '7,8,11')

grade_atypical_span <- grades |>
  filter(!Grade.List %in% typical_span)

is_grade_atypical_span <- grade_atypical_span$School.Code
```

In our attempt to tackle the inconsistency in the combinations of grade levels, we eliminate the non-traditional schools because they are probably not representative of schools as a whole and are likely going through some type of transition (opening/closing/expanding). We consulted Austin from Ravenswood District, and decided to include schools with the following grade spans:

* 3,4,5
* 3,4,5,6
* 3,4,5,6,7,8
* 3,4,5,6,7,8,11
* 6,7,8
* 6,7,8,11
* 7,8
* 7,8,11
 
This includes 74% of the data (7629/10258 schools) and would give more confidence that schools are more directly comparable. 

### Outcome inspection

```{r scores-missing}
grade_non_missing <- caaspp22 |>
  # select 10258 schools (1 more since there is one school only have math scores)
  filter(School.Code > 0) |>
  # select only ELA results
  filter(Test.ID==1) |>
  # select `All Students`
  filter(Student.Group.ID == 1) |>
  filter(Grade != 13) |>
  # select rows where `Score` is '' or '*'
  filter(Mean.Scale.Score != '' & Mean.Scale.Score != '*') |> 
  group_by(School.Code) |>
  summarise(Non.Missing.Grade.Count = n_distinct(Grade), 
            Non.Missing.Grade.List = paste0(unique(Grade), collapse = ','))

# find schools where there are all missing grades 
grade_all_missing <- grade_non_missing |>
  right_join(y=grades, by = c("School.Code" = "School.Code")) |>
  filter(is.na(Non.Missing.Grade.Count))

# find schools where there is only one non-missing grade and it is grade 11
grade_all_missing_but_11 <- grade_non_missing |>
  right_join(y=grades, by = c("School.Code" = "School.Code")) |>
  filter(Non.Missing.Grade.Count == 1 & Grade.Count > 1 & Non.Missing.Grade.List == '11')

is_grade_all_missing <- grade_all_missing$School.Code
is_grade_all_missing_but_11 <- grade_all_missing_but_11$School.Code
```


### Missing values 

Noticeably there are a lot of missing values. For some of the student subgroups, the number of missing values are as high as more than 50%. In average, there are more than 2 missing values for every row. The specific distribution could be seen below.

Besides `NA` value, there are also a lot of `*` in the table. These two types of missing values are generated by two different mechanism:

* `NA`: comes from the pivoting process, i.e. there is not even a row for that specific student group for that specific school
* `*`: comes from the original dataset (the original dataset has no `NA` value) It is for privacy reasons, and theoretically indicates small values. 

```{r inspect-missing}
# NA for original data
caaspp22 %>% summarise_all(~ sum(is.na(.))) # NA in each column
table(rowSums(is.na(caaspp22))) # NA in each row

# NA for pivoted 
caaspp22_ela_pivoted %>% summarise_all(~ sum(is.na(.))) # NA in each column
table(rowSums(is.na(caaspp22_ela_pivoted))) # NA in each row

count_asterick <- function(x) length(which(x == '*'))

# * for original data
apply(caaspp22, 2, count_asterick) # * in each column
table(apply(caaspp22, 1, count_asterick)) # * in each row

# * for pivoted
apply(caaspp22_ela_pivoted, 2, count_asterick) # * in each column
table(apply(caaspp22_ela_pivoted, 1, count_asterick)) # * in each row
```

The way we handle the missing values involve two steps. For all `*` cells, we replace it with n=1 to show it is an arbitrarily small value. For all `NA` cells, we replace it with 0 assuming the reason why there is missing row is because the school does not have that type of students. 

In addition to replacing missing values by imputing, we also process them by dropping observations.

* we drop all 223 schools where 0 students, or 287 schools with `*` in `All Students` column
* we drop those schools with missing `Mean.Scale.Score` data  
* we drop those schools with atypical grade scope

```{r full-data-missing-val}
caaspp22_ela_filled <- data.frame(caaspp22_ela_pivoted) # duplicate 
  # merge on grade data
  # merge on outcome data

caaspp22_ela_filled <- caaspp22_ela_filled |>
  # drop schools with too few students
  filter(All.Students != 0 & All.Students != '*') |>
  # drop those schools with missing `Mean.Scale.Score` data  
  filter(!School.Code %in% is_grade_all_missing)  |> 
  filter(!School.Code %in% is_grade_all_missing_but_11)  |> 
  # drop those schools with atypical grade scope
  filter(!School.Code %in% is_grade_atypical_span) 
  

caaspp22_ela_filled[caaspp22_ela_filled=='*'] <- '1'    # *s replaced by 1
caaspp22_ela_filled[is.na(caaspp22_ela_filled)] <- '0'  # NA replaced by 0
caaspp22_ela_filled <- caaspp22_ela_filled %>% mutate_if(is.character, as.numeric)
```


```{r rename-export}
col_names <- names(caaspp22_ela_filled)[3:length(names(caaspp22_ela_filled))]
for (col_name in col_names) {
  # rename with the prefix - Percentage
  name <- paste0('Perc.', col_name)
  # calculate the percentage (%)
  caaspp22_ela_filled[name] <- caaspp22_ela_filled[col_name] / caaspp22_ela_filled$`All.Students` * 100
}
caaspp22_ela_demo <- caaspp22_ela_filled %>% select(-all_of(col_names))

write.csv(caaspp22_ela_demo,file='data/caaspp22_ela_demo.csv')
```

```{r demographic-table}
caaspp22_demo_summary <- caaspp22 |> 
  # Grade 13 means all grades, Test.ID = 1 (ELA)
  filter(Grade==13 & Test.ID==1 & School.Code == 0) |>
  filter(County.Code == 0 | District.Code == '68999') |>
  select(all_of(c('District.Code','Student.Group.ID','Students.Tested'))) |>
  filter(Student.Group.ID %in% select_groups) |>
  # merging on student group name
  left_join(y=student_groups, by = c("Student.Group.ID" = "Demographic ID")) |>
  rename("Demographic.ID.Num" = "Demographic ID Num",
         "Demographic.Name" = "Demographic Name") |> 
  select(all_of(c('District.Code',
                  'Students.Tested',
                  'Demographic.Name'))) |> 
  pivot_wider(names_from = Demographic.Name, values_from = 'Students.Tested') |>
  mutate_if(is.character, as.numeric)

colnames(caaspp22_demo_summary) <-colnames(caaspp22_ela_demo)
col_names <- names(caaspp22_demo_summary)[3:length(names(caaspp22_demo_summary))]
for (col_name in col_names) {
  # calculate the percentage (%)
  caaspp22_demo_summary[col_name] <- caaspp22_demo_summary[col_name] / caaspp22_demo_summary$`All.Students` * 100
}

ravenswood_school_codes <- caaspp22 %>% 
            filter(District.Code=='68999') %>% 
            #filter(Type.ID==7) %>%
            select(School.Code) %>% 
            unique() %>% 
            .$School.Code

caaspp22_ela_demo <- caaspp22_ela_demo |> 
  filter(School.Code %in% ravenswood_school_codes) 
caaspp22_demo_summary <- caaspp22_demo_summary|>
  bind_rows(caaspp22_ela_demo)|>
  rename('Code'='School.Code',
         'Size' = 'All.Students') |>
  column_to_rownames(var = 'Code') |>
  mutate_if(is.numeric, round, digits=2)

caaspp22_demo_summary[, 2:length(names(caaspp22_demo_summary))]<- paste0(as.matrix(caaspp22_demo_summary[, 2:length(names(caaspp22_demo_summary))]), '%')
caaspp22_demo_summary <- t(caaspp22_demo_summary)

#library(xtable)
#xtable(caaspp22_demo_summary)
```

## Outcome Encoding

```{r scores-distance}
#df by grade
caaspp_ela_outcome <- caaspp22 |>
  select(School.Code, Student.Group.ID, Test.ID, Grade, Students.Enrolled, Students.Tested, Mean.Scale.Score) |>
  # select 10258 schools (1 more since there is one school only have math scores)
  filter(School.Code > 0) |>
  # select only ELA results
  filter(Test.ID==1) |>
  # select `All Students`
  filter(Student.Group.ID == 1) |>
  # only for grade 3 - 8
  filter(Grade == 3 | Grade == 4 | Grade == 5 | Grade == 6 | Grade == 7 | Grade == 8) |>
  select(School.Code, Grade, Students.Enrolled, Students.Tested, Mean.Scale.Score)

#Grade 3 distance from score (2432)
grade_3 <- caaspp_ela_outcome |>
  filter(Grade == 3) |>
  mutate(distance_from_score_ela = as.numeric(Mean.Scale.Score) - 2432)

#Grade 4 distance from score (2473)
grade_4 <- caaspp_ela_outcome |>
  filter(Grade == 4) |>
  mutate(distance_from_score_ela = as.numeric(Mean.Scale.Score) - 2473)

#Grade 5 distance from score (2502)
grade_5 <- caaspp_ela_outcome |>
  filter(Grade == 5) |>
  mutate(distance_from_score_ela = as.numeric(Mean.Scale.Score) - 2502)

#Grade 6 distance from score (2618)
grade_6 <- caaspp_ela_outcome |>
  filter(Grade == 6) |>
  mutate(distance_from_score_ela = as.numeric(Mean.Scale.Score) - 2618)

#Grade 7 distance from score (2649)
grade_7 <- caaspp_ela_outcome |>
  filter(Grade == 7) |>
  mutate(distance_from_score_ela = as.numeric(Mean.Scale.Score) - 2649)

#Grade 8 distance from score (2668)
grade_8 <- caaspp_ela_outcome |>
  filter(Grade == 8) |>
  mutate(distance_from_score_ela = as.numeric(Mean.Scale.Score) - 2668)

caaspp_ela_score_distance <- rbind(grade_3, grade_4, grade_5, grade_6, grade_7, grade_8) |>
  na.omit()

caaspp_ela_aggregate <- caaspp_ela_score_distance |>
  group_by(School.Code) |>
  summarise_at(vars(distance_from_score_ela), list(name = mean)) |>
  rename("Avg.Score.Dist"="name")

```
## Output
```{r merge-output}

# output a file attaching score results to demo data
caaspp22_ela_score <- caaspp22_ela_demo |> left_join(y=caaspp_ela_aggregate, by = c("School.Code" = "School.Code"))

write.csv(caaspp22_ela_score, file='data/caaspp22_ela_score.csv')
```

```{r merge-both}
caaspp22_math_score <- read_csv('data/caaspp22_math_score.csv') |>
  select(c('School.Code','Avg.Score.Dist.Math'))
caaspp22_ela_score <- read_csv('data/caaspp22_ela_score.csv')
caaspp22_score <- caaspp22_ela_score |>
  inner_join(y=caaspp22_math_score, by = c("School.Code" = "School.Code")) |>
  select(-1)

write.csv(caaspp22_score, file='data/caaspp22_score.csv')
```

# Modeling


## Pre-processing
Before running models, we would like to conduct the proper preprocessing. 

```{r importing-preprocessing}

# 7291 Schools with the non-trivial size, proper grade span, and non-missing grades
caaspp22_final <- read_csv('data/caaspp22_ela_demo.csv') |>
    # drop the column to set School.Code as index
    column_to_rownames(var = 'School.Code') |>
    select(-'...1') |> 
    # standardize all columns
    mutate(across(where(is.numeric), scale))

# no missing values
# caaspp22_ela[is.na(caaspp22_ela), ]

caaspp22_score <- read_csv('data/caaspp22_score.csv')
```
Since scaling a variable is a linear transformation and it will not change the distribution of the variable so it does not matter if the variable has a non-normal distribution, we don't have to concern about the outlier School, City of Angels at Los Angeles Unified, which is way greater than all other schools (n=6421>>3763, the second large). 

## Nearest Neighbor

For NN, we use the `nn2` function from `RANN` package, which finds the k nearest neighbours for every point in a given dataset in O(N log N) time using Arya and Mount's ANN library (v1.1.3). The distance is computed using the L2 (Euclidean) metric. For alternative metric, one could see package 'RANN.L1' for the same
functionality using the L1 (Manhattan, taxicab) metric.

Notice that the built-in knn function runs a k-nearest neighbour classification for test set from training set, which does not apply for our situation. 

```{r modeling}
#NN search: https://search.r-project.org/CRAN/refmans/RANN/html/nn2.html

# k: The maximum number of nearest neighbours to compute. 
# to generate top 10 nearest neighbours we set k = 11 to offset oneself
nearest <- nn2(caaspp22_final, caaspp22_final, k=11)
```

```{r modeling-loop-all}
ravenswood_school_codes <- caaspp22 %>% 
                                filter(District.Code=='68999') %>% 
                                filter(Type.ID==7 | Type.ID==9 | Type.ID==10) %>%
                                select(School.Code) %>% 
                                unique() %>% 
                                .$School.Code

top10_ravenswood <- data.frame(Rank=numeric(0), 
                               Dist=numeric(0),
                               Center=numeric(0),
                               School.Code=numeric(0))

for (school in ravenswood_school_codes) {
  
  # each ravenswood school has ten rows, i.e. top ten neighbors 
  row_idx <- which(rownames(caaspp22_final) == as.character(school)) 
  top10_id <- nearest$nn.idx[row_idx,]
  top10_ds <-nearest$nn.dists[row_idx,]
  
  # each row a neighbor with info about rank，dist, school.code
  for (rk in 1:11) {
    nn_idx <- top10_id[rk]
    neighbor_code <- rownames(caaspp22_final)[nn_idx]
    top10_ravenswood[nrow(top10_ravenswood) + 1, ] <- c(
      rk-1, top10_ds[rk], school, neighbor_code
    )  
  }
  
}

top10_ravenswood <- top10_ravenswood |> 
  mutate_if(is.character, as.numeric) |>
  left_join(y=entities, by = c("School.Code" = "School Code")) |>
  left_join(y=caaspp22_score, by = c("School.Code" = "School.Code"))

write.csv(top10_ravenswood, file='data/top10_ravenswood.csv')

```

## KNN regression
```{r KNN_regression}

#KNN regression: https://www.rdocumentation.org/packages/caret/versions/6.0-86/topics/knnreg
#knnreg {caret}: https://search.r-project.org/CRAN/refmans/caret/html/knnreg.html

df <- readr::read_csv('data/caaspp22_score.csv') |>
    # drop the column to set School.Code as index
    tibble::column_to_rownames(var = 'School.Code') |>
    # scale every column
    dplyr::mutate(across(2:length(df)-2, scale))
#df[c(length(df)-1, length(df))] <- lapply(df[c(length(df)-1, length(df))], log)

# ravenswood results as test set 
test <- df |>
  dplyr::filter(rownames(df) %in% ravenswood_school_codes)
test_x = test[, -c(1, length(df)-1, length(df))]
test_ela = test[, length(df)-1]
test_math = test[, length(df)]

# create train and validation set 
df <- df |> 
  dplyr::filter(!rownames(df) %in% ravenswood_school_codes) 
train <- df |> 
  dplyr::sample_frac(0.85)
val  <- dplyr::anti_join(df, train, by ='...1')

train_x = train[, -c(1, length(df)-1, length(df))]
train_ela = train[, length(df)-1]
train_math = train[, length(df)]

val_x = val[, -c(1, length(df)-1, length(df))]
val_ela = val[, length(df)-1]
val_math = val[, length(df)]
```

```{r KNN_regression_k}

# K-fold validation for optimal k: https://towardsdatascience.com/beginners-guide-to-k-nearest-neighbors-in-r-from-zero-to-hero-d92cd4074bdb

nfold = 10
set.seed(1)

# cut() divides the range into several intervals
folds = seq.int(nrow(df)) %>%
     cut(breaks = nfold, labels=FALSE) %>%  
     sample

x = df[, -c(1, length(df)-1, length(df))]
y_ela = df[, length(df)-1]
y_math = df[, length(df)]

do.chunk <- function(chunkid, folddef, Xdat, Ydat, k) { 
    train = (folddef!=chunkid)# training index
    Xtr = Xdat[train,] # training set by the index
    Ytr = Ydat[train] # true label in training set
    Xvl = Xdat[!train,] # test set
    Yvl = Ydat[!train] # true label in test set
    knnmodel = knnreg(Xtr, Ytr, k=k)
    pred_y = predict(knnmodel, data.frame(Xvl))
    rmse = caret::RMSE(Yvl, pred_y)
    data.frame(fold=chunkid, val.error = rmse) # test error per fold
}

# set error.folds to save validation errors
error.folds.ela=NULL
error.folds.math=NULL
# create a sequence of data for k
kvec = seq(5,80,2)

for (j in kvec) {
   # apply do.function to each fold
   tmp_ela = ldply(1:nfold, do.chunk, 
               folddef=folds, Xdat=x, Ydat=y_ela, k=j) # required arguments
   tmp_math = ldply(1:nfold, do.chunk, 
               folddef=folds, Xdat=x, Ydat=y_math, k=j)
   # track each value of neighbors 
   tmp_ela$neighbors = j
   tmp_math$neighbors = j
   error.folds.ela = rbind(error.folds.ela, tmp_ela)
   error.folds.math = rbind(error.folds.math, tmp_math)
}
```

```{r }

val.error.means = error.folds.ela %>%
    #group the selected data by neighbors
    group_by(neighbors) %>%
    #cacluate CV error for each k
    summarise_each(funs(mean), val.error) %>%
    #remove existing grouping
    ungroup() %>% 
    filter(val.error==min(val.error))
# Best number of neighbors
# if there is a tie, pick larger number of neighbors for simpler model
numneighbor = max(val.error.means$neighbors)
numneighbor

val.error.means = error.folds.math %>%
    #group the selected data by neighbors
    group_by(neighbors) %>%
    #cacluate CV error for each k
    summarise_each(funs(mean), val.error) %>%
    #remove existing grouping
    ungroup() %>% 
    filter(val.error==min(val.error))
# Best number of neighbors
# if there is a tie, pick larger number of neighbors for simpler model
numneighbor = max(val.error.means$neighbors)
numneighbor
```

```{r KNN_regression-ELA}
### ELA

# train the KNN model 
knnmodel = knnreg(x, y_ela, k=15)

# predict on ravenswood schools
pred_y = predict(knnmodel, data.frame(test_x))
knn_ela <- data.frame(test_ela, pred_y, test_ela-pred_y)

mse = mean((test_ela - pred_y)^2)
mae = caret::MAE(test_ela, pred_y)
rmse = caret::RMSE(test_ela, pred_y)

cat("MSE: ", mse, "MAE: ", mae, " RMSE: ", rmse)

```

```{r KNN_Math}
### Math

# train the KNN model 
knnmodel_math = knnreg(x, y_math, k=15)

# predict on ravenswood schools
pred_y = predict(knnmodel_math, data.frame(test_x))
knn_math <- data.frame(test_math, pred_y, test_math-pred_y)

mse = mean((test_math - pred_y)^2)
mae = caret::MAE(test_math, pred_y)
rmse = caret::RMSE(test_math, pred_y)

#cat("MSE: ", mse, "MAE: ", mae, " RMSE: ", rmse)
```

## OLS & LASSO

One limitation of KNN is its lack in explanability. To complement it, we use other methods to help understand. 

First we conduct regular multivariate linear regression.

```{r OLS}

# multivariate OLS
model_ela <- lm(Avg.Score.Dist ~ ., data = df[, -c(1, length(df))])
summary(model_ela)

model_math <- lm(Avg.Score.Dist.Math ~ ., data = df[, -c(1, length(df)-1)])
summary(model_math)
```

Then we use a LASSO regression to conduct feature selection. 

```{r LASSO_ELA}
mod_cv <- cv.glmnet(x=as.matrix(df[, -c(1, length(df)-1, length(df))]), 
                    y=df[, length(df)-1], 
                    family='gaussian',
                    alpha=1)

# cvm : The mean cross-validated error - a vector of length length(lambda)
# lambda.min : the λ at which the minimal MSE is achieved.
# lambda.1se : the largest λ at which the MSE is within one standard error of the minimal MSE.
   
plot(log(mod_cv$lambda), mod_cv$cvm)
plot(mod_cv) 
coef(mod_cv, c(mod_cv$lambda.min,
               mod_cv$lambda.1se))
print(paste(mod_cv$lambda.min,
            log(mod_cv$lambda.min)))
print(paste(mod_cv$lambda.1se,
            log(mod_cv$lambda.1se)))


best_model_ela <- glmnet(x=as.matrix(df[, -c(1, length(df)-1, length(df))]), 
                    y=df[, length(df)-1], 
                    family='gaussian',
                    lambda = mod_cv$lambda.1se,
                    alpha=1)

plot(mod_cv$glmnet.fit, 
     "lambda", label=TRUE)
```
```{r LASSO_plot_ELA}
X = as.matrix(df[, -c(1, length(df)-1, length(df))])
colnames(X) <- c("Size", "%Female", "%FRPL", "%Black", "%American.Indian",
                 "%Asian", "%Hispanic", "%White", "%Not.HS", "%IEP", 
                 "%2+R", "%ELL", "%Filipino", "%Native.H")  
mod <- glmnet(x=X, 
                    y=df[, length(df)-1], 
                    family='gaussian',
                    alpha=1)

library(plotmo) # for plot_glmnet
#png('lasso.png', res=300)
plot_glmnet(mod, label=5)
#dev.off()

```

```{r LASSO_MATH}
mod_cv <- cv.glmnet(x=as.matrix(df[, -c(1, length(df)-1, length(df))]), 
                    y=df[, length(df)], 
                    family='gaussian',
                    alpha=1)

# cvm : The mean cross-validated error - a vector of length length(lambda)
# lambda.min : the λ at which the minimal MSE is achieved.
# lambda.1se : the largest λ at which the MSE is within one standard error of the minimal MSE.
   
plot(log(mod_cv$lambda), mod_cv$cvm)
plot(mod_cv) 
coef(mod_cv, c(mod_cv$lambda.min,
               mod_cv$lambda.1se))
print(paste(mod_cv$lambda.min,
            log(mod_cv$lambda.min)))
print(paste(mod_cv$lambda.1se,
            log(mod_cv$lambda.1se)))
```


```{r LASSO_plot_MATH}
X = as.matrix(df[, -c(1, length(df)-1, length(df))])
colnames(X) <- c("Size", "%Female", "%FRPL", "%Black", "%American.Indian",
                 "%Asian", "%Hispanic", "%White", "%Not.HS", "%IEP", 
                 "%2+R", "%ELL", "%Filipino", "%Native.H")  
mod <- glmnet(x=X, 
                    y=df[, length(df)], 
                    family='gaussian',
                    alpha=1)

library(plotmo) # for plot_glmnet
#png('lasso.png', res=300)
plot_glmnet(mod, label=5)
#dev.off()

```

```{r KNN_reduced}
drop_by_LASSO <- c('Perc.Hispanic.or.Latino', 'Perc.White', 'Perc.Two.or.more.races', 'Perc.Filipino')
caaspp22_final_drop <- caaspp22_final |>
  select(-drop_by_LASSO)
```

## K-means Clustering
```{r pca}
# scale. = FALSE/TRUE will return different proportion of variance explained
# but the contributions of each variable look alike
res.pca <- prcomp(caaspp22_final, center = FALSE, scale. = FALSE)
summary(res.pca)

pcaCharts <- function(x) {
    x.var <- x$sdev ^ 2
    x.pvar <- x.var/sum(x.var)
    print("proportions of variance:")
    print(x.pvar)
    
    par(mfrow=c(2,2))
    plot(x.pvar,xlab="Principal component", ylab="Proportion of variance explained", ylim=c(0,1), type='b')
    plot(cumsum(x.pvar),xlab="Principal component", ylab="Cumulative Proportion of variance explained", ylim=c(0,1), type='b')
    screeplot(x)
    screeplot(x,type="l")
    par(mfrow=c(1,1))
    corrplot(res.pca$rotation, is.corr=TRUE)
    corrplot(res.pca$rotation, is.corr=FALSE)
}

pcaCharts(res.pca)

fviz_pca_var(res.pca,
             col.var = "contrib", # Color by contributions to the PC
             gradient.cols = c("snow3", "darkgoldenrod2", "#8c1415ff"),
             repel = TRUE,    # Avoid text overlapping
             )

ggsave("pca.png", width = 6, height = 6)

caaspp22_pca = as.data.frame(-res.pca$x[,1:2])

fviz_cos2(res.pca, choice = "var", axes = 1:3)
fviz_contrib(res.pca, choice = "var", axes = 1:3, top = 10)
```


```{r WORKING}
#fig.width=5,fig.height=5

peer_school_codes <- top10_ravenswood %>% 
  filter(Rank != 0) %>% 
  filter(Center == 6044309) %>%
  .$School.Code

#cols <- c("snow3", "#55AD89", "#EF6F6A")
cols <- c("snow3", "darkgoldenrod2", "#8c1415ff")

caaspp22_pca <- read_csv('data/caaspp22_ela_demo.csv') |>
    # drop the column to set School.Code as index
    filter('School.Code' != '1996115') |>
    column_to_rownames(var = 'School.Code') |>
    select(-'...1') |>
    mutate(across(-c('All.Students'), scale))
  
normalize <- function(x) {
  x <- log(x + 1e-10)
  return ((x - min(x)) / (max(x) - min(x)) * 100) 
}
caaspp22_pca['All.Students'] <- lapply(caaspp22_pca['All.Students'], normalize)
    

caaspp22_pca_biplot <- caaspp22_pca |>
  rownames_to_column('School.Code') %>%
  mutate(School.Group = ifelse(
                 School.Code==6044309, "Ravenswood - Belle Haven",
                 ifelse(School.Code %in% peer_school_codes, "Peer School", "Other")))

#1
ggplot(caaspp22_pca_biplot, 
       aes(x = Perc.Economically.disadvantaged, 
           y = Perc.Not.a.high.school.graduate,
           color = School.Group
          )) +
  geom_point(#color="cornflowerblue", 
             #size = 1, 
             aes(alpha = School.Group,
                 size= School.Group)) +
  geom_text(aes(label=ifelse(School.Code==6044309,'','')),
            hjust=.2,
            vjust=.2) + 
  scale_color_manual(values = cols) + 
  scale_size_discrete(c(0.2, 0.5, 1.5)) + 
  scale_alpha_discrete(c(0.1, 0.5, 1)) +
  labs(x = "Perc.Economically.disadvantaged",
       y = "Perc.Not.a.high.school.graduate") +
  theme_minimal()

#2
ggplot(caaspp22_pca_biplot, 
       aes(x = All.Students, 
           y = Perc.EL..English.learner.,
           color = School.Group
          )) +
  geom_point(#color="cornflowerblue", 
             #size = 1, 
             aes(alpha = School.Group,
                 size= School.Group)) +
  geom_text(aes(label=ifelse(School.Code==6044309,'','')),
            hjust=.2,
            vjust=.2) + 
  scale_color_manual(values = cols) + 
  scale_size_discrete(c(0.2, 0.5, 1.5)) + 
  scale_alpha_discrete(c(0.1, 0.5, 1)) +
  labs(x = "All.Students",
       y = "Perc.EL..English.learner.") +
  theme_minimal()

#4
ggplot(caaspp22_pca_biplot, 
       aes(x = Perc.Hispanic.or.Latino,                                                            
           y = Perc.White,
           color = School.Group
          )) +
  geom_point(#color="cornflowerblue", 
             #size = 1, 
             aes(alpha = School.Group,
                 size= School.Group)) +
  geom_text(aes(label=ifelse(School.Code==6044309,'','')),
            hjust=.2,
            vjust=.2) + 
  scale_color_manual(values = cols) +
  scale_size_discrete(c(0.2, 0.5, 1.5)) + 
  scale_alpha_discrete(c(0.1, 0.5, 1)) +
  labs(x = "Perc.Hispanic.or.Latino",
       y = "Perc.White") +
  theme_minimal()
```

```{r kmeans_nclusters}
fviz_nbclust(caaspp22_pca, 
             kmeans, 
             method='wss' # the elbow shows up when k = 4
             #method='silhouette' # similarly suggest k = 3
             )
```

```{r kmeans}
res.km3 <- kmeans(caaspp22_pca, 3)
res.km4 <- kmeans(caaspp22_pca, 4)
#concern: not exactly same size
#print(res.km4)
fviz_cluster(res.km4, data=caaspp22_final, 
             geom="point")

#It’s possible to compute the mean of each variables by clusters using the original data:
#aggregate(USArrests, by=list(cluster=km.res$cluster), mean)
```

```{r kmeans_eclust}
res.km3.eclust <- eclust(caaspp22_final, "kmeans", k = 3,
                 nstart = 25, graph = FALSE)
res.km4.eclust <- eclust(caaspp22_final, "kmeans", k = 4,
                 nstart = 25, graph = FALSE)

fviz_cluster(res.km3.eclust,
            geom = "point",
            ellipse.type = "norm",
            ellipse.level = 0.95,
            ellipse.alpha = 0.2,
            pointsize = 0.8,
            ggtheme=theme_bw())
             
fviz_cluster(res.km4.eclust,
            geom = "point",
            ellipse.type = "norm",
            ellipse.level = 0.95,
            ellipse.alpha = 0.2,
            pointsize = 0.8,
            ggtheme=theme_bw())
             
fviz_silhouette(res.km3.eclust)
fviz_silhouette(res.km4.eclust)
```