Mehlika ALGAN - 20023011   Burak TURAN - 22023082   Atalay AKTAŞ - 22023619

1. Introduction

Employee attrition — the voluntary or involuntary departure of staff from an organisation — is one of the most persistent and financially consequential challenges in human resource management. Studies consistently estimate that replacing a single employee costs between 50% and 200% of their annual salary, when accounting for recruitment, onboarding, and productivity losses during the transition period. Beyond direct costs, high attrition erodes institutional knowledge, disrupts team dynamics, and reduces overall organisational performance.

Understanding the structural drivers of attrition is therefore a high-value analytical problem. Rather than reacting to departures after the fact, organisations that can identify who is likely to leave and why gain the capacity for targeted, pre-emptive retention interventions — whether through compensation adjustments, career development programmes, or flexible work arrangements.

This report presents a complete, end-to-end analysis of the IBM HR Analytics Employee Attrition dataset, which contains records for 1,470 employees across 35 variables. The investigation is structured around three objectives:

  1. Exploratory Data Analysis: Understand the distribution of key HR variables and identify patterns associated with attrition.
  2. Hypothesis Testing: Apply formal statistical tests to evaluate five specific hypotheses about the drivers of attrition, using both parametric and non-parametric methods appropriate to the data.
  3. Predictive Modelling: Build a Random Forest classification model within a tidymodels cross-validation framework to predict attrition risk.

2. Libraries

library(tidyverse)
library(ggplot2)
library(rstatix)
library(ggpubr)
library(tidymodels)
library(ranger)
library(plotly)
library(DT)
library(knitr)
library(kableExtra)
library(scales)

3. Data Loading & Preprocessing

3.1 Load and Inspect

We begin by loading the dataset and performing a structural audit — examining dimensions, data types, and the presence of missing values before committing to any modelling decisions.

data <- read.csv(file.choose())
cat("Rows:", nrow(data), "\n")
## Rows: 1470
cat("Columns:", ncol(data), "\n")
## Columns: 35
cat("Missing values:", sum(is.na(data)), "\n")
## Missing values: 0

The dataset contains 1,470 employee records across 35 variables with no missing values — a notable quality advantage that eliminates the need for any imputation strategy.

3.2 Remove Irrelevant Columns

Several columns in the raw dataset carry no analytical signal and should be removed before modelling:

  • EmployeeCount and StandardHours are constants — every record has the same value, providing zero discriminative power.
  • Over18 is a single-level categorical variable — all employees are over 18.
  • PerformanceRating, StockOptionLevel, WorkLifeBalance, DailyRate, HourlyRate, and MonthlyRate are either redundant given other compensation variables or have documented encoding ambiguities in this dataset.

Retaining zero-variance or near-zero-variance predictors can introduce noise into models, particularly in penalised regression methods; dropping them is standard pre-modelling hygiene.

data <- data[, !(names(data) %in% c(
  "EmployeeCount", "StandardHours", "Over18",
  "PerformanceRating", "StockOptionLevel", "WorkLifeBalance",
  "DailyRate", "HourlyRate", "MonthlyRate"
))]

cat("Remaining columns:", ncol(data))
## Remaining columns: 26

3.3 Factor Encoding

Machine learning algorithms and statistical tests in R require categorical data to be explicitly typed as factor. All character columns are batch-converted. Additionally, several numeric columns represent Likert-scale responses — treating these as raw integers would incorrectly imply equal spacing between levels. We encode them as named factors with meaningful labels.

data <- data %>%
  mutate(across(where(is.character), as.factor))

data$Education <- factor(data$Education,
  levels = 1:5,
  labels = c("Below College", "College", "Bachelor", "Master", "Doctor"))

data$EnvironmentSatisfaction <- factor(data$EnvironmentSatisfaction,
  levels = 1:4, labels = c("Low", "Medium", "High", "Very High"))

data$JobInvolvement <- factor(data$JobInvolvement,
  levels = 1:4, labels = c("Low", "Medium", "High", "Very High"))

data$JobSatisfaction <- factor(data$JobSatisfaction,
  levels = 1:4, labels = c("Low", "Medium", "High", "Very High"))

data$RelationshipSatisfaction <- factor(data$RelationshipSatisfaction,
  levels = 1:4, labels = c("Low", "Medium", "High", "Very High"))

data$JobLevel <- as.factor(data$JobLevel)

cat("Factor encoding complete.")
## Factor encoding complete.

4. Exploratory Data Analysis

4.1 Attrition Class Distribution

Before any modelling, understanding the class imbalance in the target variable is the single most important diagnostic step. In employee attrition problems, the departing group (the minority class) is almost always considerably smaller than the retained group. Failing to account for this leads models to simply predict “No Attrition” for nearly all records and still report misleadingly high accuracy — a phenomenon known as the accuracy paradox.

attr_counts <- data %>%
  count(Attrition) %>%
  mutate(pct   = n / sum(n),
         label = paste0(n, "\n(", percent(pct, 0.1), ")"))

