
In this lab you will respond to a set of prompts for two parts.
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:
Change your outcome to students’ final exam performance (note: check the data dictionary for a pointer!).
Using the same data (and testing and training data sets), recipe, and workflow as you used in the case study, change the mode of your model from classification to regression and change the engine from a glm to an lm model.
Interpret your regression machine learning model in terms of three regression machine learning model metrics: MAE, MSE, and RMSE. Read about these metrics here. Similar to how we interpreted the classification machine learning metrics, focus on the substantive meaning of these statistics.
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:
MAE: A MAE of 195.16 indicates that the predicted scores are about 195 points off the actual on average.
MSE: An MSE of 65089.45 indicates the average squared differences between actual and predicted
RMSE: An RMSE of 255.12 indicates the average magnitude of errors between the predicted and actual. None of these results indicate a good model unfortunately.
Complete the following steps to knit and publish your work:
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.
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.
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!
