In this lab you will respond to a set of prompts for two parts.

Part I: Data Product

For the data product, you will interpret a different type of model – a model in a regression mode.

So far, we have specified and interpreted a classification model: one predicting a dichotomous outcome (i.e., whether students pass a course). In many cases, however, we are interested in predicting a continuous outcome (e.g., students’ number of points in a course or their score on a final exam).

While many parts of the machine learning process are the same for a regression machine learning model, one key part that is relevant to this lab is different: their interpretation. The confusion matrix we created to parse the predictive strength of our classification model does not pertain to regression machine learning models. Different metrics are used. For this lab, you will specify and interpret a regression machine learning model.

The requirements are as follows:

Please use the code chunk below for your code:

#clear the memory
rm(list = ls())

setwd("/Users/nickponzi/lab/data")
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.5.0     ✔ purrr   1.0.2
## ✔ tibble  3.2.1     ✔ dplyr   1.1.4
## ✔ tidyr   1.3.1     ✔ stringr 1.5.0
## ✔ readr   2.1.3     ✔ forcats 0.5.2
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(tidymodels)
## Warning: package 'tidymodels' was built under R version 4.2.3
## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
## ✔ broom        1.0.5      ✔ rsample      1.2.1 
## ✔ dials        1.2.1      ✔ tune         1.2.0 
## ✔ infer        1.0.7      ✔ workflows    1.1.4 
## ✔ modeldata    1.3.0      ✔ workflowsets 1.1.0 
## ✔ parsnip      1.2.1      ✔ yardstick    1.3.1 
## ✔ recipes      1.0.10
## Warning: package 'dials' was built under R version 4.2.3
## Warning: package 'scales' was built under R version 4.2.3
## Warning: package 'infer' was built under R version 4.2.3
## Warning: package 'modeldata' was built under R version 4.2.3
## Warning: package 'parsnip' was built under R version 4.2.3
## Warning: package 'recipes' was built under R version 4.2.3
## Warning: package 'rsample' was built under R version 4.2.3
## Warning: package 'tune' was built under R version 4.2.3
## Warning: package 'workflows' was built under R version 4.2.3
## Warning: package 'workflowsets' was built under R version 4.2.3
## Warning: package 'yardstick' was built under R version 4.2.3
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ scales::discard() masks purrr::discard()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ recipes::fixed()  masks stringr::fixed()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step()   masks stats::step()
## • Search for functions across packages at https://www.tidymodels.org/find/
library(janitor)
## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(yardstick)

# Load the assessments file
assessments <- read_csv("oulad-assessments.csv")
## Rows: 173912 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): code_module, code_presentation, assessment_type
## dbl (7): id_assessment, id_student, date_submitted, is_banked, score, date, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#Load the students file
students <- read.csv("oulad-students.csv")

# Count the number of assessments of different types
assessments %>% 
  count(assessment_type)
## # A tibble: 3 × 2
##   assessment_type     n
##   <chr>           <int>
## 1 CMA             70527
## 2 Exam             4959
## 3 TMA             98426
# Count assessments per course
assessments %>% 
  count(assessment_type, code_module, code_presentation)
## # A tibble: 41 × 4
##    assessment_type code_module code_presentation     n
##    <chr>           <chr>       <chr>             <int>
##  1 CMA             BBB         2013B              5049
##  2 CMA             BBB         2013J              6416
##  3 CMA             BBB         2014B              4493
##  4 CMA             CCC         2014B              3920
##  5 CMA             CCC         2014J              5846
##  6 CMA             DDD         2013B              5252
##  7 CMA             FFF         2013B              6681
##  8 CMA             FFF         2013J              8847
##  9 CMA             FFF         2014B              5549
## 10 CMA             FFF         2014J              8915
## # ℹ 31 more rows
# Summarize assignment due dates
assessments %>% 
  summarize(mean_date = mean(date, na.rm = TRUE),
            median_date = median(date, na.rm = TRUE),
            sd_date = sd(date, na.rm = TRUE),
            min_date = min(date, na.rm = TRUE),
            max_date = max(date, na.rm = TRUE))
## # A tibble: 1 × 5
##   mean_date median_date sd_date min_date max_date
##       <dbl>       <dbl>   <dbl>    <dbl>    <dbl>
## 1      131.         129    78.0       12      261
# Find the first quantile value for assignment due dates
code_module_dates <- assessments %>% 
  group_by(code_module, code_presentation) %>% 
  summarize(quantile_cutoff_date = quantile(date, probs = .25, na.rm = TRUE))
## `summarise()` has grouped output by 'code_module'. You can override using the
## `.groups` argument.
# Join assessments and code_module_dates
assessments_joined <- left_join(assessments, code_module_dates, by = c("code_module", "code_presentation"))

# Filter assessments to include only assignments before the cutoff date
assessments_filtered <- assessments_joined %>% 
  filter(date < quantile_cutoff_date)

# Summarize average score for each student
assessments_summarized <- assessments_filtered %>% 
  mutate(weighted_score = score * weight) %>% 
  group_by(id_student) %>% 
  summarize(mean_weighted_score = mean(weighted_score))

# Process students data
students <- students %>% 
  mutate(pass = ifelse(final_result == "Pass", 1, 0)) %>% 
  mutate(pass = as.factor(pass))

students <- students %>% 
  mutate(imd_band = factor(imd_band, levels = c("0-10%", "10-20%", "20-30%", "30-40%", "40-50%", "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"))) %>% 
  mutate(imd_band = as.integer(imd_band))

# Join students and assessments_summarized
students_and_assessments <- left_join(students, assessments_summarized, by = "id_student")