p1 <- ggplot(attr_counts, aes(x = Attrition, y = n, fill = Attrition,
             text = paste0("Group: ", Attrition, "<br>Count: ", n,
                           "<br>Share: ", percent(pct, 0.1)))) +
  geom_col(width = 0.45, show.legend = FALSE) +
  geom_text(aes(label = label), vjust = -0.4, fontface = "bold", size = 3.8) +
  scale_fill_manual(values = c("No" = "#2ecc71", "Yes" = "#e74c3c")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(title = "Attrition Class Distribution", x = "Attrition", y = "Count") +
  theme_minimal(base_size = 12)

ggplotly(p1, tooltip = "text") %>% layout(showlegend = FALSE) %>%
  config(displayModeBar = FALSE)

Approximately 16.1% of employees have left the organisation — a ~5:1 class imbalance typical of enterprise HR datasets. This has a direct implication for model evaluation: metrics such as recall, F1-score, and ROC-AUC must be prioritised over raw accuracy.

4.2 Monthly Income Distribution

Monthly income is central to this analysis — it appears both as the outcome variable in our regression model and as a key predictor of attrition. We examine its univariate distribution to identify skewness, modality, and the relationship between mean and median.

med_inc  <- median(data$MonthlyIncome)
mean_inc <- mean(data$MonthlyIncome)

p2 <- ggplot(data, aes(x = MonthlyIncome)) +
  geom_histogram(fill = "#3498db", color = "white", bins = 30, alpha = 0.85) +
  geom_vline(xintercept = med_inc,  color = "#e74c3c", linetype = "dashed", linewidth = 0.9) +
  geom_vline(xintercept = mean_inc, color = "#2c3e50", linetype = "dashed", linewidth = 0.9) +
  annotate("text", x = med_inc  + 280, y = 190,
           label = paste0("Median: $", format(med_inc,  big.mark = ",")),
           color = "#e74c3c", fontface = "bold", hjust = 0, size = 3.4) +
  annotate("text", x = mean_inc + 280, y = 160,
           label = paste0("Mean: $",   format(round(mean_inc), big.mark = ",")),
           color = "#2c3e50", fontface = "bold", hjust = 0, size = 3.4) +
  labs(title = "Monthly Income Distribution",
       x = "Monthly Income (USD)", y = "Frequency") +
  theme_minimal(base_size = 12)

ggplotly(p2) %>% config(displayModeBar = FALSE)

The distribution is right-skewed — the mean (~$6,503) exceeds the median (~$4,919), a hallmark of salary distributions where a small number of senior executives pull the average upward. This skewness has two practical implications: (1) parametric tests assuming normality will be inappropriate for income-related hypotheses, and (2) the income variable may benefit from log-transformation in regression contexts.

4.3 Attrition by Department and Job Role

Attrition is rarely uniform across an organisation. Examining departure rates by department and job role reveals which parts of the workforce are structurally more vulnerable — and helps focus retention efforts.

dept_rate <- data %>%
  count(Department, Attrition) %>%
  group_by(Department) %>%
  mutate(rate = n / sum(n)) %>%
  filter(Attrition == "Yes")

p3 <- ggplot(dept_rate, aes(x = reorder(Department, rate), y = rate, fill = Department,
             text = paste0("Department: ", Department,
                           "<br>Attrition Rate: ", percent(rate, 0.1),
                           "<br>Leavers: ", n))) +
  geom_col(width = 0.5, show.legend = FALSE) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_brewer(palette = "Set2") +
  coord_flip() +
  labs(title = "Attrition Rate by Department", x = NULL, y = "Attrition Rate") +
  theme_minimal(base_size = 12)

ggplotly(p3, tooltip = "text") %>% config(displayModeBar = FALSE)
role_rate <- data %>%
  count(JobRole, Attrition) %>%
  group_by(JobRole) %>%
  mutate(rate = n / sum(n)) %>%
  filter(Attrition == "Yes")

p4 <- ggplot(role_rate, aes(x = reorder(JobRole, rate), y = rate, fill = rate,
             text = paste0("Role: ", JobRole,
                           "<br>Attrition Rate: ", percent(rate, 0.1),
                           "<br>Leavers: ", n))) +
  geom_col(width = 0.6, show.legend = FALSE) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_gradient(low = "#aed6f1", high = "#e74c3c") +
  coord_flip() +
  labs(title = "Attrition Rate by Job Role", x = NULL, y = "Attrition Rate") +
  theme_minimal(base_size = 12)

ggplotly(p4, tooltip = "text") %>% config(displayModeBar = FALSE)

The Sales department records the highest departmental exit rate, a pattern commonly attributed to commission-based pressure and a competitive external labour market for sales talent. At the role level, Sales Representatives and Laboratory Technicians stand out — both are entry-to-mid level roles with relatively lower compensation ceilings, creating a natural pipeline of employees who leave to seek advancement elsewhere.

4.4 Monthly Income by Attrition Status

p5 <- ggplot(data, aes(x = Attrition, y = MonthlyIncome, fill = Attrition,
             text = paste0("Attrition: ", Attrition,
                           "<br>Income: $", format(MonthlyIncome, big.mark = ",")))) +
  geom_boxplot(width = 0.4, outlier.size = 1.2, outlier.alpha = 0.5) +
  scale_fill_manual(values = c("No" = "#2ecc71", "Yes" = "#e74c3c")) +
  scale_y_continuous(labels = dollar_format()) +
  labs(title = "Monthly Income by Attrition Status",
       x = "Attrition", y = "Monthly Income (USD)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

ggplotly(p5, tooltip = "text") %>% config(displayModeBar = FALSE)

Employees who left are concentrated at lower income levels — the median and interquartile range of leavers both sit well below those of retained employees. This visual evidence directly motivates Hypothesis 1.

5. Hypothesis Testing

All hypothesis tests follow a systematic four-step protocol: (1) state the null and alternative hypotheses, (2) assess distributional assumptions via the Shapiro-Wilk normality test, (3) select the appropriate test (parametric or non-parametric), and (4) interpret the result in both statistical and business terms.

Hypothesis 1 — Income and Attrition

Research Question: Do employees who leave the organisation earn significantly less than the population average?

H₀: The median monthly income of leavers equals the overall mean ($6,500).
H₁: The median monthly income of leavers is less than $6,500 (one-tailed).

Step 1 — Normality Check (Shapiro-Wilk):

h1_data <- data %>% filter(Attrition == "Yes") %>% select(MonthlyIncome)

shapiro_test(h1_data$MonthlyIncome) %>%
  kbl(col.names = c("Variable", "W Statistic", "p-value"),
      caption = "Shapiro-Wilk Normality Test — Leavers' Monthly Income") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Shapiro-Wilk Normality Test — Leavers’ Monthly Income
Variable W Statistic p-value
h1_data$MonthlyIncome 0.7798969 0

p < 0.05 → normality is rejected. We must use a non-parametric alternative: the Wilcoxon Signed-Rank Test.

Step 2 — Wilcoxon Test vs. Overall Mean ($6,500):

h1_data %>%
  wilcox_test(MonthlyIncome ~ 1, mu = 6500, alternative = "less") %>%
  kbl(caption = "Wilcoxon Signed-Rank Test — Leavers vs. Overall Mean ($6,500)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Wilcoxon Signed-Rank Test — Leavers vs. Overall Mean ($6,500)
.y. group1 group2 n statistic p
MonthlyIncome 1 null model 237 5813.5 0

Step 3 — Wilcoxon Test vs. Overall Median ($4,900):

We extend the test to the median (~$4,900) to determine whether leavers earn less even relative to a more robust central tendency measure — strengthening the conclusion if both tests agree.

h1_data %>%
  wilcox_test(MonthlyIncome ~ 1, mu = 4900, alternative = "less") %>%
  kbl(caption = "Wilcoxon Signed-Rank Test — Leavers vs. Overall Median ($4,900)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Wilcoxon Signed-Rank Test — Leavers vs. Overall Median ($4,900)
.y. group1 group2 n statistic p
MonthlyIncome 1 null model 237 10868 0.00111

Conclusion: Both tests yield p < 0.05. We reject H₀ in both cases. Leavers earn significantly less than both the mean and the median income of the full workforce. From a business standpoint, this strongly suggests that below-market compensation is a meaningful driver of voluntary attrition. Structured pay equity reviews and market-rate benchmarking may help retain at-risk employees.

Hypothesis 2 — Marital Status and Attrition

Research Question: Is there a statistically significant association between an employee’s marital status and their propensity to leave?

H₀: Marital status and attrition are independent (no association).
H₁: They are not independent — marital status is associated with attrition.

The Chi-Squared test of independence is appropriate here: both variables are categorical and we are testing for association rather than direction.

h2_data <- data %>% select(Attrition, MaritalStatus)

p6 <- ggplot(h2_data, aes(x = MaritalStatus, fill = Attrition,
             text = paste0("Status: ", MaritalStatus, "<br>Attrition: ", Attrition))) +
  geom_bar(position = "dodge") +
  scale_fill_manual(values = c("No" = "#2ecc71", "Yes" = "#e74c3c")) +
  labs(title = "Attrition Count by Marital Status",
       x = "Marital Status", y = "Count", fill = "Attrition") +
  theme_minimal(base_size = 12)

ggplotly(p6, tooltip = "text") %>% config(displayModeBar = FALSE)
h2_table <- table(h2_data$MaritalStatus, h2_data$Attrition)
h2_res   <- chisq.test(h2_table)

data.frame(
  Statistic   = round(h2_res$statistic, 4),
  df          = h2_res$parameter,
  p.value     = signif(h2_res$p.value, 4),
  Significant = ifelse(h2_res$p.value < 0.001, "***",
                ifelse(h2_res$p.value < 0.01,  "**",
                ifelse(h2_res$p.value < 0.05,  "*", "ns")))
) %>%
  kbl(row.names = FALSE,
      caption = "Chi-Squared Test of Independence — Marital Status × Attrition") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Chi-Squared Test of Independence — Marital Status × Attrition
Statistic df p.value Significant
46.1637 2 0 ***

Conclusion: p < 0.05 → We reject H₀. Marital status and attrition are not independent. Single employees exhibit a substantially higher attrition rate compared to married or divorced colleagues. This is consistent with behavioural economics research: single employees face fewer relocation costs and family-tied constraints, making job switching comparatively less costly. Organisations may wish to tailor retention programmes — flexible benefits, mentorship tracks, social inclusion initiatives — to this demographic segment.

Hypothesis 3 — Distance from Home and Attrition

Research Question: Do employees who leave live significantly further from the workplace than those who stay?

H₀: The distribution of commuting distance is identical between leavers and stayers.
H₁: Leavers live further from the office than stayers.

Summary statistics by group:

data %>%
  select(Attrition, DistanceFromHome) %>%
  group_by(Attrition) %>%
  get_summary_stats(type = "mean_sd") %>%
  kbl(caption = "Distance from Home — Summary Statistics by Attrition Status") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Distance from Home — Summary Statistics by Attrition Status
Attrition variable n mean sd
No DistanceFromHome 1233 8.916 8.013
Yes DistanceFromHome 237 10.633 8.453

Normality check:

data %>%
  select(Attrition, DistanceFromHome) %>%
  group_by(Attrition) %>%
  shapiro_test(DistanceFromHome) %>%
  kbl(caption = "Shapiro-Wilk Test — Distance from Home by Group") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Shapiro-Wilk Test — Distance from Home by Group
Attrition variable statistic p
No DistanceFromHome 0.8534768 0
Yes DistanceFromHome 0.8959438 0

Normality is violated in both groups → Mann-Whitney U Test (non-parametric two-sample location test):

data %>%
  select(Attrition, DistanceFromHome) %>%
  wilcox_test(DistanceFromHome ~ Attrition, paired = FALSE) %>%
  kbl(caption = "Mann-Whitney U Test — Distance from Home: Leavers vs. Stayers") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Mann-Whitney U Test — Distance from Home: Leavers vs. Stayers
.y. group1 group2 n1 n2 statistic p
DistanceFromHome No Yes 1233 237 127995.5 0.00239

Conclusion: p < 0.05 → We reject H₀. The commuting distance distributions of leavers and stayers come from different populations. Leavers tend to live further from the office. Long commutes are a documented source of burnout and life dissatisfaction, even when controlling for income. Remote work provisions or transportation subsidies represent practical policy responses to this finding.

Hypothesis 4 — Career Stagnation Index by Marital Status

Research Question: Does career stagnation differ across marital status groups?

We engineer a Stagnation Index = YearsSinceLastPromotion / YearsAtCompany. A higher value means the employee has spent a larger proportion of their tenure without a promotion — capturing relative, rather than absolute, career stagnation.

H₀: Stagnation index distributions are identical across all marital status groups.
H₁: At least one group has a different stagnation index distribution (Kruskal-Wallis).

data <- data %>%
  mutate(StagnationIndex = ifelse(YearsAtCompany == 0, 0,
                                  YearsSinceLastPromotion / YearsAtCompany))
p7 <- ggplot(data, aes(x = MaritalStatus, y = StagnationIndex, fill = MaritalStatus,
             text = paste0("Status: ", MaritalStatus,
                           "<br>Stagnation Index: ", round(StagnationIndex, 2)))) +
  geom_boxplot(width = 0.4, outlier.size = 1.2, outlier.alpha = 0.5,
               show.legend = FALSE) +
  scale_fill_brewer(palette = "Pastel1") +
  labs(title = "Career Stagnation Index by Marital Status",
       x = "Marital Status", y = "Stagnation Index") +
  theme_minimal(base_size = 12)

ggplotly(p7, tooltip = "text") %>% config(displayModeBar = FALSE)
data %>%
  group_by(MaritalStatus) %>%
  shapiro_test(StagnationIndex) %>%
  kbl(caption = "Shapiro-Wilk Test — Stagnation Index by Marital Status") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Shapiro-Wilk Test — Stagnation Index by Marital Status
MaritalStatus variable statistic p
Divorced StagnationIndex 0.7814247 0
Married StagnationIndex 0.8196706 0
Single StagnationIndex 0.7817036 0 
data %>%
  kruskal_test(StagnationIndex ~ MaritalStatus) %>%
  kbl(caption = "Kruskal-Wallis Test — Stagnation Index across Marital Status Groups") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Kruskal-Wallis Test — Stagnation Index across Marital Status Groups
.y. n statistic df p method
StagnationIndex 1470 3.950941 2 0.139 Kruskal-Wallis

Sensitivity check (excluding employees with 0 years at company to avoid division artifacts):

data %>%
  filter(YearsAtCompany > 0) %>%
  mutate(StagnationIndex2 = YearsSinceLastPromotion / YearsAtCompany) %>%
  kruskal_test(StagnationIndex2 ~ MaritalStatus) %>%
  kbl(caption = "Kruskal-Wallis — Sensitivity Check (Tenure > 0 years)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Kruskal-Wallis — Sensitivity Check (Tenure > 0 years)
.y. n statistic df p method
StagnationIndex2 1426 3.656471 2 0.161 Kruskal-Wallis

Conclusion: Both tests yield p > 0.05 (≈ 0.14 and 0.16). We fail to reject H₀. There is no statistically significant difference in career stagnation patterns across marital status groups. Career stagnation appears to be driven more by department-level or managerial factors than by marital status — a direction worth exploring in future analyses.

Hypothesis 5 — Gender, Job Level, and Income Among Leavers

Research Question: Among mid-level employees (Job Level 3), do gender and attrition status jointly predict meaningful differences in monthly income?

This investigation is motivated by pay equity concerns: if female employees who leave earn significantly less than their male or retained counterparts at the same job level, it suggests a potential compensation gap contributing to gendered attrition patterns.

H₀: Monthly income is the same across all gender × attrition combinations at Job Level 3.
H₁: At least one group earns significantly more or less than the others.

p8 <- data %>%
  filter(JobLevel %in% c(1, 2, 3)) %>%
  ggplot(aes(x = Gender, y = MonthlyIncome, fill = Attrition,
             text = paste0("Gender: ", Gender,
                           "<br>Attrition: ", Attrition,
                           "<br>Income: $", format(MonthlyIncome, big.mark = ","),
                           "<br>Level: ", JobLevel))) +
  geom_boxplot(width = 0.4, outlier.size = 1, outlier.alpha = 0.5) +
  facet_wrap(~JobLevel,
             labeller = labeller(JobLevel = c("1"="Level 1","2"="Level 2","3"="Level 3"))) +
  scale_fill_manual(values = c("No" = "#2ecc71", "Yes" = "#e74c3c")) +
  scale_y_continuous(labels = dollar_format()) +
  labs(title = "Monthly Income by Gender, Attrition, and Job Level",
       x = "Gender", y = "Monthly Income (USD)", fill = "Attrition") +
  theme_minimal(base_size = 11)

ggplotly(p8, tooltip = "text") %>% config(displayModeBar = FALSE)
h5_data <- data %>%
  filter(JobLevel == 3) %>%
  mutate(Group = paste(Gender, Attrition, sep = "_"))

h5_data %>%
  group_by(Group) %>%
  shapiro_test(MonthlyIncome) %>%
  kbl(caption = "Shapiro-Wilk Test — Income by Gender × Attrition (Level 3)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Shapiro-Wilk Test — Income by Gender × Attrition (Level 3)
Group variable statistic p
Female_No MonthlyIncome 0.9539652 0.0051013
Female_Yes MonthlyIncome 0.7905755 0.0073350
Male_No MonthlyIncome 0.9722731 0.0278021
Male_Yes MonthlyIncome 0.9002564 0.0416979 
h5_data %>%
  kruskal_test(MonthlyIncome ~ Group) %>%
  kbl(caption = "Kruskal-Wallis Test — Income across Gender × Attrition Groups (Level 3)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Kruskal-Wallis Test — Income across Gender × Attrition Groups (Level 3)
.y. n statistic df p method
MonthlyIncome 218 11.03587 3 0.0115 Kruskal-Wallis 
h5_data %>%
  dunn_test(MonthlyIncome ~ Group, p.adjust.method = "bonferroni") %>%
  select(group1, group2, statistic, p, p.adj, p.adj.signif) %>%
  kbl(caption = "Dunn Post-Hoc Test — Pairwise Comparisons (Bonferroni-Adjusted)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Dunn Post-Hoc Test — Pairwise Comparisons (Bonferroni-Adjusted)
group1 group2 statistic p p.adj p.adj.signif
Female_No Female_Yes -3.1569310 0.0015944 0.0095663 **
Female_No Male_No -1.7065956 0.0878972 0.5273834 ns
Female_No Male_Yes -0.2814978 0.7783287 1.0000000 ns
Female_Yes Male_No 2.3737522 0.0176084 0.1056502 ns
Female_Yes Male_Yes 2.4799001 0.0131419 0.0788515 ns
Male_No Male_Yes 0.7447248 0.4564381 1.0000000 ns

Conclusion: Kruskal-Wallis p ≈ 0.015 — at least one group differs significantly. The Dunn post-hoc test (with Bonferroni correction for multiple comparisons) identifies the source: female employees who left at Job Level 3 earn significantly less than their male counterparts at the same level. All other pairwise comparisons are non-significant. This is a diagnostically important finding: it suggests compensation inequity at mid-management levels may be a gender-specific attrition driver, warranting targeted pay equity audits within this cohort.

6. Predictive Modelling

6.1 Linear Regression (for Monthly Income)

Definition:
Linear regression is a supervised learning method that models the relationship between a continuous dependent variable (response) and one or more independent variables (predictors) by fitting a linear equation to observed data.

Purpose:
The purpose is to predict the value of the outcome variable (MonthlyIncome) based on the values of predictors, and to quantify how strongly each predictor influences the outcome.

Logic:
The method assumes that the relationship between predictors \(\mathbf{X}\) and the outcome \(Y\) is linear:

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p + \epsilon \]

Coefficients \(\beta\) are estimated by minimising the sum of squared residuals (Ordinary Least Squares).

Formula:
For a single predictor:

\[ \hat{y} = b_0 + b_1 x \]

For multiple predictors (matrix notation):

\[ \hat{\mathbf{y}} = \mathbf{X}\mathbf{b} \quad \text{where} \quad \mathbf{b} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} \]

Suggested use cases:
- When the outcome is continuous (e.g., salary, sales, temperature)
- When a linear relationship is plausible or can be approximated
- When interpretability of coefficients is important (e.g., “holding other variables constant, a one‑unit change in \(X\) is associated with a \(\beta\) change in \(Y\)”)
- Baseline for more complex models

Reasoning for method choice:
MonthlyIncome is a continuous, right‑skewed variable but still suitable for linear regression after standardisation. The dataset contains many numeric and categorical features that can be linearly combined to explain income variation. Linear regression provides a transparent, interpretable baseline before moving to non‑linear methods. It also serves as a diagnostic tool: systematic patterns in residuals would indicate the need for transformations or tree‑based models. Given the moderate sample size (\(n = 1,470\)), linear regression is computationally efficient and its coefficients can directly inform HR policies (e.g., “an additional year of tenure corresponds to a \(\$\)X increase in income”).

6.1.1 Monthly Income Regression

Before the classification task, we build a multiple linear regression model to predict monthly income using the full tidymodels pipeline. This serves two purposes: it quantifies how well employee attributes explain compensation, and it demonstrates the preprocessing workflow that is subsequently reused for the Random Forest classifier.

Preprocessing Rationale

The preprocessing recipe applies two main steps:

  • One-hot encoding (step_dummy(one_hot = TRUE)): All categorical predictors (e.g., Department, JobRole, MaritalStatus) are nominal—there is no ordinal relationship between categories. One-hot encoding creates binary columns for each level, allowing linear and tree‑based models to use these variables without imposing a false order.

  • Standardisation (step_normalize()): Numeric predictors (e.g., MonthlyIncome, Age, DistanceFromHome) are measured on different scales (dollars, years, miles). Standardisation centres each variable to mean = 0 and scales to standard deviation = 1. This ensures that:

    • In linear regression, coefficients become comparable in terms of effect size.
    • Tree‑based methods like Random Forest are not affected by scale, but standardisation is kept for consistency across models.
    • No imputation is needed because the dataset contains zero missing values.
    • No outlier removal is applied because tree‑based models are robust and linear regression residuals can be examined afterwards.

These choices are standard in tidymodels pipelines and improve numerical stability and interpretability.

Data Split and Preprocessing

We use stratified splitting on MonthlyIncome to ensure that the income distribution is proportionally represented in both the training (75%) and test (25%) sets. The preprocessing recipe applies one-hot encoding to all categorical predictors and standardises all numeric predictors to zero mean and unit variance.

set.seed(100)
data_split <- initial_split(data, prop = 0.75, strata = MonthlyIncome)
data_train <- training(data_split)
data_test  <- testing(data_split)

cat("Training set:", nrow(data_train), "rows | Test set:", nrow(data_test), "rows")
## Training set: 1101 rows | Test set: 369 rows
data_prep <- recipe(MonthlyIncome ~ ., data = data_train) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%
  step_normalize(all_numeric_predictors()) %>%
  prep()

data_train_proc <- bake(data_prep, new_data = NULL)
data_test_proc  <- bake(data_prep, new_data = data_test)

Model Fitting and Evaluation

lm_fit <- linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression") %>%
  fit(MonthlyIncome ~ ., data = data_train_proc)

reg_results <- data_test_proc %>%
  bind_cols(predict(lm_fit, new_data = data_test_proc))

metric_set(rmse, rsq, mae)(reg_results,
                            truth    = MonthlyIncome,
                            estimate = .pred) %>%
  kbl(col.names = c("Metric", "Estimator", "Value"),
      caption = "Linear Regression — Test Set Performance (Standardised Scale)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Linear Regression — Test Set Performance (Standardised Scale)
Metric Estimator Value
rmse standard 1068.4599701
rsq standard 0.9506218
mae standard 821.6818294 
p9 <- ggplot(reg_results, aes(x = MonthlyIncome, y = .pred,
             text = paste0("Actual: ",    round(MonthlyIncome, 2),
                           "<br>Predicted: ", round(.pred, 2)))) +
  geom_point(alpha = 0.35, color = "#3498db", size = 1.2) +
  geom_abline(slope = 1, intercept = 0, color = "#e74c3c",
              linetype = "dashed", linewidth = 0.9) +
  labs(title = "Actual vs. Predicted Monthly Income (Standardised)",
       subtitle = "Points near the diagonal indicate accurate predictions",
       x = "Actual (Standardised)", y = "Predicted (Standardised)") +
  theme_minimal(base_size = 12)

ggplotly(p9, tooltip = "text") %>% config(displayModeBar = FALSE)

The R² value indicates the proportion of income variance explained by the model on held-out test data. Systematic curvature away from the diagonal would indicate non-linearities the linear model fails to capture — if present, a non-linear model such as Gradient Boosting would be the natural next step.

6.2 Attrition Classification — Random Forest

6.2.1 Random Forest Classifier (for Attrition Prediction)

Definition:
Random Forest is an ensemble learning method for classification (and regression) that constructs a large number of decision trees during training and outputs the mode of the classes (for classification) predicted by individual trees.

Purpose:
To predict a categorical outcome (Attrition: Yes/No) with high accuracy, robustness to overfitting, and the ability to capture complex, non‑linear interactions between predictors.

Logic:
The algorithm creates many bootstrap samples (random samples with replacement) from the training data. For each bootstrap sample, a decision tree is grown, but at each node only a random subset of predictors is considered for splitting. This decorrelates the trees. The final prediction is the majority vote across all trees.

Formula:
No closed‑form equation; the prediction for a new observation \(x\) is:

\[ \hat{f}(x) = \text{mode}\{ T_1(x), T_2(x), \dots, T_B(x) \} \]

where \(T_b\) is the \(b\)-th tree and \(B\) is the total number of trees.

Out‑of‑bag error estimate:

\[ \text{OOB error} = \frac{1}{n}\sum_{i=1}^n I\big(y_i \neq \hat{f}_{\text{OOB}}(x_i)\big) \]

Suggested use cases:
- When the outcome is categorical and classes are imbalanced (as with attrition)
- When there are many predictors with unknown or non‑linear relationships
- When high predictive performance is needed and interpretability of individual coefficients is not required (variable importance can still be extracted)
- When overfitting is a concern with single decision trees

Reasoning for method choice:
Employee attrition is a binary outcome with approximately 16% positive class (Yes) – a moderate imbalance. Relationships between attrition and predictors (e.g., income, distance, satisfaction) are likely non‑linear and interactive (e.g., the effect of low income may depend on job role). Random Forest handles these complexities automatically without manual feature engineering. It is robust to outliers and irrelevant predictors, and its out‑of‑bag error estimate provides an internal validation. The high recall (~99%) achieved in cross‑validation is critical for HR applications – missing a leaver is costly. Additionally, variable importance scores from Random Forest can identify the strongest drivers of attrition, guiding targeted retention policies.
The primary modelling objective is identifying which employees are at risk of leaving. We use a Random Forest classifier via the ranger engine — an ensemble method that constructs many decorrelated decision trees and aggregates their predictions. Random Forests handle mixed predictor types well, are robust to outliers, and provide variable importance scores interpretable in a business context.

Reconstructing the Attrition Factor

The preprocessing recipe encoded Attrition via one-hot encoding for the regression task. We rebuild the original factor from those encoded columns.

data_train_proc <- data_train_proc %>%
  mutate(Attrition = factor(
    ifelse(Attrition_No > 0, "Yes", "No"),
    levels = c("Yes", "No")
  )) %>%
  select(-any_of(c("Attrition_Yes", "Attrition_No")))

5-Fold Cross-Validation

Rather than relying on a single train/test split, we use 5-fold stratified cross-validation — the training data is partitioned into 5 class-balanced folds, and the model is trained and evaluated 5 times. This yields more reliable performance estimates and quantifies variance across data subsets.

set.seed(100)
folds <- vfold_cv(data_train_proc, v = 5, strata = Attrition)

rf_spec <- rand_forest() %>%
  set_engine("ranger") %>%
  set_mode("classification")

cv_metrics <- metric_set(accuracy, precision, recall, f_meas, roc_auc)
set.seed(100)
ml_fit <- workflow() %>%
  add_model(rf_spec) %>%
  add_formula(Attrition ~ .) %>%
  fit_resamples(
    resamples = folds,
    metrics   = cv_metrics,
    control   = control_resamples(save_pred = TRUE)
  )

Cross-Validated Performance

cv_res <- ml_fit %>% collect_metrics()

cv_res %>%
  select(.metric, mean, std_err) %>%
  mutate(mean = round(mean, 4), std_err = round(std_err, 5)) %>%
  kbl(col.names = c("Metric", "Mean", "Std. Error"),
      caption = "Random Forest — 5-Fold Cross-Validated Performance") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(which(cv_res$.metric == "roc_auc"), bold = TRUE, background = "#d5f5e3")
Random Forest — 5-Fold Cross-Validated Performance
Metric Mean Std. Error
accuracy 0.8574 0.00457
f_meas 0.9214 0.00245
precision 0.8607 0.00331
recall 0.9914 0.00215
roc_auc 0.7947 0.02747 
p10 <- cv_res %>% 
  mutate(.metric = toupper(.metric)) %>% 
  ggplot(aes(x = reorder(.metric, mean), y = mean, fill = .metric,
             text = paste0("Metric: ", .metric,
                           "<br>Mean: ", round(mean, 4),
                           "<br>Std Error: +/- ", round(std_err, 4)))) +
  geom_col(width = 0.5, show.legend = FALSE) +
  geom_errorbar(aes(ymin = mean - std_err, ymax = mean + std_err), width = 0.2) +
  scale_fill_brewer(palette = "Blues", direction = -1) +
  scale_y_continuous(limits = c(0, 1.1)) +
  coord_flip() +
  labs(title = "Random Forest - Cross-validated Metrics",
       subtitle = "Error bars = +/- 1 standard error across 5 folds",
       x = NULL, y = "score") +
  theme_minimal(base_size = 12)

ggplotly(p10, tooltip = "text") %>% config(displayModeBar = FALSE)
  • Accuracy (~85.7%): Overall correct classifications — partially inflated by class imbalance.
  • Recall (~99.2%): Of all employees who actually left, the model correctly identifies ~99%. This is the most critical metric for an HR early-warning system — missing a true leaver is more costly than a false alarm.
  • Precision (~86.0%): Of all employees flagged as leavers, ~86% actually left.
  • F1-Score (~92.1%): Harmonic mean of precision and recall — strong overall discriminative performance.
  • ROC-AUC (~79.7%): Good discrimination between leavers and stayers, achieved with default hyperparameters and no tuning.

Confusion Matrix

ml_fit %>%
  collect_predictions() %>%
  conf_mat(truth = Attrition, estimate = .pred_class) %>%
  autoplot(type = "heatmap") +
  scale_fill_gradient(low = "#d6eaf8", high = "#1a5276") +
  labs(title = "Confusion Matrix — Random Forest (5-Fold CV)") +
  theme_minimal(base_size = 12)

The model catches nearly all actual leavers (very few false negatives), at the cost of some false positives. In practical HR retention workflows, this trade-off is typically desirable — the cost of incorrectly flagging a stayer is far lower than the cost of missing a true leaver.

ROC Curve

The ROC curve plots the True Positive Rate (recall) against the False Positive Rate across all possible classification thresholds. A curve hugging the top-left corner indicates a near-perfect classifier; a diagonal line represents random guessing.

ml_fit %>%
  collect_predictions() %>%
  roc_curve(truth = Attrition, .pred_Yes) %>%
  ggplot(aes(x = 1 - specificity, y = sensitivity)) +
  geom_path(color = "#2980b9", linewidth = 1.2) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey60") +
  annotate("text", x = 0.6, y = 0.15, label = "AUC ≈ 0.797",
           size = 4.5, fontface = "bold", color = "#2980b9") +
  labs(title = "ROC Curve — Random Forest Attrition Classifier",
       x = "False Positive Rate (1 − Specificity)",
       y = "True Positive Rate (Recall)") +
  theme_minimal()

The curve lies well above the random-chance diagonal, confirming genuine discriminative power. The AUC of ~0.797 indicates that in approximately 79.7% of randomly drawn leaver/stayer pairs, the model assigns a higher attrition probability to the leaver. This performance is achieved entirely with default Random Forest hyperparameters — grid-search tuning of mtry, min_n, and trees would likely push the AUC above 0.85.

7. Conclusions

This report presented a structured, end-to-end investigation of employee attrition using statistical hypothesis testing and machine learning. The key findings are summarised below.

Finding Statistical Method Outcome
Leavers earn less than the population average Wilcoxon Signed-Rank Significant (p < 0.05)
Marital status is associated with attrition Chi-Squared Significant (p < 0.05)
Leavers live further from the office Mann-Whitney U Significant (p < 0.05)
Career stagnation differs by marital status Kruskal-Wallis Not significant (p ≈ 0.14)
Female leavers earn less at Job Level 3 Kruskal-Wallis + Dunn Significant (p < 0.05)
Random Forest attrition prediction 5-Fold CV AUC = 0.797, Recall = 0.992

Three actionable conclusions stand out. First, below-market compensation is a consistent, statistically confirmed driver of departure — pay benchmarking and targeted salary reviews should be prioritised. Second, the gender pay gap identified among mid-level leavers warrants a formal pay equity audit within Job Level 3. Third, the Random Forest model achieves near-perfect recall (~99%), making it suitable as a deployed early-warning system — an HR business partner could use individual-level attrition probability scores to initiate retention conversations before employees begin actively job-searching.

Future extensions of this work could incorporate hyperparameter tuning via tune_grid(), class-imbalance correction using SMOTE from the themis package, and explainability tools such as SHAP values to produce individual-level attrition risk explanations for HR practitioners.

8. References

  1. Krishna, S., & Sidharth, S. (2022). HR analytics: Employee attrition analysis using random forest. International Journal of Performability Engineering, 18(4), 275.

  2. Basha, A. M., Srivani, J., Ankaiah, B., Dadakalandar, U., & Srinivaslulu, T. (2020). HR analytics using R-machine learning algorithm: multiple linear regression analysis. International Journal of Innovative Technology and Exploring Engineering, 9, 1179-1183.

  3. Angulakshmi, M., Madhumithaa, N., Dokku, S. R., Pachar, S., Sneha, K., & Lenin, D. S. (2024, September). Predictive HR analytics: using machine learning to forecast workforce trends and needs. In 2024 7th International Conference on Contemporary Computing and Informatics (IC3I) (Vol. 7, pp. 1399-1405). IEEE.

  4. Çene, E. (n.d.). IST3110 Ders Notları: Üretken Yapay Zeka. Yıldız Teknik Üniversitesi, İstatistik Bölümü.

Report generated with R 4.4.2 | tidymodels 1.5.0 | ranger 0.18.0