high_school_proj.Rmd

---
title: "High School Performance"
author: "Xin Bu"
date: "11/14/2023"
output: html_document
html_document:
    toc: true
    toc_depth: 6
    number_sections: true
    toc_float: true
    code_folding: hide
    theme: flatly
    code_download: true
---
-------------

### 1. Introduction

The data was collected from three high schools in the US, consisting of information on the students' performance in math, reading, and writing, which are continuous outcome variables. The students' demographic information including gender, race/ethnicity, parental level of education, lunch type, test preparation course are categorical variables as well as predictors in this project. The purpose of this project is to visualize the data, build a multi-linear regression model, apply a 10-fold cross-validation resampling method, and evaluate LM, SVM, and KNN model performance.


```{r setup, include=FALSE}
library(tidyverse,warn.conflicts = F)
library(tidymodels,warn.conflicts = F)
library(lubridate,warn.conflicts = F)
library(dplyr)
library(here)
library(table1)
library(kknn)
library(kernlab)
library(gcookbook)
library(patchwork)

knitr::opts_chunk$set(echo = TRUE)

```

#### 1.1. Load data
```{r model-fit1, message=F}
hs <- read_csv("high_school.csv")
head(hs)
```

#### 1.2. Rename the columns
```{r message=F,warning=F}

hs <- hs %>% 
    rename("gender" = "gender", "Eth" = "race/ethnicity", "Par_edu" = "parental level of education", "Lunch" = "lunch", "Test_prep" = "test preparation course", "Math" = "math score", "Reading" = "reading score", "Writing" = "writing score")

```


#### 1.3. Reorder the categorical variable Par_edu

```{r}

hs$Par_edu <- factor(hs$Par_edu, levels=c("some high school","high school","some college","associate's degree","bachelor's degree","master's degree"))

```


### 2. Data Visualization

#### 2.1. Visualize the percentage of parental education level by race/ethnicity
```{r}
ggplot(hs, mapping = aes(x = Eth, fill = Par_edu))+
  geom_bar(position = "fill")+
  scale_fill_brewer(palette="Set3")+
  scale_y_continuous(labels=percent)+
  labs(x = "Race/ethnicity", y = " ", title= "Parental Education Level by Race/Ethnicity")+
  theme(plot.title = element_text(hjust = 0.5, vjust = 3)) 

```

#### 2.2. Visualize math, reading, and writing by gender

```{r}

plot1 <- ggplot(hs,aes(x= Math, fill = gender, color = gender))+
  geom_density(alpha = 0.5,mapping=(aes(y = after_stat(scaled))))+
  scale_fill_manual(values = c("#E46726","#6D9EC1"))

plot2 <-ggplot(hs,aes(x= Reading, fill = gender, color = gender))+
  geom_density(alpha = 0.5,mapping=(aes(y = after_stat(scaled))))+
  scale_fill_manual(values = c("#E46726","#6D9EC1"))

plot3 <-ggplot(hs,aes(x= Writing, fill = gender, color = gender))+
  geom_density(alpha = 0.5,mapping=(aes(y = after_stat(scaled))))+
  scale_fill_manual(values = c("#E46726","#6D9EC1"))

plot1+plot2+plot3+
  plot_layout(ncol = 2)
```

#### 2.3. Visualize the relationship between math and reading by parental education level and lunch type

```{r}

ggplot(data = hs) +
  geom_point(mapping = aes(x = Math, y = Reading, color = Lunch, shape = Lunch), size=2)+ 
  scale_color_hue(l=45, c=60)+
  facet_wrap(~Par_edu) 

```


#### 2.4. Visualize math, reading, and writing for those who completed the test preparation
```{r}
test_prep <- hs %>%
  filter(Test_prep == "completed") %>%
  select(Math,Reading,Writing)

pairs(test_prep,
        col = c("#E46726","#6D9EC1"),
        pch = 21,
        main = "Pairwise Correlation")
```

#### 2.5. Visualize the count of test preparation (completed or none) by race/ethnicity
```{r}

ggplot(hs, aes(x = Eth, fill = Test_prep)) + 
  geom_bar(alpha = 1, position = "dodge2")+
  labs(x="Race/ethnicity", y=" ", title = "Test Preparation by Race/Ethnicity")+
  scale_fill_brewer(palette=c("YlOrRd"))+
  theme(plot.title = element_text(hjust = 0.5, vjust = 3)) 

```

#### 2.6. Visualize math by lunch type and gender based on parental education level
```{r}
ggplot(data = hs) +
  geom_boxplot(mapping = aes(x = reorder(Par_edu, Math, FUN = median), y = Math, fill=gender)) +
  facet_wrap(~Lunch)+ 
  coord_flip()+
  scale_fill_brewer(palette="RdPu")+
  labs(x="", y="Math", title = "Math by gender and Parental Education Level")+
  theme(plot.title = element_text(hjust = 0.5, vjust = 3)) 
```

