Despite the AI hype, many tech companies still rely heavily on machine learning to power critical applications, from personalized recommendations to fraud detection.\u00a0<\/p>\n
I\u2019ve seen firsthand how undetected drifts can result in significant costs \u2014 missed fraud detection, lost revenue, and suboptimal business outcomes, just to name a few. So, it\u2019s crucial to have robust monitoring in place if your company has deployed or plans to deploy machine learning models into production.<\/p>\n
Undetected Model Drift<\/a> can lead to significant financial losses, operational inefficiencies, and even damage to a company\u2019s reputation. To mitigate these risks, it\u2019s important to have effective model monitoring, which involves:<\/p>\n A well-implemented monitoring system can help identify issues early, saving considerable time, money, and resources.<\/p>\n In this comprehensive guide, I\u2019ll provide a framework on how to think about and implement effective Model Monitoring<\/a>, helping you stay ahead of potential issues and ensure stability and reliability of your models in production.<\/p>\n Score drift refers to a gradual change in the distribution of model scores. If left unchecked, this could lead to a decline in model performance<\/strong>, making the model less accurate over time.<\/p>\n On the other hand, feature drift occurs when one or more features experience changes in the distribution. These changes in feature values can affect the underlying relationships that the model has learned, and ultimately lead to inaccurate model predictions.<\/p>\n To model real-world fraud detection challenges, I created a synthetic dataset with five financial transaction features.<\/p>\n The reference dataset<\/strong> represents the original distribution, while the production dataset<\/strong> introduces shifts to simulate an increase in high-value transactions without PIN verification on newer accounts,<\/strong> indicating an increase in fraud.<\/p>\n Each feature has different underlying distributions:<\/p>\n To approximate model scores, I randomly assigned weights to these features and applied a sigmoid function to constrain predictions between 0 to 1. This mimics how a logistic regression fraud model generates risk scores.<\/p>\n As shown in the plot below:<\/p>\n This setup allows us to analyze how feature drift impacts model scores in production.<\/p>\n To monitor model scores, I used population stability index (PSI) to measure how much model score distribution has shifted over time.<\/p>\n PSI works by binning continuous model scores and comparing the proportion of scores in each bin between the reference and production datasets. It compares the differences in proportions and their logarithmic ratios to compute a single summary statistic to quantify the drift.<\/p>\n Python<\/a> implementation:<\/strong><\/p>\n Below is a summary of how to interpret PSI values:<\/p>\n The PSI value of 0.6374<\/strong> suggests a significant drift between our reference and production datasets. This aligns with the histogram of model score distributions, which visually confirms the shift towards higher scores in production \u2014 indicating an increase in risky transactions.<\/strong><\/p>\n Kolmogorov-Smirnov test for numeric features<\/strong><\/p>\n The Kolmogorov-Smirnov (K-S) test is my preferred method for detecting drift in numeric features, because it is non-parametric,<\/strong> meaning it doesn\u2019t assume a normal distribution.<\/p>\n The test compares a feature\u2019s distribution in the reference and production datasets by measuring the maximum difference between the empirical cumulative distribution functions (ECDFs). The resulting K-S statistic ranges from 0 to 1:<\/p>\n Python implementation:<\/strong><\/p>\n Below are ECDF charts of the four numeric features in our dataset:<\/p>\n Let\u2019s look at the account age feature as an example: the x-axis represents account age (0-50 months), while the y-axis shows the ECDF for both reference and production datasets. The production dataset skews towards newer accounts, as it has a larger proportion of observations with lower account ages.<\/p>\n Chi-Square test for categorical features<\/strong><\/p>\n To detect shifts in categorical and boolean features, I like to use the Chi-Square test.<\/p>\n This test compares the frequency distribution of a categorical feature in the reference and production datasets, and returns two values:<\/p>\n Python implementation:<\/strong><\/p>\n The Chi-Square statistic of 57.31 with a p-value of 3.72e-14 confirms a large shift in our categorical feature, Spearman Correlation<\/a> for shifts in pairwise interactions<\/strong><\/p>\n In addition to monitoring individual feature shifts, it\u2019s important to track shifts in relationships or interactions between features<\/strong>, known as multivariate shifts. Even if the distributions of individual features remain stable, multivariate shifts can signal meaningful differences in the data.<\/p>\n By default, Pandas\u2019 To capture this, we use Spearman correlation<\/strong> to measure monotonic relationships<\/strong> between features. It captures whether features change together<\/strong> in a consistent direction, even if their relationship isn\u2019t strictly linear.<\/p>\n To assess shifts in feature relationships, we compare:<\/p>\n Python implementation:<\/strong><\/p>\n Example: Change in correlation<\/strong><\/p>\n Now, let\u2019s look at the correlation between The absolute difference in correlation indicates a shift<\/strong> in the relationship between There used to be no relationship between these two features, but the production dataset indicates that there is now a weak negative correlation, meaning that newer accounts have higher transaction amounts. This is spot on!<\/p>\n Autoencoder for complex, high-dimensional multivariate shifts<\/strong><\/p>\n In addition to monitoring pairwise interactions, we can also look for shifts across more dimensions in the data.<\/p>\n Autoencoders are powerful tools for detecting high-dimensional multivariate shifts<\/strong>, where multiple features collectively change in ways that may not be apparent from looking at individual feature distributions or pairwise correlations.<\/p>\n An autoencoder is a neural network that learns a compressed representation of data through two components:<\/p>\n To detect shifts, we compare the reconstructed output<\/strong> to the original input<\/strong> and compute the reconstruction loss<\/strong>.<\/p>\n Unlike traditional drift metrics that focus on individual features or pairwise relationships<\/strong>, autoencoders capture complex, non-linear dependencies<\/strong> across multiple variables simultaneously.<\/p>\n Python implementation:<\/strong><\/p>\n The charts below show the distribution of reconstruction loss between both datasets.<\/p>\n The production dataset has a higher mean reconstruction error than that of the reference dataset, indicating a shift in the overall data. This aligns with the changes in the production dataset with a higher number of newer accounts with high-value transactions.<\/p>\n Model monitoring is an essential, yet often overlooked, responsibility for data scientists and machine learning engineers.<\/p>\n All the statistical methods led to the same conclusion, which aligns with the observed shifts in the data: they detected a trend in production towards newer accounts making higher-value transactions<\/strong>. This shift resulted in higher model scores, signaling an increase in potential fraud.<\/p>\n In this post, I covered techniques for detecting drift on three different levels:<\/p>\n These are just a few of the techniques I rely on for comprehensive monitoring \u2014 there are plenty of other equally valid statistical methods that can also detect drift effectively.<\/p>\n Detected shifts often point to underlying issues that warrant further investigation. The root cause could be as serious as a data collection bug, or as minor as a time change like daylight savings time adjustments.<\/p>\n There are also fantastic python packages, like evidently.ai<\/a>, that automate many of these comparisons. However, I believe there\u2019s significant value in deeply understanding the statistical techniques behind drift detection, rather than relying solely on these tools.<\/p>\n What\u2019s the model monitoring process like at places you\u2019ve worked?<\/p>\n Want to build your AI skills?<\/strong><\/p>\n\n
What\u2019s the difference between feature drift and score drift?<\/h2>\n
Simulating score shifts<\/h2>\n
\n
\n
![\"\"](\"https:\/\/i0.wp.com\/towardsdatascience.com\/wp-content\/uploads\/2025\/03\/Model-drift-1.png?resize=768%2C863&ssl=1\")
\n
![\"\"](\"https:\/\/i0.wp.com\/towardsdatascience.com\/wp-content\/uploads\/2025\/03\/Model-drift-2.png?resize=564%2C646&ssl=1\")
\n
![\"\"](\"https:\/\/i0.wp.com\/towardsdatascience.com\/wp-content\/uploads\/2025\/03\/Model-drift-3.png?resize=564%2C646&ssl=1\")
Detecting model score drift using PSI<\/h3>\n
# Define function to calculate PSI given two datasets\ndef calculate_psi(reference, production, bins=10):\n # Discretize scores into bins\n min_val, max_val = 0, 1\n bin_edges = np.linspace(min_val, max_val, bins + 1)\n\n # Calculate proportions in each bin\n ref_counts, _ = np.histogram(reference, bins=bin_edges)\n prod_counts, _ = np.histogram(production, bins=bin_edges)\n\n ref_proportions = ref_counts \/ len(reference)\n prod_proportions = prod_counts \/ len(production)\n \n # Avoid division by zero\n ref_proportions = np.clip(ref_proportions, 1e-8, 1)\n prod_proportions = np.clip(prod_proportions, 1e-8, 1)\n\n # Calculate PSI for each bin\n psi = np.sum((ref_proportions - prod_proportions) * np.log(ref_proportions \/ prod_proportions))\n\n return psi\n \n# Calculate PSI\npsi_value = calculate_psi(ref_data['model_score'], prod_data['model_score'], bins=10)\nprint(f\"PSI Value: {psi_value}\")<\/code><\/pre>\n\n
![\"\"](\"https:\/\/i0.wp.com\/towardsdatascience.com\/wp-content\/uploads\/2025\/03\/Model-drift-4.png?resize=848%2C544&ssl=1\")
Detecting feature drift<\/h2>\n
\n
# Create an empty dataframe\nks_results = pd.DataFrame(columns=['Feature', 'KS Statistic', 'p-value', 'Drift Detected'])\n\n# Loop through all features and perform the K-S test\nfor col in numeric_cols:\n ks_stat, p_value = ks_2samp(ref_data[col], prod_data[col])\n drift_detected = p_value < 0.05\n\t\t\n\t\t# Store results in the dataframe\n ks_results = pd.concat([\n ks_results,\n pd.DataFrame({\n 'Feature': [col],\n 'KS Statistic': [ks_stat],\n 'p-value': [p_value],\n 'Drift Detected': [drift_detected]\n })\n ], ignore_index=True)\n\n<\/code><\/pre>\n![\"\"](\"https:\/\/i0.wp.com\/towardsdatascience.com\/wp-content\/uploads\/2025\/03\/Model-drift-5-1024x681.png?resize=1024%2C681&ssl=1\")
\n
# Create empty dataframe with corresponding column names\nchi2_results = pd.DataFrame(columns=['Feature', 'Chi-Square Statistic', 'p-value', 'Drift Detected'])\n\nfor col in categorical_cols:\n # Get normalized value counts for both reference and production datasets\n ref_counts = ref_data[col].value_counts(normalize=True)\n prod_counts = prod_data[col].value_counts(normalize=True)\n\n # Ensure all categories are represented in both\n all_categories = set(ref_counts.index).union(set(prod_counts.index))\n ref_counts = ref_counts.reindex(all_categories, fill_value=0)\n prod_counts = prod_counts.reindex(all_categories, fill_value=0)\n\n # Create contingency table\n contingency_table = np.array([ref_counts * len(ref_data), prod_counts * len(prod_data)])\n\n # Perform Chi-Square test\n chi2_stat, p_value, _, _ = chi2_contingency(contingency_table)\n drift_detected = p_value < 0.05\n\n # Store results in chi2_results dataframe\n chi2_results = pd.concat([\n chi2_results,\n pd.DataFrame({\n 'Feature': [col],\n 'Chi-Square Statistic': [chi2_stat],\n 'p-value': [p_value],\n 'Drift Detected': [drift_detected]\n })\n ], ignore_index=True)<\/code><\/pre>\nEntered PIN<\/code>. This finding aligns with the histogram below, which visually illustrates the shift:<\/p>\n![\"\"](\"https:\/\/i0.wp.com\/towardsdatascience.com\/wp-content\/uploads\/2025\/03\/Model-drift-6.png?resize=375%2C433&ssl=1\")
Detecting multivariate shifts<\/h2>\n
.corr()<\/code> function calculates Pearson correlation, which only captures linear relationships between variables. However, relationships between features are often non-linear<\/strong> yet still follow a consistent trend.<\/p>\n\n
ref_corr<\/code>): Captures historical feature relationships in the reference dataset.<\/li>\nprod_corr<\/code>): Captures new feature relationships in production.<\/li>\n# Calculate correlation matrices\nref_corr = ref_data.corr(method='spearman')\nprod_corr = prod_data.corr(method='spearman')\n\n# Calculate correlation difference\ncorr_diff = abs(ref_corr - prod_corr)<\/code><\/pre>\ntransaction_amount<\/code> and account_age_in_months<\/code>:<\/p>\n\n
ref_corr<\/code>, the correlation is 0.00095, indicating a weak relationship between the two features.<\/li>\nprod_corr<\/code>, the correlation is -0.0325, indicating a weak negative correlation.<\/li>\ntransaction_amount<\/code> and account_age_in_months<\/code>.<\/p>\n\n
\n
ref_features = ref_data[numeric_cols + categorical_cols]\nprod_features = prod_data[numeric_cols + categorical_cols]\n\n# Normalize the data\nscaler = StandardScaler()\nref_scaled = scaler.fit_transform(ref_features)\nprod_scaled = scaler.transform(prod_features)\n\n# Split reference data into train and validation\nnp.random.shuffle(ref_scaled)\ntrain_size = int(0.8 * len(ref_scaled))\ntrain_data = ref_scaled[:train_size]\nval_data = ref_scaled[train_size:]\n\n# Build autoencoder\ninput_dim = ref_features.shape[1]\nencoding_dim = 3 \n# Input layer\ninput_layer = Input(shape=(input_dim, ))\n# Encoder\nencoded = Dense(8, activation=\"relu\")(input_layer)\nencoded = Dense(encoding_dim, activation=\"relu\")(encoded)\n# Decoder\ndecoded = Dense(8, activation=\"relu\")(encoded)\ndecoded = Dense(input_dim, activation=\"linear\")(decoded)\n# Autoencoder\nautoencoder = Model(input_layer, decoded)\nautoencoder.compile(optimizer=\"adam\", loss=\"mse\")\n\n# Train autoencoder\nhistory = autoencoder.fit(\n train_data, train_data,\n epochs=50,\n batch_size=64,\n shuffle=True,\n validation_data=(val_data, val_data),\n verbose=0\n)\n\n# Calculate reconstruction error\nref_pred = autoencoder.predict(ref_scaled, verbose=0)\nprod_pred = autoencoder.predict(prod_scaled, verbose=0)\n\nref_mse = np.mean(np.power(ref_scaled - ref_pred, 2), axis=1)\nprod_mse = np.mean(np.power(prod_scaled - prod_pred, 2), axis=1)<\/code><\/pre>\n![\"\"](\"https:\/\/i0.wp.com\/towardsdatascience.com\/wp-content\/uploads\/2025\/03\/Model-drift-7-1024x418.png?resize=1024%2C418&ssl=1\")
Summarizing<\/h2>\n
\n
\n
![\"\ud83d\udc49\ud83c\udffb\"](\"https:\/\/i0.wp.com\/s.w.org\/images\/core\/emoji\/15.0.3\/72x72\/1f449-1f3fb.png?ssl=1\")
I run the