Top Performers Model

Introduction

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.

Set-up

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.csv is generated by the checked-in generate-demo-data.R with a fixed seed, invented parameters and obvious SYNTH_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>

Data Preparation

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"

Tune only within the development population

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.

Fitting the selected model

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.

Evaluating the model

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:

  • The first row of the matrix represents predictions for the negative class.
  • The second row of the matrix represents predictions for the positive class.
  • The first column represents actual observations in the negative class.
  • The second column represents actual observations in the positive class.

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)
  • True Positives (TP): These are cases in which we predicted yes (positive), and the actual was also yes (positive).
  • True Negatives (TN): We predicted no (negative), and the actual was also no (negative).
  • False Positives (FP): We predicted yes (positive), but the actual was no (negative). Also known as “Type I error”.
  • False Negatives (FN): We predicted no (negative), but the actual was yes (positive). Also known as “Type II error”.

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.

Variable importance

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.

  1. 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.

  2. 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.

Calculating and Visualising Feature Importance

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.