### 3. Descriptive statistics 
#### 3.1. Math, reading, and writing related to each categorical variable
```{r}
table1::label(hs$gender) <= "gender"
table1::label(hs$Eth)<="Race/Ethnicity"
table1::label(hs$Lunch) <="Lunch"
table1::label(hs$Par_edu)<="Parental Level of Education"
table1::label(hs$Test_prep) <= "Test Preparation"
table1::label(hs$Math)<= "Math"
table1::label(hs$Reading) <= "Reading"
table1::label(hs$Writing) <= "Writing"

table1::table1(~Math+Reading+Writing|Eth, data=hs)
table1::table1(~Math+Reading+Writing|Lunch, data=hs)
table1::table1(~Math+Reading+Writing|gender,data=hs)
table1::table1(~Math+Reading+Writing|Test_prep, data=hs)
table1::table1(~Math+Reading+Writing|Par_edu, data=hs)

```

### 4. Model Building
#### 4.1. Build a model predicting math based on the categorical variables

```{r model-fit2, message=F}


lm.1 <- lm(Math ~ Eth + Lunch + Par_edu + Test_prep + gender, data=hs)

mod.summary <- summary(lm.1)

coefs <- mod.summary$coefficients[,1]
coef.se <- mod.summary$coefficients[,2] 
rse <- mod.summary$sigma


for.table <- tibble(term=rownames(mod.summary$coefficients),
                    coefficient=coefs,
                    se=coef.se)

knitr::kable(for.table,digits=2)

```

#### 4.2. Split a training dataset and build the model recipe
```{r}
set.seed(20231116)

splits <- initial_split(hs,strata=Par_edu)
cv_folds <- vfold_cv(training(splits))

high_school.rec <- 
  recipe(Math ~ Eth + Lunch + gender + Test_prep,
         data=training(splits)) %>% 
  step_dummy(all_nominal()) %>% 
  step_normalize(all_predictors()) %>% 
  prep()

```

#### 4.3. Build a linear regression model, support vector machine model, and k nearest neighbors model
```{r}
high_school.linear <- linear_reg() %>% 
  set_engine("lm")

high_school.svm <- svm_rbf(
  cost = tune(),
  rbf_sigma = tune()) %>% 
  set_engine("kernlab") %>% 
  set_mode("regression")

high_school.knn <- nearest_neighbor(
  neighbors = tune()) %>% 
  set_engine("kknn") %>% 
  set_mode("regression")

lin.wf <- workflow() %>% 
  add_recipe(high_school.rec) %>% 
  add_model(high_school.linear)

svm.wf <- workflow() %>% 
  add_recipe(high_school.rec) %>% 
  add_model(high_school.svm)

knn.wf <- workflow() %>% 
  add_recipe(high_school.rec) %>% 
  add_model(high_school.knn)
```

#### 4.4. Create different values for tuning the model parameters and get results of tuning those models
```{r cache = T}

knn.grid <- grid_regular(
  neighbors(range = c(1,50)),
  levels=10
)

svm.grid <- grid_regular(
  cost(),
  rbf_sigma(),
  levels=4)

knn.results <- 
  knn.wf %>% 
  tune_grid(
    resamples=cv_folds,
    grid=knn.grid)

svm.results <- 
  svm.wf %>% 
  tune_grid(
    resamples=cv_folds,
    grid=svm.grid)


if(F){
  knn.results %>% collect_metrics() %>% filter(.metric=="rsq") %>% arrange(mean)
  svm.results %>% collect_metrics() %>% filter(.metric=="rmse") %>% arrange(mean)
}

```

#### 4.5. Finalize workflows, fit the training data, use each model to predict, and then compare the results
```{r}

svm.final <- svm.wf %>% 
  finalize_workflow(svm.results %>% select_best("rmse")) %>% 
  fit(training(splits))

knn.final <- knn.wf %>% 
  finalize_workflow(knn.results %>% select_best("rmse")) %>% 
  fit(training(splits))

lin.final <- lin.wf %>% 
  fit(training(splits))


high_school <- testing(splits)

high_school <- high_school %>% 
  mutate(svm_pred = predict(svm.final,new_data=high_school) %>% pull(.pred)) %>% 
  mutate(knn_pred = predict(knn.final,new_data=high_school) %>% pull(.pred)) %>% 
  mutate(lin_pred = predict(lin.final,new_data=high_school) %>% pull(.pred)) 

holder <- metrics(high_school,truth=Math,estimate=svm_pred) %>% 
  select(-.estimator) %>% 
  rename(metric = .metric,
         svm = .estimate)

holder <- holder %>% 
  mutate(knn = 
    metrics(high_school,truth=Math,estimate=knn_pred) %>% 
      pull(.estimate)
  ) %>% 
  mutate(linear = 
    metrics(high_school,truth=Math,estimate=lin_pred) %>% 
      pull(.estimate)
  )

holder %>% 
  pivot_longer(-metric) %>% 
  ggplot(aes(x=name,y=value,group=1)) + 
  geom_point() + 
  facet_wrap(~metric,scales="free_y") + 
  theme_bw() +
  labs(x="Model",y="Metric")


```