This educational example predicts a synthetically assigned label, not actual employee performance. Collaboration metrics do not establish performance or causal drivers of success. This demonstration must not be used to rank employees or make employment decisions. We use a random forest to illustrate leakage-safe model selection and evaluation.
We’ll start by loading necessary libraries and importing data, then preprocess the data, build the model, and evaluate its performance.
Let’s begin by loading the required libraries and importing the dataset.
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(tidyverse)
library(vivainsights)
library(randomForest) # For fitting random forest model and extracting stats
library(caret) # For the final confusion matrix
library(pROC)
In R, once a package is loaded with library(package),
it’s not necessary to prefix the function with the package name when
calling it. However, if you want to use a function without loading the
entire package, or if there’s a naming conflict between functions from
different packages, you can use package::function() to
explicitly call the function. For clarity in this demonstration, we’ll
use this explicit notation to show which package each function comes
from.
The next step here is to load in the dataset, and then examine the data. In our local directory, we have a demo dataset that has a similar structure to a Person Query, with an additional 5-point scale ‘performance’ attribute that represents performance scores.
Synthetic fixture provenance.
_data/Top_Performers_Dataset_v2.csvis generated by the checked-ingenerate-demo-data.Rwith a fixed seed, invented parameters and obviousSYNTH_RF_IDs. No original records, hashes or fitted distributions are used. This replaces an earlier fixture with unverified provenance; see the topic README for the synthetic-data contract. It is not de-identified employee data.
import_query() imports the demo person query data, and
performs cleaning on the variable names. An alternative to this is to
use read_csv() (from the package readr,
which is part of tidyverse), which does the same thing
of reading in the input csv file.
# Set path to direct to where the Person Query is saved
raw_data <- vivainsights::import_query("_data/Top_Performers_Dataset_v2.csv")
# Examine the data
head(raw_data)
## # A tibble: 6 × 6
## PersonId Internal_network_size Collaboration_hours weekend_collaboratio…¹
## <chr> <int> <dbl> <dbl>
## 1 SYNTH_RF_0001 37 13.8 0.57
## 2 SYNTH_RF_0001 39 7.46 0.68
## 3 SYNTH_RF_0001 34 10.4 0.02
## 4 SYNTH_RF_0001 39 8.2 0.91
## 5 SYNTH_RF_0001 34 8.91 0.33
## 6 SYNTH_RF_0001 37 13.8 0.29
## # ℹ abbreviated name: ¹weekend_collaboration_hours
## # ℹ 2 more variables: After_hours_call_hours <dbl>, performance <int>
There are typically a number of data preparation and validation procedures involved before fitting a model, such as: - Handling missing values - Changing variable types - Handling outliers and unwanted data - Splitting data into training and test sets
In this notebook, we will assume that the dataset is in decent quality, and all that is required are the standard procedures of changing variable types and splitting data into train/test sets.
It is optional, but we convert the performance variable
into a binary variable (perform_cat), so we would yield a
classification model. This step is for demo purposes as there are more
use cases where the outcome variable is binary rather than ordinal or
continuous. We retain PersonId for now, as it is needed to
split the data by person, and drop it prior to fitting the model.
clean_data <-
raw_data %>%
mutate(perform_cat = ifelse(performance >= 4, 1, 0)) %>% # Create binary variable
mutate(perform_cat = factor(perform_cat)) %>% # ensures model is classification
select(-performance) # drop unnecessary columns
head(clean_data)
## # A tibble: 6 × 6
## PersonId Internal_network_size Collaboration_hours weekend_collaboratio…¹
## <chr> <int> <dbl> <dbl>
## 1 SYNTH_RF_0001 37 13.8 0.57
## 2 SYNTH_RF_0001 39 7.46 0.68
## 3 SYNTH_RF_0001 34 10.4 0.02
## 4 SYNTH_RF_0001 39 8.2 0.91
## 5 SYNTH_RF_0001 34 8.91 0.33
## 6 SYNTH_RF_0001 37 13.8 0.29
## # ℹ abbreviated name: ¹weekend_collaboration_hours
## # ℹ 2 more variables: After_hours_call_hours <dbl>, perform_cat <fct>
We first average the observed metrics within each person, so every modelled observation and every evaluation weight represents one person. Labels must be stable within a person. This is same-window classification, not a forecast of future performance.
We reserve approximately 20% of each class’s people for final testing and 20% for validation; the rest train the tuning candidates. Validation is part of the development population, not the final test set. At least three people in each class are required so each partition contains both classes. Smaller classes fail with an actionable error instead of producing an undefined ROC.
predictors <- c("Internal_network_size", "Collaboration_hours",
"weekend_collaboration_hours", "After_hours_call_hours")
stopifnot(all(c("PersonId", "perform_cat", predictors) %in% names(clean_data)),
!anyNA(clean_data[c("PersonId", "perform_cat", predictors)]),
all(vapply(clean_data[predictors], function(x) all(is.finite(x)), logical(1))))
labels <- clean_data %>% distinct(PersonId, perform_cat)
if (anyDuplicated(labels$PersonId)) stop("Each person must have a stable outcome label.")
person_data <- clean_data %>%
group_by(PersonId, perform_cat) %>%
summarise(across(all_of(predictors), mean), .groups = "drop")
split_people <- function(data, seed = 123L) {
counts <- table(factor(data$perform_cat, levels = c("0", "1")))
if (anyDuplicated(data$PersonId) || anyNA(data[c("PersonId", "perform_cat")]) ||
!all(as.character(data$perform_cat) %in% c("0", "1")) || any(counts < 3L)) {
stop("Need unique people and at least three people in each binary class.")
}
set.seed(seed)
ids <- split(data$PersonId, data$perform_cat, drop = TRUE)
parts <- lapply(ids, function(x) {
x <- x[sample.int(length(x))]
n_holdout <- max(1L, floor(length(x) * 0.2))
list(test = x[seq_len(n_holdout)],
validation = x[n_holdout + seq_len(n_holdout)],
train = x[-seq_len(2L * n_holdout)])
})
lapply(c("train", "validation", "test"), function(p) {
unlist(lapply(parts, `[[`, p), use.names = FALSE)
}) |> setNames(c("train", "validation", "test"))
}
partition <- split_people(person_data)
stopifnot(!anyDuplicated(unlist(partition)),
setequal(unlist(partition), person_data$PersonId))
train_ids <- partition$train
validation_ids <- partition$validation
test_ids <- partition$test
model_columns <- c("perform_cat", predictors)
train_df <- person_data %>% filter(PersonId %in% train_ids) %>% select(all_of(model_columns))
validation_df <- person_data %>% filter(PersonId %in% validation_ids) %>% select(all_of(model_columns))
development_df <- person_data %>% filter(!PersonId %in% test_ids) %>% select(all_of(model_columns))
It is good practice to double check what is in your training data frame prior to running the model, to ensure that you are not including any unwanted predictors by mistake:
names(train_df)
## [1] "perform_cat" "Internal_network_size"
## [3] "Collaboration_hours" "weekend_collaboration_hours"
## [5] "After_hours_call_hours"
ntree controls the number of trees,
nodesize the minimum terminal node size, and
mtry the candidate predictors at each split. This small
joint grid selects the highest validation AUC, with deterministic ties
favouring the earlier grid row. The positive class and ROC direction are
fixed in advance, as is the final classification threshold of 0.5. No
test labels or test predictions are used here.
auc_score <- function(actual, probability) {
as.numeric(pROC::auc(pROC::roc(actual, probability,
levels = c("0", "1"), direction = "<", quiet = TRUE)))
}
grid <- expand.grid(ntree = c(100L, 300L), nodesize = c(1L, 5L),
mtry = seq_along(predictors))
tuning_results <- purrr::map_dfr(seq_len(nrow(grid)), function(i) {
set.seed(1000L + i)
candidate <- randomForest::randomForest(
perform_cat ~ ., data = train_df,
ntree = grid$ntree[i], nodesize = grid$nodesize[i], mtry = grid$mtry[i]
)
probability <- predict(candidate, newdata = validation_df, type = "prob")[, "1"]
cbind(grid[i, ], validation_auc = auc_score(validation_df$perform_cat, probability))
})
selected <- tuning_results[which.max(tuning_results$validation_auc), ]
tuning_results
## ntree nodesize mtry validation_auc
## 1 100 1 1 0.8557143
## 2 300 1 1 0.8457143
## 3 100 5 1 0.8471429
## 4 300 5 1 0.8400000
## 5 100 1 2 0.8071429
## 6 300 1 2 0.8185714
## 7 100 5 2 0.8028571
## 8 300 5 2 0.8300000
## 9 100 1 3 0.8257143
## 10 300 1 3 0.7914286
## 11 100 5 3 0.8014286
## 12 300 5 3 0.8128571
## 13 100 1 4 0.7742857
## 14 300 1 4 0.7742857
## 15 100 5 4 0.7714286
## 16 300 5 4 0.8000000
selected
## ntree nodesize mtry validation_auc
## 1 100 1 1 0.8557143
Validation AUC is a model-selection score, not an unbiased final performance estimate. With few people in a class it is especially unstable. The fixture is large enough to demonstrate the workflow, not to justify real-world accuracy.
The next step is to fit the random forest model, with
randomForest() from the randomForest
package.
With randomForest(), it is possible to specify your
variables either in the (i) formula style or by (ii) supplying data
frames of predictors and outcome. The following example shows the
formula style.
Note that randomForest() comes with many default
parameters, which you can find out more from
its official reference manual.
After selection, refit once on all development people (training plus
validation). importance = TRUE enables the interpretation
examples below.
set.seed(2026)
# Build the random forest model
rf <- randomForest(
formula = perform_cat ~ .,
data = development_df,
ntree = selected$ntree, nodesize = selected$nodesize, mtry = selected$mtry,
importance = TRUE # to allow importance to be calculated afterwards
)
rf
##
## Call:
## randomForest(formula = perform_cat ~ ., data = development_df, ntree = selected$ntree, nodesize = selected$nodesize, mtry = selected$mtry, importance = TRUE)
## Type of random forest: classification
## Number of trees: 100
## No. of variables tried at each split: 1
##
## OOB estimate of error rate: 26.09%
## Confusion matrix:
## 0 1 class.error
## 0 83 19 0.1862745
## 1 23 36 0.3898305
Here are some bullet points on how to interpret the model summary:
Type of random forest: classification: This tells us that the model is used for classification tasks, not regression.
Number of trees: 100; candidate variables per split: 1.
OOB error: 26.1%. This internal development diagnostic is not the final held-out evaluation.
OOB confusion matrix: the printed
rf$confusion has actual classes in rows and predictions in
columns; class.error reports each class’s misclassification
rate. Counts and errors are computed from this run, not a frozen
example.
The model reflects the invented signal in the generator. No level of synthetic accuracy is evidence that collaboration metrics measure employee success.
If no errors or warnings pop up, then the first iteration of the model is trained. The next step is to understand the model, and then to interpret and evaluate its outputs.
# Open the reserved test population only after the selected model is frozen.
test_df <- person_data %>% filter(PersonId %in% test_ids) %>% select(all_of(model_columns))
test_probability <- predict(object = rf, newdata = test_df, type = "prob")[, "1"]
pred <- factor(ifelse(test_probability >= 0.5, "1", "0"), levels = c("0", "1"))
# Attach this to a data frame for easy referencing
test_df_with_pred <-
test_df %>%
mutate(predictions = pred)
# Extract and assign variables
actual <- test_df_with_pred$perform_cat
predicted <- test_df_with_pred$predictions
# Create a confusion matrix
cm <- confusionMatrix(predicted, actual, positive = "1")
# Accuracy
accuracy <- cm$overall['Accuracy']
# Precision, Recall, and F1 Score
precision <- cm$byClass['Pos Pred Value']
recall <- cm$byClass['Sensitivity']
tp <- cm$table["1", "1"]
fp <- cm$table["1", "0"]
fn <- cm$table["0", "1"]
f1_score <- 2 * tp / (2 * tp + fp + fn)
test_auc <- auc_score(actual, test_probability)
# Print the metrics
data.frame(
statistic = c("Accuracy", "Precision", "Recall", "F1 Score", "AUC"),
value = c(accuracy, precision, recall, f1_score, test_auc)
)
## statistic value
## 1 Accuracy 0.6923077
## 2 Precision 0.5714286
## 3 Recall 0.5714286
## 4 F1 Score 0.5714286
## 5 AUC 0.8114286
The single test prediction vector supplies all final metrics. If no positives are predicted, precision is undefined (NA); F1 is zero when positives are missed. These are person-level, same-window metrics without uncertainty intervals. Do not tune further against these test results.
We also generated a number of metrics for assessing the model. The first four metrics below range between 0 and 1, and an idealistic perfect model would return 1, meaning that it makes no errors of the type:
Accuracy: This is the ratio of correct predictions to the total number of predictions. It’s a good measure when the target variable classes in the data are nearly balanced. However, it can be misleading if the classes are imbalanced.
Precision: Precision is the ratio of true positives (correctly predicted positive observations) to the total predicted positives. It’s a measure of a classifier’s exactness. A low precision indicates a high number of false positives (Type I errors).
Recall (Sensitivity): Recall is the ratio of true positives to the total actual positives. It’s a measure of a classifier’s completeness. A low recall indicates a high number of false negatives (Type II errors).
F1 Score: The F1 Score is the weighted average of Precision and Recall. It tries to balance the two metrics. It’s a good measure to use if you need to seek a balance between Precision and Recall and there is an uneven class distribution. It is given by:
F1 = 2 * (precision * recall) / (precision + recall)
All of these metrics can be calculated directly from a
Confusion Matrix, which is a table that is often used
to describe the performance of a classification model on a set of data
for which the true values are known. It contains information about
actual and predicted classifications done by the classifier. It’s a good
way to visualize the performance of the model. In R, this is generated
with the confusionMatrix() function from
caret.
The choice of metric depends on your business objective. For example, if the cost of having false positives is high, the strategy might be to optimize for precision; this arguably applies to a top performers use case, where it is preferred that the model predicts fewer top performers. If the cost of missing positives (having false negatives) is high, the strategy might be to optimize for recall, which could be more relevant for an attrition use case.
In R, you can call the confusion matrix with cm$table,
if you have assigned the output of confusionMatrix() to
cm:
cm$table
## Reference
## Prediction 0 1
## 0 19 6
## 1 6 8
See below for a guide on how to interpret the confusion matrix:
A confusion matrix is a table that is often used to describe the
performance of a classification model on a set of data for which the
true values are known. In binary classification, the confusion matrix is
a 2x2 matrix. caret::confusionMatrix() places predictions
in rows and actual reference classes in columns. Here’s how to interpret
it:
So, the confusion matrix looks like this:
| Actual Negative | Actual Positive | |
|---|---|---|
| Predicted Negative | True Negative (TN) | False Negative (FN) |
| Predicted Positive | False Positive (FP) | True Positive (TP) |
The diagonal elements represent the number of points for which the predicted label is equal to the true label, while off-diagonal elements are those that are mislabeled by the classifier. The higher the diagonal values of the confusion matrix, the better, indicating many correct predictions.
One of the major outputs of the Random Forest model is feature importance.
Feature importance can be calculated with
randomForest::importance(), which allows you to return two
types of calculations.
Impurity-based Feature Importance: Mean Decrease in Gini impurity (MDG): MDG is the total decrease in node impurities from splitting on the variable, averaged over all trees. For classification, the node impurity is measured by the Gini index. For regression, it is measured by residual sum of squares. This is sometimes called Mean Decrease in Impurity (MDI). This method is fast to compute and does not require a separate validation set or model re-fitting. However, it tends to inflate the importance of continuous features or high-cardinality categorical variables. It is also biased towards features with more categories.
Permutation-based Feature Importance: Mean Decrease in Accuracy (MDA): MDA is computed from permuting OOB data: For each tree, the prediction error on the out-of-bag portion of the data is recorded (error rate for classification, MSE for regression). Then the same is done after permuting each predictor variable. The difference between the two are then averaged over all trees, and normalized by the standard deviation of the differences. If the standard deviation of the differences is equal to 0 for a variable, the division is not done (but the average is almost always equal to 0 in that case). This method is more reliable and has less bias towards continuous or high-cardinality features, but it is computationally expensive as it requires re-fitting the model for each feature.
It is worth noting that both methods aim to capture the importance of features, but they focus on different aspects. MDG emphasizes impurity reduction during tree construction, while MDA directly considers the impact on model accuracy.
Here is an example of MDG, as well as how to visualise it (using
varImpPlot()):
randomForest::importance(rf, type = 2)
## MeanDecreaseGini
## Internal_network_size 23.31580
## Collaboration_hours 20.68764
## weekend_collaboration_hours 14.52797
## After_hours_call_hours 16.19708
randomForest::varImpPlot(rf, type = 2)
And here is the equivalent for MDA:
randomForest::importance(rf, type = 1)
## MeanDecreaseAccuracy
## Internal_network_size 8.156576
## Collaboration_hours 9.072448
## weekend_collaboration_hours -1.160662
## After_hours_call_hours 3.586712
randomForest::varImpPlot(rf, type = 1)
For those comparing the results between Python’s
scikit-learn library and R’s
randomForest package, note that
feature_importances_ in scikit-learn by
default computes Mean Decrease Impurity (MDI), and not Mean
Decrease Accuracy (MDA). To compute MDA,
scikit-learn uses a separate function
permutation_importance() to do so. The main difference is
that in R, this is being controlled by the importance
argument within randomForest itself.