# Split data into training and testing sets
set.seed(20230712)
train_test_split <- initial_split(students_and_assessments, prop = .50, strata = "pass")
data_train <- training(train_test_split)
data_test <- testing(train_test_split)

# Specify the recipe
my_rec <- recipe(pass ~ disability +
                   date_registration + 
                   gender +
                   code_module +
                   mean_weighted_score, 
                 data = data_train) %>% 
  step_dummy(disability) %>% 
  step_dummy(gender) %>%  
  step_dummy(code_module)

# Specify the model
my_mod <-
  logistic_reg() %>% 
  set_engine("glm") %>% 
  set_mode("classification")

# Specify the workflow
my_wf <-
  workflow() %>% 
  add_model(my_mod) %>% 
  add_recipe(my_rec)

# Specify the metric set (excluding kappa)
class_metrics <- metric_set(accuracy, sensitivity, specificity, precision, npv)

# Fit the model and collect metrics
final_fit <- last_fit(my_wf, train_test_split, metrics = class_metrics)

# Confusion matrix
conf_matrix <- collect_predictions(final_fit) %>% 
  conf_mat(.pred_class, pass)

# Access the confusion matrix table
conf_table <- conf_matrix$table

# Extract true positives, true negatives, false positives, and false negatives
TP <- conf_table[2, 2]  # True positives
TN <- conf_table[1, 1]  # True negatives
FP <- conf_table[1, 2]  # False positives
FN <- conf_table[2, 1]  # False negatives

# Calculate accuracy
accuracy <- ((TP + TN) / (TP + TN + FP + FN)) * 100

# Calculate sensitivity
sensitivity <- (TP / (TP + FN)) * 100

# Calculate specificity
specificity <- (TN / (FP + TN)) * 100

# Calculate precision (positive predictive value)
precision <- (TP / (TP + FP)) * 100

# Calculate negative predictive value
npv <- (TN / (TN + FN)) * 100

# Calculate Cohen's Kappa
kappa <- conf_matrix$kappa

# Repeat accuracy to match the length of other metrics
accuracy_vector <- rep(accuracy, length(class_metrics))


# Output all metrics
confusion_matrix_output <- data.frame(
  Statistic = c("Accuracy", "Sensitivity", "Specificity", "Precision", "Negative Predictive Value"),
  Value = c(accuracy_vector, sensitivity, specificity, precision, npv),
  Percentage = c(accuracy_vector, sensitivity, specificity, precision, npv),
  stringsAsFactors = FALSE
)


# Output confusion matrix and metrics
confusion_matrix_output
##                   Statistic    Value Percentage
## 1                  Accuracy 57.14979   57.14979
## 2               Sensitivity 39.40930   39.40930
## 3               Specificity 72.43300   72.43300
## 4                 Precision 55.18845   55.18845
## 5 Negative Predictive Value 58.11790   58.11790
# Specify the recipe for regression
my_rec_regression <- recipe(mean_weighted_score ~ disability +
                              date_registration + 
                              gender +
                              code_module,
                            data = data_train) %>% 
  step_dummy(disability) %>% 
  step_dummy(gender) %>%  
  step_dummy(code_module)

# Specify the model for regression
my_mod_regression <- 
  linear_reg() %>% 
  set_engine("lm") %>% 
  set_mode("regression")

# Specify the workflow for regression
my_wf_regression <- 
  workflow() %>% 
  add_model(my_mod_regression) %>% 
  add_recipe(my_rec_regression)

# Fit the regression model
final_fit_regression <- fit(my_wf_regression, data = data_train)

# Predict on training data
predicted_values <- predict(final_fit_regression, new_data = data_train)

# Remove rows with missing values
complete_cases <- complete.cases(data_train$mean_weighted_score, predicted_values$.pred)
data_train_complete <- data_train[complete_cases, ]
predicted_values_complete <- predicted_values[complete_cases, ]

# Calculate Mean Absolute Error (MAE) on complete cases
mae <- mean(abs(data_train_complete$mean_weighted_score - predicted_values_complete$.pred))

# Calculate Root Mean Squared Error (RMSE) on complete cases
rmse <- sqrt(mean((data_train_complete$mean_weighted_score - predicted_values_complete$.pred)^2))


# Remove rows with missing values
complete_cases <- complete.cases(data_train$mean_weighted_score, predicted_values$.pred)
data_train_complete <- data_train[complete_cases, ]
predicted_values_complete <- predicted_values[complete_cases, ]

# Calculate Mean Squared Error (MSE) on complete cases
mse <- mean((data_train_complete$mean_weighted_score - predicted_values_complete$.pred)^2)

# Output the MAE, MSE, and RMSE
mae
## [1] 195.1671
mse
## [1] 65089.45
rmse
## [1] 255.1263

Please add your interpretations here:

Part II: Reflect and Plan

  1. What is an example of an outcome related to your research interests that could be modeled using a classification machine learning model?
  1. What is an example of an outcome related to your research interests that could be modeled using a regression machine learning model?
  1. Look back to the study you identified for the first machine learning lab badge activity. Was the outcome one that is modeled using a classification or a regression machine learning model? Identify which mode(s) the authors of that paper used and briefly discuss the appropriateness of their decision.

Knit and Publish

Complete the following steps to knit and publish your work:

  1. First, change the name of the author: in the YAML header at the very top of this document to your name. The YAML header controls the style and feel for knitted document but doesn’t actually display in the final output.

  2. Next, click the knit button in the toolbar above to “knit” your R Markdown document to a HTML file that will be saved in your R Project folder. You should see a formatted webpage appear in your Viewer tab in the lower right pan or in a new browser window. Let’s us know if you run into any issues with knitting.

  3. Finally, publish your webpage on RPubs by clicking the “Publish” button located in the Viewer Pane after you knit your document. See screenshot below.

Have fun!