# 1. Import libraries
import numpy as np
import pandas as pd
from IPython.display import display
import matplotlib.pyplot as plt
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.ensemble import RandomForestClassifier
from sklearn.pipeline import Pipeline
from sklearn.metrics import (confusion_matrix, classification_report, accuracy_score, roc_curve, roc_auc_score)
from sklearn.inspection import permutation_importance
RANDOM_STATE = 42
N_TREES = 100
MTRY = 2
print('Libraries imported successfully.')
# 2. Read the supplied Employee Churn dataset
empdata = pd.read_csv('Employee Churn.csv')
print('Shape:', empdata.shape)
display(empdata.head())
print('\nData types:')
print(empdata.dtypes)
print('\nStatus counts:')
print(empdata['status'].value_counts().sort_index())
# 3. Select the variables used in the R case study
features = ['function', 'exp', 'gender', 'source']
target = 'status'
X = empdata[features].copy()
y = empdata[target].astype(int).copy()
print('Predictors:', features)
print('Target:', target)
print('\nUnique values in categorical predictors:')
for col in features:
    print(f'{col}: {X[col].unique()}')
# 4. Build the Random Forest
# churn_rf <- randomForest(
# status ~ function. + exp + gender + source, data = empdata,
# mtry = 2,
# ntree = 100, importance = TRUE, cutoff = c(0.6, 0.4)
# )
# One-hot encode categorical predictors.
categorical_features = features
preprocessor = ColumnTransformer( transformers=[
('cat', OneHotEncoder(handle_unknown='ignore'), categorical_features)
]
)
rf = RandomForestClassifier( n_estimators=N_TREES, max_features=MTRY, bootstrap=True, oob_score=True,
random_state=RANDOM_STATE,
n_jobs=-1
)
churn_rf = Pipeline([ ('preprocessor', preprocessor), ('model', rf)
])
churn_rf.fit(X, y)
print('Random Forest fitted.')
print('Number of trees:', N_TREES)
print('max_features (mtry):', MTRY)
print('OOB score:', round(rf.oob_score_, 4))
print('OOB error rate:', round(1 - rf.oob_score_, 4))
# 6. Predictions and confusion matrices
proba = churn_rf.predict_proba(X)[:, 1]
# Standard classification threshold
pred_050 = (proba >= 0.50).astype(int)
# R-style binary cutoff c(0.6, 0.4): status 1 is selected when
# P(status=1) / 0.4 > P(status=0) / 0.6, which is equivalent to P(status=1) > 0.4.
pred_040 = (proba >= 0.40).astype(int)
print('Confusion matrix — threshold 0.50')
print(confusion_matrix(y, pred_050))
print('\nAccuracy:', round(accuracy_score(y, pred_050), 4))
print('\nConfusion matrix — R-style cutoff 0.40 for status=1')
print(confusion_matrix(y, pred_040))
print('\nAccuracy:', round(accuracy_score(y, pred_040), 4))
print('\nClassification report — R-style cutoff')
print(classification_report(y, pred_040, digits=4))
oob_proba = rf.oob_decision_function_[:, 1]
oob_pred_050 = (oob_proba >= 0.50).astype(int)
oob_pred_040 = (oob_proba >= 0.40).astype(int)
print('OOB confusion matrix — threshold 0.50')
print(confusion_matrix(y, oob_pred_050))
print('\nOOB confusion matrix — R-style cutoff 0.40')
print(confusion_matrix(y, oob_pred_040))
print('\nOOB error at 0.50:', round(1 - accuracy_score(y, oob_pred_050), 4))
print('OOB error at 0.40:', round(1 - accuracy_score(y, oob_pred_040), 4))
# Aggregate impurity-based importance from one-hot encoded columns back to original variables.
encoder = churn_rf.named_steps['preprocessor'].named_transformers_['cat']
encoded_names = encoder.get_feature_names_out(categorical_features)
encoded_importance = rf.feature_importances_
importance_df = pd.DataFrame({ 'encoded_feature': encoded_names, 'importance': encoded_importance
})
importance_df['original_variable'] = np.repeat(categorical_features, [len(categories) for categories in encoder.categories_])
grouped_gini = (
importance_df.groupby('original_variable', as_index=False)['importance']
.sum()
.sort_values('importance', ascending=False)
)
display(grouped_gini)
plt.figure(figsize=(8, 5))
plt.barh(grouped_gini['original_variable'], grouped_gini['importance'])
plt.gca().invert_yaxis()
plt.xlabel('Aggregated impurity importance')
plt.ylabel('Predictor')
plt.title('Random Forest Variable Importance')
plt.tight_layout()
plt.show()
# Permutation importance on the ORIGINAL dataframe columns.
# The pipeline accepts the four original categorical columns, so this gives
# a variable-level importance measure rather than one importance per dummy variable.
perm = permutation_importance(
churn_rf, X, y, n_repeats=30,
random_state=RANDOM_STATE, scoring='accuracy',
n_jobs=-1
)
perm_df = pd.DataFrame({ 'variable': X.columns,
'mean_decrease_accuracy': perm.importances_mean, 'std': perm.importances_std
}).sort_values('mean_decrease_accuracy', ascending=False)
display(perm_df)
plt.figure(figsize=(8, 5))
plt.barh(perm_df['variable'], perm_df['mean_decrease_accuracy'], xerr=perm_df['std'])
plt.gca().invert_yaxis()
plt.xlabel('Mean decrease in accuracy after permutation')
plt.ylabel('Predictor')
plt.title('Permutation Importance — Mean Decrease Accuracy Analogue')
plt.tight_layout()
plt.show()
# ROC curve using the fitted model's status=1 probabilities.
fpr, tpr, thresholds = roc_curve(y, proba)
auc_value = roc_auc_score(y, proba)
print('AUC =', round(auc_value, 7))
plt.figure(figsize=(7, 6))
plt.plot(fpr, tpr, label=f'Random Forest (AUC = {auc_value:.4f})')
plt.plot([0, 1], [0, 1], linestyle='--', label='Random classifier')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve — Employee Churn Random Forest')
plt.legend()
plt.tight_layout()
plt.show()
new_employee = pd.DataFrame({ 'function': ['CS'],
'exp': ['<3'],
'gender': ['M'],
'source': ['external']
})
new_probability = churn_rf.predict_proba(new_employee)[0, 1]
new_prediction = int(new_probability >= 0.40) # R-style cutoff
print('Probability of leaving within 18 months:', round(new_probability, 4))
print('Predicted status using R-style cutoff:', new_prediction)
print('Interpretation:', 'Likely to leave' if
new_prediction == 1 else 'Likely to stay')
prediction_output = empdata.copy()
prediction_output['prob_status_1'] = proba
prediction_output['predicted_status_0.40'] = pred_040
display(prediction_output.head(20))
