Final Project
2024-05-20
Introduction
The aim of this project is to predict which variables best predict whether an individual earns more than $50,000 per year. I will be using the Adult Census Income dataset, pulled from the 1994 US Census Database by Barry Becker. It’s made up of information such as occupation, age, native country, race, capital gain, capital loss, education, work class, and more. The focal point of this project is to predict which variables most effectively determine whether an individual earns more than $50,000 per year, framing it as a classification problem where the outcome is binary: earning <=50K or >50K per year.
Various machine learning models will be explored and evaluated to handle the classification task, including logistic regression, decision trees, and more complex approaches like random forests and elastic net models. I start with essential steps such as handling missing data, normalizing data, and encoding categorical variables to prepare the data for effective model training and then finish with using ROC AUC as a metric to assess model performance, focusing on the model’s ability to discriminate between the two income classes accurately.
ROC AUC is a widely used metric to evaluate the performance of binary classification models. It stands for “Receiver Operating Characteristic - Area Under the Curve.” The curve is created by plotting the True Positive Rate (TPR) against the False Positive Rate (FPR) at various threshold settings.
Data Citation
Jain, Aruna. Adult Income Census Dataset. Kaggle, https://www.kaggle.com/datasets/jainaru/adult-income-census-dataset. Accessed 28 May 2024.
Codebook
age Type: Numeric Description: This variable represents the age of the individual. The values range from a minimum of 17 to a maximum of 90 years, with a mean age of approximately 38.58 years.
workclass Type: Categorical (character) Description: This variable represents the work classification of the individual. It includes categories such as ‘Private’, ‘Self-emp-not-inc’, ‘Local-gov’, etc. This variable has various class labels indicating the sector of employment.
fnlwgt Type: Numeric Description: The final weight is the number of units in the target population that the responding unit represents. It is used to adjust the dataset to be representative of the general population.
education Type: Categorical (character) Description: This variable represents the highest level of education attained by the individual. Categories include ‘Bachelors’, ‘HS-grad’, ‘11th’, ‘Masters’, ‘9th’, ‘Some-college’, ‘Assoc-acdm’, ‘Assoc-voc’, etc.
education.num Type: Ordered Categorical (integer) Description: This variable provides a numerical representation of the educational attainment. It ranges from 1 to 16, with 1 representing the lowest level of education and 16 representing the highest.
marital.status Type: Categorical (character) Description: This variable indicates the marital status of the individual. Categories include ‘Married-civ-spouse’, ‘Divorced’, ‘Never-married’, ‘Separated’, ‘Widowed’, ‘Married-spouse-absent’, and ‘Married-AF-spouse’.
occupation Type: Categorical (character) Description: This variable describes the type of occupation of the individual. Categories include ‘Tech-support’, ‘Craft-repair’, ‘Other-service’, ‘Sales’, ‘Exec-managerial’, ‘Prof-specialty’, ‘Handlers-cleaners’, etc.
relationship Type: Categorical (character) Description: This variable indicates the individual’s relationship within the household. Categories include ‘Wife’, ‘Own-child’, ‘Husband’, ‘Not-in-family’, ‘Other-relative’, and ‘Unmarried’.
race Type: Categorical (character) Description: This variable represents the race of the individual. Categories include ‘White’, ‘Asian-Pac-Islander’, ‘Amer-Indian-Eskimo’, ‘Other’, and ‘Black’.
sex Type: Categorical (character) Description: This variable indicates the sex of the individual. The categories are ‘Male’ and ‘Female’.
capital.gain Type: Numeric Description: This variable represents the capital gains earned by the individual from investment sources. The values range from 0 to 99,999.
capital.loss Type: Numeric Description: This variable represents the capital losses incurred by the individual from investment sources. The values range from 0 to 4,356.
hours.per.week Type: Numeric Description: This variable indicates the number of hours the individual works per week. The values range from 1 to 99 hours, with a mean of 40.44 hours per week.
native.country Type: Categorical (character) Description: This variable represents the native country of the individual. Categories include ‘United-States’, ‘Cambodia’, ‘England’, ‘Puerto-Rico’, ‘Canada’, ‘Germany’, ‘Outlying-US(Guam-USVI-etc)’, etc.
income Type: Categorical (character) Description: This is the target variable indicating whether the individual’s income is above or below $50,000 per year. The categories are ‘>50K’ and ‘<=50K’.
We start by loading in all the necessary libraries for the remainder of our project.
library(readr)
library(tidymodels)## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──## ✔ broom 1.0.5 ✔ recipes 1.0.10
## ✔ dials 1.2.1 ✔ rsample 1.2.1
## ✔ dplyr 1.1.4 ✔ tibble 3.2.1
## ✔ ggplot2 3.5.1 ✔ tidyr 1.3.1
## ✔ infer 1.0.7 ✔ tune 1.2.1
## ✔ modeldata 1.3.0 ✔ workflows 1.1.4
## ✔ parsnip 1.2.1 ✔ workflowsets 1.1.0
## ✔ purrr 1.0.2 ✔ yardstick 1.3.1## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ 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(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ scales::col_factor() masks readr::col_factor()
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ stringr::fixed() masks recipes::fixed()
## ✖ dplyr::lag() masks stats::lag()
## ✖ yardstick::spec() masks readr::spec()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(dplyr)
library(tidyr)
library(ggplot2)
library(corrplot)## corrplot 0.92 loadedlibrary(forcats)
library(naniar)
library(yardstick) # for area under ROC curve for log model
library(themis) # for upsampling
library(vip) #for important predictors##
## Attaching package: 'vip'
##
## The following object is masked from 'package:utils':
##
## vilibrary(ranger)
library(discrim)##
## Attaching package: 'discrim'
##
## The following object is masked from 'package:dials':
##
## smoothnesslibrary(poissonreg)
library(broom)
library(parsnip)
library(ranger)
tidymodels_prefer() #for the select function
income_data <- read_csv("C:/Users/roxy/Documents/school/Spring Qrtr -2024/pstat 131/pstat 131 final projecgt/adult_income_cleaned.csv")## Rows: 32561 Columns: 15
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (9): workclass, education, marital.status, occupation, relationship, rac...
## dbl (6): age, fnlwgt, education.num, capital.gain, capital.loss, hours.per.week
##
## ℹ 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.Exploratory Data Analysis
We will poke around the data set and see if we can identify any patterns. We are mainly trying to display basic information and summary statistics.
# Display basic information and summary statistics
income_data %>%
head()## # A tibble: 6 × 15
## age workclass fnlwgt education education.num marital.status occupation
## <dbl> <chr> <dbl> <chr> <dbl> <chr> <chr>
## 1 90 <NA> 77053 HS-grad 9 Widowed <NA>
## 2 82 Private 132870 HS-grad 9 Widowed Exec-manager…
## 3 66 <NA> 186061 Some-college 10 Widowed <NA>
## 4 54 Private 140359 7th-8th 4 Divorced Machine-op-i…
## 5 41 Private 264663 Some-college 10 Separated Prof-special…
## 6 34 Private 216864 HS-grad 9 Divorced Other-service
## # ℹ 8 more variables: relationship <chr>, race <chr>, sex <chr>,
## # capital.gain <dbl>, capital.loss <dbl>, hours.per.week <dbl>,
## # native.country <chr>, income <chr>colnames(income_data)## [1] "age" "workclass" "fnlwgt" "education"
## [5] "education.num" "marital.status" "occupation" "relationship"
## [9] "race" "sex" "capital.gain" "capital.loss"
## [13] "hours.per.week" "native.country" "income"dim(income_data)## [1] 32561 15str(income_data)## spc_tbl_ [32,561 × 15] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ age : num [1:32561] 90 82 66 54 41 34 38 74 68 41 ...
## $ workclass : chr [1:32561] NA "Private" NA "Private" ...
## $ fnlwgt : num [1:32561] 77053 132870 186061 140359 264663 ...
## $ education : chr [1:32561] "HS-grad" "HS-grad" "Some-college" "7th-8th" ...
## $ education.num : num [1:32561] 9 9 10 4 10 9 6 16 9 10 ...
## $ marital.status: chr [1:32561] "Widowed" "Widowed" "Widowed" "Divorced" ...
## $ occupation : chr [1:32561] NA "Exec-managerial" NA "Machine-op-inspct" ...
## $ relationship : chr [1:32561] "Not-in-family" "Not-in-family" "Unmarried" "Unmarried" ...
## $ race : chr [1:32561] "White" "White" "Black" "White" ...
## $ sex : chr [1:32561] "Female" "Female" "Female" "Female" ...
## $ capital.gain : num [1:32561] 0 0 0 0 0 0 0 0 0 0 ...
## $ capital.loss : num [1:32561] 4356 4356 4356 3900 3900 ...
## $ hours.per.week: num [1:32561] 40 18 40 40 40 45 40 20 40 60 ...
## $ native.country: chr [1:32561] "United-States" "United-States" "United-States" "United-States" ...
## $ income : chr [1:32561] "<=50K" "<=50K" "<=50K" "<=50K" ...
## - attr(*, "spec")=
## .. cols(
## .. age = col_double(),
## .. workclass = col_character(),
## .. fnlwgt = col_double(),
## .. education = col_character(),
## .. education.num = col_double(),
## .. marital.status = col_character(),
## .. occupation = col_character(),
## .. relationship = col_character(),
## .. race = col_character(),
## .. sex = col_character(),
## .. capital.gain = col_double(),
## .. capital.loss = col_double(),
## .. hours.per.week = col_double(),
## .. native.country = col_character(),
## .. income = col_character()
## .. )
## - attr(*, "problems")=<externalptr>summary(income_data)## age workclass fnlwgt education
## Min. :17.00 Length:32561 Min. : 12285 Length:32561
## 1st Qu.:28.00 Class :character 1st Qu.: 117827 Class :character
## Median :37.00 Mode :character Median : 178356 Mode :character
## Mean :38.58 Mean : 189778
## 3rd Qu.:48.00 3rd Qu.: 237051
## Max. :90.00 Max. :1484705
## education.num marital.status occupation relationship
## Min. : 1.00 Length:32561 Length:32561 Length:32561
## 1st Qu.: 9.00 Class :character Class :character Class :character
## Median :10.00 Mode :character Mode :character Mode :character
## Mean :10.08
## 3rd Qu.:12.00
## Max. :16.00
## race sex capital.gain capital.loss
## Length:32561 Length:32561 Min. : 0 Min. : 0.0
## Class :character Class :character 1st Qu.: 0 1st Qu.: 0.0
## Mode :character Mode :character Median : 0 Median : 0.0
## Mean : 1078 Mean : 87.3
## 3rd Qu.: 0 3rd Qu.: 0.0
## Max. :99999 Max. :4356.0
## hours.per.week native.country income
## Min. : 1.00 Length:32561 Length:32561
## 1st Qu.:40.00 Class :character Class :character
## Median :40.00 Mode :character Mode :character
## Mean :40.44
## 3rd Qu.:45.00
## Max. :99.00There are 15 columns. That means 15 different variables. These variables are: ‘age’, ‘fnlwgt’, ‘education.num’, ‘marital.status’, ‘occupation’, ‘relationship’, ‘race’, ‘sex’, ‘capital.gain’, ‘capital.loss’, ‘hours.per.week’, ‘native.country’, and lastly our predictor variable: ‘income’.
Converting variables into factors is crucial for ensuring proper data handling during analysis. For income, it ensures that the variable is correctly interpreted as categorical with two distinct groups, crucial for classification tasks. For education.num, turning it into an ordered factor acknowledges the inherent order in educational levels, which can be important for analyses that explore trends in education in relation to other variables. This approach enhances the accuracy and interpretability of your statistical models and analysis results.
# turning income and education.num into a factor variable
income_data$income <- factor(income_data$income, levels = c(">50K", "<=50K"))
income_data$education.num <- factor(income_data$education.num,
levels = 1:16,
ordered = TRUE)
#check levels to confirm
class(income_data$income)## [1] "factor"class(income_data$education.num)## [1] "ordered" "factor"levels(income_data$income)## [1] ">50K" "<=50K"levels(income_data$education.num)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15"
## [16] "16"Correlation
Now, let’s check for correlation. Are any of our variables correlated?
cor_income_data1 <- income_data %>%
select(where(is.numeric)) %>% # selecting numeric columns
cor(use = "pairwise.complete.obs") %>% # handling missing data
corrplot(type = "lower", diag = FALSE) # printing lower half of matrix
From the matrix, it appears that most variables show very light-colored circles, which suggests that there are no strong correlations between the variables.
When we’re working with a categorical outcome, we are often interested in whether the levels of our outcome are balanced, meaning whether the count of observations at one level is approximately equal to the count at the other level(s).
income_data %>%
ggplot(aes(income)) +
geom_bar()
Since >50k has 7,841 observations and <=50k has 24,720 observations, this means we will later upsample our outcome variable.
Most of my variables are categorical variables, but I plotted by one numeric variable using a histogram to see the distribution of age.
ggplot(income_data, aes(x = age)) +
geom_histogram(bins = 30, fill = "blue", color = "black") +
labs(title = "Histogram of Age", x = "Age", y = "Count") +
theme_minimal()
The histogram shows a “skewed left” distribution, where a larger number
of individuals are younger, and fewer individuals are older, with the
frequency decreasing as age increases.
I made bar charts for some categorical variables to understand the frequency of each category and how they’re related to the predictor variables.
Graph 1 shows the income distribution within various workclasses,
highlighting significant income inequality between sectors such as the
Federal government and private work. Graph 2 illustrates gender
disparities in income, with a visible difference in income levels
between males and females across the two main income categories. Graph 3
visualizes a clear trend of higher education levels correlating with
higher income, indicating the strong impact of educational attainment on
earning potential. Graph 4 suggests differences in income distribution
across marital statuses, potentially indicating economic advantages or
disadvantages associated with marital ties. Graph 5 displays the income
distribution across various occupations, identifying which professions
are more likely to earn above the $50K threshold. Graph 6 is an updated
graph on workclass and income which shows what percentage of each
workclass makes above or below 50k.
Now that we got a good feel for our data, let’s move on
Missing data
Let’s start with counting how much of my data is missing.
vis_miss(income_data)
This shows us 0.9% of the entire data set is missing. In workclass, 6% of the observations are missing, 6% in occupation, and 2% in native.country.
The 6% in workclass and occupation is significant, but not to the point we’d need to complicated of methods to deal with it.
I decided to go with imputation to deal with my missing data. Although simple, this approach is particularly beneficial because it avoids the complications and potential biases associated with listwise deletion (dropping rows with any missing values) or more complex imputation methods that require detailed assumptions about the nature of the data and the mechanisms behind the missing data.
Imputation is a method used to fill in missing data with substituted values. ‘Impute_mean’ is the process of substituting missing values in a numerical dataset with the mean (average) of the non-missing values in the same dataset. ‘Impute_mode’ involves replacing missing values with the mode, which is the most frequently occurring value in the dataset. Using the mean or mode helps in maintaining the overall distribution of the dataset, especially if the missing data is a small fraction of the total data.
We do this all in one step in our following recipe.
Setting up Models
Train/Test Split
Before building our model, we split our dataset into training, testing, and validation sets to ensure accurate evaluation and to prevent overfitting. We opted for a 75/25 split between training and testing sets. This allows the model to learn from a substantial amount of data while reserving a significant portion for unbiased evaluation.
Additionally, a validation set comprising 75% of the training data is used for model tuning. This step is crucial to refine the model without using the test set, which remains untouched to provide a final measure of the model’s performance.
We set a random seed to ensure that our data split is reproducible, and we stratify the split by the income variable to maintain consistent income category proportions across the datasets. This approach ensures that our model training and testing are both rigorous and reproducible.
set.seed(2003)
income_split <- initial_split(income_data, strata = "income", prop = 0.75)
income_train <- training(income_split)
income_test <- testing(income_split)
# Setting Up Cross-Validationt
income_valid <- validation_split(income_train, prop = .75, strata = income)## Warning: `validation_split()` was deprecated in rsample 1.2.0.
## ℹ Please use `initial_validation_split()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.Class imbalance is significant because it can lead models to exhibit bias toward the majority class. In predictive modeling, if one class significantly outnumbers the other, the model might become very good at predicting the majority class but fail to accurately predict the minority class. This can result in poor model performance, especially when the predictive accuracy of the minority class is crucial.
Let’s see the balance of each class in our dataset.
income_props <- income_train %>%
group_by(income) %>%
summarise(prop = n()/(dim(income_train)[1])) %>%
ungroup()
count(income_train)## # A tibble: 1 × 1
## n
## <int>
## 1 24420count(income_test)## # A tibble: 1 × 1
## n
## <int>
## 1 8141head(income_train)## # A tibble: 6 × 15
## age workclass fnlwgt education education.num marital.status occupation
## <dbl> <chr> <dbl> <chr> <ord> <chr> <chr>
## 1 90 <NA> 77053 HS-grad 9 Widowed <NA>
## 2 82 Private 132870 HS-grad 9 Widowed Exec-manager…
## 3 66 <NA> 186061 Some-college 10 Widowed <NA>
## 4 41 Private 264663 Some-college 10 Separated Prof-special…
## 5 29 Private 77009 11th 7 Separated Sales
## 6 61 Private 29059 HS-grad 9 Divorced Sales
## # ℹ 8 more variables: relationship <chr>, race <chr>, sex <chr>,
## # capital.gain <dbl>, capital.loss <dbl>, hours.per.week <dbl>,
## # native.country <chr>, income <fct>
Our data is imbalanced. It looks like 75% of the data is <=50k and
25% is >50k. To address this imbalance, I use upsampling within my
data preprocessing recipe. Upsampling involves artificially inflating
the minority class (>50K in this case) by replicating its
observations so that the class proportions are more balanced.
Recipe Building
To ensure our models are robust and well-prepared, we’ve developed a
standardized preprocessing recipe for our income_data.
After setting up these preprocessing steps, we check for any remaining
missing values to ensure data integrity before moving forward with
modeling.
income_recipe <- recipe(income ~ ., data = income_train) %>%
step_impute_mean(all_numeric_predictors(), -all_outcomes()) %>% # Impute missing values for numeric predictors
step_impute_mode(all_nominal_predictors(), -all_outcomes()) %>% # Impute missing values for categorical predictors
step_dummy(all_nominal_predictors()) %>%
step_center(all_predictors()) %>%
step_scale(all_predictors()) %>%
step_upsample(income, over_ratio = 0.5, skip = TRUE) The income_recipe is prepared with prep() and applied to income_train using bake(). This processes the entire dataset based on the steps defined in the recipe, such as imputation, normalization, and encoding. The processed data is then grouped by the income category, and the number of records in each income group is counted using summarise(count = n()). This helps assess the distribution of data across different income classes post-processing. Finally, the script checks for any remaining missing values in the processed data with sum(is.na(processed_data)) to ensure data integrity and completeness after preprocessing.
processed_data <- prep(income_recipe) %>% bake(new_data = income_train) %>%
group_by(income) %>%
summarise(count = n())
# Prepare the recipe with the training data
prepared_recipe <- prep(income_recipe, training = income_train)
sum(is.na(processed_data))## [1] 0The models
The general workflow of model building is setting up the model by specifying what type it is, set it to ‘classification’ for our case, and then set up the workflow.
For most of my models I then set up a tuning grid with the parameters I want tuned, tune the model, and then select the most accurate. I will then finalize the workdflow with these tuning paramters. We will then fit that model with our workflow to our training dataset.
Specify an Engine, then create a workflow
- Logistic Regression Workflow
log_reg <- logistic_reg()%>%
set_engine("glm") %>%
set_mode("classification")
#workflow
income_log_wkflow <- workflow() %>%
add_model(log_reg) %>%
add_recipe(income_recipe)
control <- control_resamples(save_pred = TRUE, save_workflow = TRUE)This workflow employs a logistic regression model configured for classification with the glm engine. It is ideal for this dataset as it effectively models binary outcomes (like income levels) and provides probabilities for each class, helping interpretation of results.
- Elastic Net
en_model <- logistic_reg(mixture = tune(),
penalty = tune()) %>%
set_mode("classification") %>%
set_engine("glmnet")
en_workflow <- workflow() %>%
add_recipe(income_recipe) %>%
add_model(en_model)The Elastic Net model uses logistic regression with both lasso and ridge penalties to handle data with multicollinearity and to enhance model selection capabilities. The logistic regression framework is adapted to tune both the mixture and penalty parameters, using the glmnet engine for optimization.
- K-Nearest Neighbors Model Workflow
## K NEAREST NEIGHBOR MODEL
knn_income_model <- nearest_neighbor(neighbors = tune()) %>%
set_mode("classification") %>%
set_engine("kknn")
knn_income_wkflow <- workflow() %>%
add_model(knn_income_model) %>%
add_recipe(income_recipe)
income_fold <- vfold_cv(income_train, v = 10)We’ll use 10 folds to perform stratified cross validation.
The KNN model, set to optimize the number of neighbors, is particularly useful for the dataset as it classifies individuals based on similarity in feature space, which can capture non-linear relationships that might exist between demographic variables and income.
- Decision Tree Model & Workflow
tree_spec <- decision_tree(cost_complexity = tune()) %>%
set_engine("rpart") %>%
set_mode("classification")
tree_wf <- workflow() %>%
add_model(tree_spec) %>%
add_recipe(income_recipe)A Decision Tree model, with tunable complexity, is incorporated for its ability to handle categorical data effectively and provide interpretable results, making it suitable for the dataset where interpretability regarding how different features impact income levels is crucial.
- Random Forest Model
rf_spec <- rand_forest(mtry = tune(),
trees = tune(),
min_n = tune()) %>%
set_engine("ranger", importance = "impurity") %>%
set_mode("classification")
rf_wf <- workflow() %>%
add_model(rf_spec) %>%
add_recipe(income_recipe)This Random Forest workflow utilizes multiple decision trees to enhance predictive accuracy and control overfitting, ideal for complex datasets with numerous predictors and relationships like the census income dataset. Its ensemble approach is robust against overfitting and capable of handling high-dimensional data effectively.
Model Tuning
#Create grids for KNN parameter (neighbors)
neighbors_grid <- grid_regular(neighbors(range = c(1, 10)), levels = 10)
neighbors_grid## # A tibble: 10 × 1
## neighbors
## <int>
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10This line creates a grid of hyperparameters for the K-Nearest Neighbors model. The grid specifies the number of neighbors to explore, ranging from 1 to 10. The grid_regular function is used to generate a regular sequence of values within this range, distributing 10 levels evenly across the range.
#tuning grind for elastic net regularization
en_grid <- grid_regular(penalty(range = c(0, 1),
trans = identity_trans()),
mixture(range = c(0, 1)),
levels = 10)For the Elastic Net model, this grid defines two parameters: penalty for the regularization strength and mixture for the balance between L1 and L2 regularization. Both parameters are varied from 0 to 1, and the identity_trans function indicates that these parameters are used as-is without transformation.
# Tuning grid for decision tree
param_grid <- grid_regular(
cost_complexity(c(-3, 1)),
levels =10
)This sets up a grid for tuning the cost complexity parameter of a decision tree, ranging from -3 to 1. This parameter controls the trade-off between tree complexity and fitting accuracy, helping to prevent overfitting.
#Tuning for Random Forests
rf_grid <- grid_regular(mtry(range = c(1, 6)),
trees(range = c(200, 600)),
min_n(range = c(10, 20)),
levels = 5)
rf_grid## # A tibble: 125 × 3
## mtry trees min_n
## <int> <int> <int>
## 1 1 200 10
## 2 2 200 10
## 3 3 200 10
## 4 4 200 10
## 5 6 200 10
## 6 1 300 10
## 7 2 300 10
## 8 3 300 10
## 9 4 300 10
## 10 6 300 10
## # ℹ 115 more rowsFor the Random Forest model, this grid tunes three parameters: mtry (the number of variables to consider for each split), trees (the number of trees in the forest), and min_n (the minimum number of samples in nodes). The ranges provide a variety of combinations to explore optimal settings for model complexity and performance.
The Tuning
“Tuning” refers to the process of finding the optimal hyperparameters for a model. This involves adjusting these parameters to find the combination that offers the best performance, typically evaluated on a validation set or through cross-validation. The objective is to optimize the model’s ability to generalize well to new data, not just to fit the existing data.
The tuning process is a systematic search through a predefined set of model configurations (the grid) to find the one that performs best according to the ROC_AUC. This is crucial for ensuring that the model not only fits the training data well but also generalizes effectively to new, unseen data, thus providing reliable predictions when deployed.
# Knn Fitting
tune_res <- tune_grid(
object = knn_income_wkflow,
resamples = income_fold,
grid = neighbors_grid,
control = control_grid(verbose = TRUE)
)
#Logistic Fitting
fit_res_log <- fit_resamples(
object = income_log_wkflow,
resamples = income_fold,
control = control_resamples(save_pred = TRUE, verbose = TRUE)
)
#Elastic Net Fitting
tune_res_en <- tune_grid(
en_workflow,
resamples = income_fold,
grid = en_grid
)
# Decision tree fitting
tune_res_DT <- tune_grid(
tree_wf,
resamples = income_fold,
grid = param_grid
)
#Random Forest Fitting
tune_rf <- tune_grid(
rf_wf,
resamples = income_fold,
grid = rf_grid
)
# Saving
save(tune_res, file = "knn.rda")
save(fit_res_log, file = "log.rda")
save(tune_res_DT, file = "DT.rda")
save(tune_rf, file = "RF.rda")
save(tune_res_en, file = "EN.rda")Model Results
One of the most useful tools for visualizing the results of models that have been tuned is the autoplot function in r. This will visualize the effects that the change in certain parameters has on our metric of choice, roc_auc.
I’m now going to load these tuning results and analyse some performance metrics. Then, I will visualize the tuning results, which helps in understanding how different hyperparameters impact model performance.
load("knn.rda")
load("log.rda")
load("EN.rda")
load("DT.rda")
load("RF.rda")
collect_metrics(tune_res)## # A tibble: 30 × 7
## neighbors .metric .estimator mean n std_err .config
## <int> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 1 accuracy binary 0.782 10 0.00292 Preprocessor1_Model01
## 2 1 brier_class binary 0.218 10 0.00292 Preprocessor1_Model01
## 3 1 roc_auc binary 0.717 10 0.00261 Preprocessor1_Model01
## 4 2 accuracy binary 0.782 10 0.00292 Preprocessor1_Model02
## 5 2 brier_class binary 0.205 10 0.00262 Preprocessor1_Model02
## 6 2 roc_auc binary 0.752 10 0.00271 Preprocessor1_Model02
## 7 3 accuracy binary 0.782 10 0.00292 Preprocessor1_Model03
## 8 3 brier_class binary 0.193 10 0.00236 Preprocessor1_Model03
## 9 3 roc_auc binary 0.779 10 0.00282 Preprocessor1_Model03
## 10 4 accuracy binary 0.782 10 0.00292 Preprocessor1_Model04
## # ℹ 20 more rowscollect_metrics(fit_res_log)## # A tibble: 3 × 6
## .metric .estimator mean n std_err .config
## <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 accuracy binary 0.842 10 0.00232 Preprocessor1_Model1
## 2 brier_class binary 0.108 10 0.00102 Preprocessor1_Model1
## 3 roc_auc binary 0.903 10 0.00194 Preprocessor1_Model1collect_metrics(tune_res_en)## # A tibble: 300 × 8
## penalty mixture .metric .estimator mean n std_err .config
## <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0 0 accuracy binary 0.840 10 0.00179 Preprocessor1_Mo…
## 2 0 0 brier_class binary 0.111 10 0.000905 Preprocessor1_Mo…
## 3 0 0 roc_auc binary 0.899 10 0.00150 Preprocessor1_Mo…
## 4 0.111 0 accuracy binary 0.836 10 0.00229 Preprocessor1_Mo…
## 5 0.111 0 brier_class binary 0.115 10 0.000864 Preprocessor1_Mo…
## 6 0.111 0 roc_auc binary 0.893 10 0.00161 Preprocessor1_Mo…
## 7 0.222 0 accuracy binary 0.834 10 0.00266 Preprocessor1_Mo…
## 8 0.222 0 brier_class binary 0.120 10 0.000831 Preprocessor1_Mo…
## 9 0.222 0 roc_auc binary 0.890 10 0.00173 Preprocessor1_Mo…
## 10 0.333 0 accuracy binary 0.833 10 0.00243 Preprocessor1_Mo…
## # ℹ 290 more rowscollect_metrics(tune_res_DT)## # A tibble: 30 × 7
## cost_complexity .metric .estimator mean n std_err .config
## <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0.001 accuracy binary 0.847 10 0.00302 Preprocessor1_Mod…
## 2 0.001 brier_class binary 0.108 10 0.00157 Preprocessor1_Mod…
## 3 0.001 roc_auc binary 0.890 10 0.00323 Preprocessor1_Mod…
## 4 0.00278 accuracy binary 0.838 10 0.00456 Preprocessor1_Mod…
## 5 0.00278 brier_class binary 0.114 10 0.00181 Preprocessor1_Mod…
## 6 0.00278 roc_auc binary 0.869 10 0.00340 Preprocessor1_Mod…
## 7 0.00774 accuracy binary 0.835 10 0.00402 Preprocessor1_Mod…
## 8 0.00774 brier_class binary 0.116 10 0.00163 Preprocessor1_Mod…
## 9 0.00774 roc_auc binary 0.862 10 0.00344 Preprocessor1_Mod…
## 10 0.0215 accuracy binary 0.841 10 0.00244 Preprocessor1_Mod…
## # ℹ 20 more rowscollect_metrics(tune_rf)## # A tibble: 375 × 9
## mtry trees min_n .metric .estimator mean n std_err .config
## <int> <int> <int> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 1 200 10 accuracy binary 0.795 10 0.00252 Preprocessor1_…
## 2 1 200 10 brier_class binary 0.150 10 0.000716 Preprocessor1_…
## 3 1 200 10 roc_auc binary 0.883 10 0.00243 Preprocessor1_…
## 4 2 200 10 accuracy binary 0.838 10 0.00267 Preprocessor1_…
## 5 2 200 10 brier_class binary 0.125 10 0.000699 Preprocessor1_…
## 6 2 200 10 roc_auc binary 0.893 10 0.00178 Preprocessor1_…
## 7 3 200 10 accuracy binary 0.847 10 0.00234 Preprocessor1_…
## 8 3 200 10 brier_class binary 0.115 10 0.000828 Preprocessor1_…
## 9 3 200 10 roc_auc binary 0.900 10 0.00182 Preprocessor1_…
## 10 4 200 10 accuracy binary 0.849 10 0.00238 Preprocessor1_…
## # ℹ 365 more rowsVisualizing Results
## Warning in ggplot2::scale_x_continuous(trans = trans): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.
## # A tibble: 24,420 × 7
## .pred_class `.pred_>50K` `.pred_<=50K` id .row income .config
## <fct> <dbl> <dbl> <chr> <int> <fct> <chr>
## 1 >50K 0.691 0.309 Fold01 2 <=50K Preprocessor1_Mod…
## 2 >50K 0.713 0.287 Fold01 8 <=50K Preprocessor1_Mod…
## 3 >50K 0.846 0.154 Fold01 16 <=50K Preprocessor1_Mod…
## 4 <=50K 0.161 0.839 Fold01 26 <=50K Preprocessor1_Mod…
## 5 >50K 0.599 0.401 Fold01 28 <=50K Preprocessor1_Mod…
## 6 >50K 0.604 0.396 Fold01 73 <=50K Preprocessor1_Mod…
## 7 >50K 0.710 0.290 Fold01 82 <=50K Preprocessor1_Mod…
## 8 <=50K 0.469 0.531 Fold01 92 <=50K Preprocessor1_Mod…
## 9 >50K 0.571 0.429 Fold01 93 <=50K Preprocessor1_Mod…
## 10 >50K 0.809 0.191 Fold01 112 <=50K Preprocessor1_Mod…
## # ℹ 24,410 more rows## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.903
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.903
The show_best() function is part of the tidymodels suite and is used specifically within the tuning results context to extract the top-performing models according to a specified metric, which in this case is the ROC AUC.
Choosing the Best
# Getting best KNN Model
best_knn_income <- show_best(tune_res, metric = "roc_auc")
best_knn_income## # A tibble: 5 × 7
## neighbors .metric .estimator mean n std_err .config
## <int> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 10 roc_auc binary 0.835 10 0.00255 Preprocessor1_Model10
## 2 9 roc_auc binary 0.831 10 0.00256 Preprocessor1_Model09
## 3 8 roc_auc binary 0.828 10 0.00232 Preprocessor1_Model08
## 4 7 roc_auc binary 0.823 10 0.00272 Preprocessor1_Model07
## 5 6 roc_auc binary 0.817 10 0.00245 Preprocessor1_Model06# Getting best log Model
best_lm_income <- show_best(fit_res_log, metric = "roc_auc")
best_lm_income## # A tibble: 1 × 6
## .metric .estimator mean n std_err .config
## <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 roc_auc binary 0.903 10 0.00194 Preprocessor1_Model1# Getting best EN Model
best_en_income <- show_best(tune_res_en, metric = "roc_auc")
best_en_income## # A tibble: 5 × 8
## penalty mixture .metric .estimator mean n std_err .config
## <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0 0.889 roc_auc binary 0.905 10 0.00154 Preprocessor1_Model081
## 2 0 0.667 roc_auc binary 0.905 10 0.00154 Preprocessor1_Model061
## 3 0 0.778 roc_auc binary 0.905 10 0.00154 Preprocessor1_Model071
## 4 0 1 roc_auc binary 0.905 10 0.00155 Preprocessor1_Model091
## 5 0 0.111 roc_auc binary 0.905 10 0.00153 Preprocessor1_Model011# Getting best DT Model
best_dt_income <- show_best(tune_res_DT, metric = "roc_auc")
best_dt_income## # A tibble: 5 × 7
## cost_complexity .metric .estimator mean n std_err .config
## <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0.001 roc_auc binary 0.890 10 0.00323 Preprocessor1_Model01
## 2 0.00278 roc_auc binary 0.869 10 0.00340 Preprocessor1_Model02
## 3 0.00774 roc_auc binary 0.862 10 0.00344 Preprocessor1_Model03
## 4 0.0215 roc_auc binary 0.844 10 0.00306 Preprocessor1_Model04
## 5 0.0599 roc_auc binary 0.817 10 0.00508 Preprocessor1_Model05# Getting best RF Model
best_rf_income <- show_best(tune_rf, metric = "roc_auc")
best_rf_income## # A tibble: 5 × 9
## mtry trees min_n .metric .estimator mean n std_err .config
## <int> <int> <int> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 6 200 20 roc_auc binary 0.911 10 0.00164 Preprocessor1_Model1…
## 2 6 400 20 roc_auc binary 0.910 10 0.00152 Preprocessor1_Model1…
## 3 6 600 20 roc_auc binary 0.910 10 0.00161 Preprocessor1_Model1…
## 4 6 500 20 roc_auc binary 0.910 10 0.00157 Preprocessor1_Model1…
## 5 6 300 12 roc_auc binary 0.910 10 0.00157 Preprocessor1_Model0…best_rf_model <- select_best(tune_rf, metric = "roc_auc")The Random Forest model#105 has the highest recorded ROC AUC value at approximately 0.9105, suggesting it is the best performer among the models evaluated based on my data for discriminating between the binary classes effectively.
The best-performing configuration in terms of ROC AUC is when the Random Forest is configured with 200 trees and 20 as the minimum node size. This configuration not only achieved the highest mean ROC AUC but also maintains a low standard error, enhancing confidence in the model’s stability and reliability.
Final Model
final_rf_model <- finalize_workflow(rf_wf,
best_rf_model)
final_rf_model <- fit(final_rf_model,
data = income_test)final_rf_model %>% extract_fit_parsnip() %>%
vip() +
theme_minimal()
The three most important predictors was marital status with a civilian
spouse, age, and never being married.
Let’s look at the testing ROC AUC
Final Model Results
# Ensure that predictions are generated correctly
predictions_rf <- predict(final_rf_model, income_test, type = "prob")
# Combine predictions with the true labels
predictions_rf <- bind_cols(income_test %>% select(income), predictions_rf)
# Check the structure to verify the combined dataframe
str(predictions_rf)## tibble [8,141 × 3] (S3: tbl_df/tbl/data.frame)
## $ income : Factor w/ 2 levels ">50K","<=50K": 2 2 2 2 1 2 1 1 1 1 ...
## $ .pred_>50K : num [1:8141] 0.123 0.19 0.153 0.214 0.686 ...
## $ .pred_<=50K: num [1:8141] 0.877 0.81 0.847 0.786 0.314 ...# Evaluate model performance using ROC AUC
roc_results <- roc_auc(predictions_rf, truth = income, `.pred_>50K`)
# Print ROC AUC results
print(roc_results)## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.937roc_data <- roc_curve(predictions_rf, truth = income, `.pred_>50K`)
autoplot(roc_data) +
labs(title = "ROC Curve", x = "1 - Specificity", y = "Sensitivity") +
theme_minimal()
Conclusion
An ROC AUC of 0.9363215 is quite high, indicating that my model has a strong ability to distinguish between individuals earning more than $50,000 and those earning less than or equal to $50,000. In other words, there is a 93.63% chance that my model will correctly differentiate between a randomly chosen positive instance (income > 50K) and a randomly chosen negative instance (income <= 50K).
n this project, we aimed to predict whether an individual earns more than $50,000 per year using the Adult Census Income dataset. We explored a variety of machine learning models, including logistic regression, elastic net, k-nearest neighbors, decision trees, and random forests. The data was carefully preprocessed to handle missing values, normalize numerical features, and encode categorical variables. Class imbalance was addressed through upsampling.
We split the dataset into training, testing, and validation sets to ensure robust model evaluation. Hyperparameter tuning was conducted for each model using cross-validation and the ROC AUC metric to select the best performing configurations.
After extensive tuning and evaluation, the Random Forest model emerged as the best performer with a ROC AUC of approximately 0.9363. This high ROC AUC value indicates a strong ability of the model to distinguish between individuals earning more than $50,000 and those earning less. The most important predictors identified were marital status with a civilian spouse, age, and never being married.
The final Random Forest model was validated on the test set, achieving a ROC AUC of 0.9363, confirming its effectiveness and robustness. This high performance suggests that the model is well-suited for predicting income levels based on demographic and socioeconomic features from the census data.
In summary, this project successfully demonstrates the application of machine learning techniques to predict income levels, highlighting the importance of careful data preprocessing, model selection, and hyperparameter tuning in achieving high predictive accuracy. The findings underscore the significant impact of marital status, age, and marital status on income, providing valuable insights for socioeconomic analyses.