Final Project - Group 4¶
Healthcare Fraud Detection Project¶
1. Case Description¶
Motivation and topic. Healthcare fraud is one of the largest and fastest growing forms of financial crime, estimated to cost payers worldwide hundreds of billions of dollars every year. Fraudulent insurance claims (for example through upcoding, billing for services that were never provided, or inflating claim amounts) are hard to detect because they often look very similar to legitimate claims and are hidden among millions of routine submissions. Manual auditing can only review a small fraction of all claims, which makes data driven approaches attractive: machine learning can surface subtle behavioural and financial patterns that are difficult to catch by fixed rules, and it can help an insurer focus limited investigation resources on the claims that matter most.
Interesting questions. This setting raises several questions that we found interesting from a business analytics perspective: Are there natural groups of claims or providers that behave differently, and do some of these segments carry a higher fraud risk? Can the financial outcome of a claim (such as the claim amount) be predicted from its characteristics? And can fraudulent claims be flagged automatically, ideally early enough to act before the claim is paid out? These three questions map onto the three analytical methods we apply in this project.
Dataset. We use the Healthcare Fraud Detection Dataset published on Kaggle by Nudrat Abbas (see Related Work, reference [1]). In contrast to the well known provider level Medicare datasets, this dataset is organized at the level of the individual claim: each row is a single medical insurance claim, and the binary target variable Is_Fraud indicates whether that claim is fraudulent (1) or legitimate (0). The data contains a realistic mix of variable types, including patient demographics (age, gender, state), clinical coding (ICD-10 diagnosis codes and CPT procedure codes), financial information (claimed and approved amounts), provider information (specialty, monthly claim volume), temporal information (submission date and the delay between service and submission), and the claim status. The original public version contains 10,000 claims. For this project we work with an extended and cleaned version prepared by our team during the data preparation step, comprising roughly 18,800 claims with a fraud rate of about 8.7 percent. The prepared data file is included with this submission so that the notebook can be run from top to bottom and reproduced. A full list of variables and their meaning is given in the data preparation section (Step 2).
Research questions tackled. Following the data analytics lifecycle, each team member applies one method to this dataset:
Clustering. (Segmentation). How can healthcare claims be segmented into meaningful groups based on cost structure, patient characteristics, provider behavior, and treatment patterns?
Regression. To what extent can the approved reimbursement amount of a healthcare claim be predicted from claim, patient, provider, and treatment-related characteristics by applying and comparing multiple regression models?
Classification. Can a claim be predicted to be fraudulent from its features, and (as detailed in Related Work) can this be done honestly at submission time, before the insurer adjudicates the claim, rather than relying on information that only becomes available after adjudication?
2. Related Work¶
The Healthcare Fraud Detection Dataset [1] is recent (published on Kaggle in 2026), so there is no peer reviewed work on this exact dataset and only a handful of community notebooks. Most established research instead uses the older provider level Kaggle dataset [2], which labels fraud per provider rather than per claim and is therefore structurally different from ours. We reviewed the most upvoted notebooks on our dataset and, more briefly, the academic literature on the three methods we apply.
2.1 Notebooks on the same dataset. Of the eleven published notebooks, we examined the three most upvoted ones. robiulhasanjisan [3] is methodologically the cleanest: it splits before preprocessing, applies SMOTE only on the training set, and uses SHAP for explanation (Random Forest and XGBoost, ROC-AUC about 0.999). The dataset author's own notebook [4] compares Logistic Regression, Random Forest, and XGBoost (all about 0.999) and engineers a claim to approval ratio and a provider level mean of the target. artheon [5] uses a soft voting ensemble of XGBoost, LightGBM, and CatBoost with class weighting (ROC-AUC about 0.999, fraud F1 about 0.95).
What we observed is that all three reach near perfect scores (ROC-AUC between roughly 0.998 and 0.9997), which is implausible for real fraud detection and points to target leakage. The dominant predictors are features that are known only after a claim has been adjudicated, namely the approved amount, the claim status, and ratios derived from them (in [3] the claim to approval ratio alone has a feature importance of about 0.79). None of the notebooks flag this. We therefore build a leakage aware model that uses only information available at the time a claim is submitted.
2.2 Academic literature. For classification, supervised studies on health insurance claims use logistic regression and tree based models and find financial variables to be the strongest fraud predictors [6], with gradient boosting on similarly imbalanced data reaching an F1 of about 0.85 while reducing false positives [7]. For clustering, fraud is often framed as unsupervised anomaly detection, for example k-means peer grouping of providers [9], association rule mining combined with isolation forest and related detectors [8], and broader surveys of unsupervised methods [10]. For regression, healthcare cost and claim amount studies compare linear models against tree ensembles [11], [13] and show that penalized regression (lasso) improves on plain OLS [12]. The cross cutting challenge throughout is class imbalance, for which SMOTE [14] and ROC based evaluation are the standard references.
2.3 How our work differs. We apply all three methods to the same claim level dataset and, throughout, prioritize avoiding the leakage that inflates the existing notebooks. Each team member owns one method.
Clustering. Our approach is most similar to Massi et al. [9], who cluster hospitals by billing behavior to identify outliers, and De Meulemeester et al. [10], who combine unsupervised anomaly detections with SHAP explanations to flag practitioners with atypical resource use. While both studies aggregate individual claims to analyze data at the provider level, our work differs in two main ways. First, we cluster individual claims instead of providers, which means each observation already carries the full mix of cost, patient, and treatment information without aggregation. Second, rather than using clustering just for outlier detection, we use it to group claims into interpretable segments, and then validate the segments post hoc against the fraud label, similar to the pattern-matching in [8]. Following the insights from [9] and [10] that processing matters more than algorithm choice, we purposefulle use standard K-Means to keept the method simple and interpretable.
Regression. The most directly relevant studies are Langenberger et al. [11] and the claim amount prediction study [13], both of which benchmark linear models against tree-based ensembles for predicting healthcare costs from patient and provider features, as well as Kan et al. [12], which demonstrates that penalized regression such as lasso improves on plain OLS in healthcare cost prediction. These studies treat predictive accuracy of the cost or claim amount as their main objective. This work takes a different angle in two respects. Rather than aiming for maximum predictive performance, regression is used here as a diagnostic tool by training each model both with and without the claimed amount as a feature, it becomes possible to isolate how much of the approved amount stems from the claim itself versus from patient, provider, and treatment characteristics. Picking up on the observation in [12] that penalized regression supports interpretability, lasso is employed not purely for performance but as a feature-selection mechanism to verify which predictors truly carry weight.
Classification. Building on the train/test discipline and SHAP explanation of [3] and the ensembling of [5], but avoiding their leakage, we treat classification as a leakage aware task defined at submission time (see the classification research question in Section 1). Our contribution is to compare a naive full feature model against a pre adjudication model under the same evaluation protocol, and to interpret which legitimate features (above all the time between service and submission) carry the signal, rather than reproducing inflated post adjudication scores.
References
[1] N. Abbas, Healthcare Fraud Detection Dataset (medical insurance claims with ICD-10, CPT codes), Kaggle. https://www.kaggle.com/datasets/nudratabbas/healthcare-fraud-detection-dataset
[2] Rohit Anand Gupta, Healthcare Provider Fraud Detection Analysis, Kaggle. https://www.kaggle.com/datasets/rohitrox/healthcare-provider-fraud-detection-analysis
[3] R. H. Jisan, Fraud Detection Healthcare, Kaggle notebook. https://www.kaggle.com/code/robiulhasanjisan/fraud-detection-healthcare
[4] N. Abbas, Detecting Revenue Loss Before Claims Are Paid, Kaggle notebook. https://www.kaggle.com/code/nudratabbas/detecting-revenue-loss-before-claims-are-paid
[5] artheon, Detecting Medical Fraud: A Trinity Ensemble, Kaggle notebook. https://www.kaggle.com/code/artheon/detecting-medical-fraud-a-trinity-ensemble
[6] "Classification of Health Insurance Fraud Risk with Machine Learning," IEEE Conference Publication, 2024. https://ieeexplore.ieee.org/document/10699052
[7] R. Y. Gupta, S. S. Mudigonda, P. K. Baruah, and P. K. Kandala, "Markov model with machine learning integration for fraud detection in health insurance," arXiv:2102.10978, 2021. https://arxiv.org/abs/2102.10978
[8] Z. Hamid, F. Khalique, S. Mahmood, A. Daud, A. Bukhari, and B. Alshemaimri, "Healthcare insurance fraud detection using data mining," BMC Medical Informatics and Decision Making, 2024. doi:10.1186/s12911-024-02512-4
[9] "Data mining application to healthcare fraud detection: a two-step unsupervised clustering method," BMC Medical Informatics and Decision Making, 2020. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7362640/
[10] "Explainable unsupervised anomaly detection for healthcare insurance data," 2024. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11720628/
[11] B. Langenberger, T. Schulte, and O. Groene, "The application of machine learning to predict high-cost patients using healthcare claims data," PLOS One, 2023. doi:10.1371/journal.pone.0279540
[12] H. J. Kan et al., "A comparison of standard and penalized linear regression models in predicting health care costs in older adults," PLOS One, 2019. doi:10.1371/journal.pone.0213258
[13] "Predicting Health Insurance Claim Amount through Machine Learning Algorithms," IEEE Conference Publication, 2024. https://ieeexplore.ieee.org/document/10625132
[14] N. V. Chawla, K. W. Bowyer, L. O. Hall, and W. P. Kegelmeyer, "SMOTE: Synthetic Minority Over-sampling Technique," Journal of Artificial Intelligence Research, vol. 16, pp. 321-357, 2002. doi:10.1613/jair.953
Step 1: Load Data¶
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
df = pd.read_csv(
"healthcare_fraud_detection.csv",
parse_dates=["Claim_Submission_Date"] # convert date pandas object
)
print(f"Dataset shape: {df.shape}")
df.head()
Dataset shape: (18764, 20)
| Provider_ID | Claim_ID | Patient_Age | Patient_Gender | Diagnosis_Code | Procedure_Code | Claim_Amount | Approved_Amount | Insurance_Type | Claim_Submission_Date | Days_Between_Service_and_Claim | Number_of_Claims_Per_Provider_Monthly | Provider_Specialty | Patient_State | Claim_Status | Is_Fraud | Length_of_Stay | Visit_Type | Chronic_Condition_Flag | Prior_Visits_12m | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | P0052 | C0000000 | 37 | Male | I25.10 | 36415 | 443.51 | 393.16 | Medicaid | 2024-09-01 | 13 | 70 | Cardiology | NY | Approved | 0 | 0 | Outpatient | 1 | 2.0 |
| 1 | P0121 | C0000001 | 21 | Female | E11.9 | 99213 | 467.50 | 461.33 | Self-Pay | 2022-09-05 | 5 | 62 | General Practice | IL | Pending | 0 | 5 | Inpatient | 1 | 2.0 |
| 2 | P0140 | C0000002 | 78 | Female | J06.9 | 93000 | 591.69 | 530.06 | Medicaid | 2022-04-11 | 29 | 60 | Cardiology | IL | Pending | 0 | 5 | Inpatient | 1 | 3.0 |
| 3 | P0202 | C0000003 | 65 | Male | I10 | 93000 | 235.15 | 189.11 | Private | 2023-10-11 | 22 | 70 | General Practice | TX | Approved | 0 | 0 | Emergency | 0 | 5.0 |
| 4 | P0135 | C0000004 | 36 | Male | M54.5 | 85025 | 487.96 | 369.91 | Private | 2023-09-05 | 21 | 67 | Pulmonology | PA | Approved | 0 | 5 | Inpatient | 0 | 4.0 |
Upon loading the dataset we can observe that it consists of 18,764 entries and 20 attributes. The dataset contains a mix of different data types, including numerical, categorical, datetime and binary variables, covering patient demographics, claim details, provider information and the binary target variable Is_Fraud.
Step 2: Prepare and Investigate the Data¶
# Overview of all variables
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 18764 entries, 0 to 18763 Data columns (total 20 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Provider_ID 18764 non-null object 1 Claim_ID 18764 non-null object 2 Patient_Age 18764 non-null int64 3 Patient_Gender 18764 non-null object 4 Diagnosis_Code 18764 non-null object 5 Procedure_Code 18764 non-null int64 6 Claim_Amount 18764 non-null float64 7 Approved_Amount 18764 non-null float64 8 Insurance_Type 18164 non-null object 9 Claim_Submission_Date 18764 non-null datetime64[ns] 10 Days_Between_Service_and_Claim 18764 non-null int64 11 Number_of_Claims_Per_Provider_Monthly 18764 non-null int64 12 Provider_Specialty 18101 non-null object 13 Patient_State 18764 non-null object 14 Claim_Status 18764 non-null object 15 Is_Fraud 18764 non-null int64 16 Length_of_Stay 18764 non-null int64 17 Visit_Type 18764 non-null object 18 Chronic_Condition_Flag 18764 non-null int64 19 Prior_Visits_12m 17936 non-null float64 dtypes: datetime64[ns](1), float64(3), int64(7), object(9) memory usage: 2.9+ MB
The dataset includes a mix of different data types, namely 7 numerical variables such as Patient_Age, Claim_Amount and Length_of_Stay, 9 categorical variables such as Patient_Gender, Insurance_Type and Claim_Status, 1 datetime variable (Claim_Submission_Date) and 1 binary target variable (Is_Fraud) indicating whether a claim is fraudulent (1) or not (0). A first notable observation is that not all columns are complete. Three attributes contain missing values: Insurance_Type with 600, Provider_Specialty with 663 and Prior_Visits_12m with 828 missing values. All remaining 17 columns are fully populated, and these missing values will need to be addressed in the data preparation step.
To get a better understanding of the data, let us take a closer look at the distribution of the key variables in the dataset.
import matplotlib.pyplot as plt
import seaborn as sns
numeric_cols = ["Claim_Amount", "Approved_Amount", "Days_Between_Service_and_Claim", #distrubtion of numeric vari.
"Number_of_Claims_Per_Provider_Monthly", "Length_of_Stay",
"Patient_Age", "Prior_Visits_12m"]
fig, axes = plt.subplots(3, 3, figsize=(15, 12))
axes = axes.flatten()
for i, col in enumerate(numeric_cols):
sns.histplot(data=df, x=col, ax=axes[i], color="steelblue", kde=True)
axes[i].set_title(col)
axes[i].set_xlabel("")
plt.suptitle("Distribution of Key Numerical Variables", fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
The histograms reveal several interesting patterns across the key numerical variables. Both Claim_Amount and Approved_Amount show a strongly right-skewed distribution, with the majority of claims concentrated at lower values and a long tail towards higher amounts. This is expected in a healthcare fraud context as most claims are of moderate value while a small number of fraudulent or unusual claims tend to involve significantly higher amounts. Days_Between_Service_and_Claim shows a relatively uniform distribution across the range of 0 to 29 days, with a slight peak at both extremes, suggesting that claims are submitted either very quickly after the service or close to the maximum allowed time window. Number_of_Claims_Per_Provider_Monthly follows an approximately normal distribution centered around 65-70 claims per month, which indicates a relatively consistent workload across providers with some variation. Length_of_Stay shows a uniform discrete distribution across 0 to 5 days, meaning all lengths of stay are roughly equally common in the dataset. Patient_Age follows a roughly normal distribution centered around 45-55 years, which is typical for a healthcare dataset where middle-aged and older patients tend to be more frequent users of healthcare services. Prior_Visits_12m is right-skewed with a large spike at 0, which reflects our imputation decision to fill missing values with 0, representing patients with no recorded prior visits.
Now lets compare fraud vs non fraud data.
fig, axes = plt.subplots(3, 3, figsize=(15, 12))
axes = axes.flatten()
for i, col in enumerate(numeric_cols):
sns.boxplot(data=df, x="Is_Fraud", y=col, ax=axes[i],
hue="Is_Fraud", palette={0: "steelblue", 1: "tomato"},
legend=False)
axes[i].set_title(col)
axes[i].set_xlabel("Is_Fraud (0=No, 1=Yes)")
plt.suptitle("Fraud vs. Non-Fraud Comparison", fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
The Fraud vs. Non-Fraud comparison reveals several important patterns. The most striking difference can be observed in Days_Between_Service_and_Claim, where fraudulent claims show a significantly lower median and a much narrower distribution compared to non-fraudulent claims. This suggests that fraudulent claims tend to be submitted very quickly after the service, which is a strong indicator of fraud. Claim_Amount and Approved_Amount show that fraudulent claims have more extreme outliers towards higher values, confirming that unusually high claim amounts are associated with fraudulent behavior. Number_of_Claims_Per_Provider_Monthly shows a slightly higher median for fraudulent claims, indicating that providers with higher claim volumes are somewhat more likely to be associated with fraud. Length_of_Stay, Patient_Age and Prior_Visits_12m show very similar distributions between fraudulent and non-fraudulent claims, suggesting that these variables alone are not strong indicators of fraud. This is consistent with the low correlations we observed earlier. Overall the EDA confirms that Days_Between_Service_and_Claim and Claim_Amount are the most promising features for fraud detection, while variables such as Length_of_Stay and Patient_Age are likely to contribute less to the predictive models.
Lets also have a look on the correlation matrix.
plt.figure(figsize=(15, 12))
corr_matrix = df[numeric_cols + ["Is_Fraud"]].corr()
sns.heatmap(corr_matrix, annot=True, fmt=".2f", cmap="coolwarm", center=0)
plt.title("Correlation Matrix")
plt.tight_layout()
plt.show()
The correlation matrix highlights three key variables that show a meaningful relationship with Is_Fraud. Days_Between_Service_and_Claim shows the strongest correlation at -0.34, indicating that fraudulent claims are submitted much faster after the service. Number_of_Claims_Per_Provider_Monthly shows a small positive correlation of 0.10, suggesting that providers with higher claim volumes are slightly more likely to be associated with fraud. Claim_Amount shows a weak positive correlation of 0.09, meaning higher claim amounts are marginally more common in fraudulent cases. These three variables will therefore be the most important predictors in our subsequent fraud detection models.
Missing values¶
Ok now lets observe the missing data that we could see earlier
# Missing values overview
missing = df.isnull().sum()
missing_pct = (missing / len(df) * 100).round(2)
missing_df = pd.DataFrame({
"missing_count": missing,
"missing_pct": missing_pct
}).query("missing_count > 0")
print(missing_df)
missing_count missing_pct Insurance_Type 600 3.20 Provider_Specialty 663 3.53 Prior_Visits_12m 828 4.41
# Check whether the missing values appear in the same rows
for col in ["Insurance_Type", "Provider_Specialty", "Prior_Visits_12m"]:
mask = df[col].isnull()
print(f"Columns where {col} is missing: {mask.sum()}")
for other_col in ["Insurance_Type", "Provider_Specialty", "Prior_Visits_12m"]:
if other_col != col:
pct = (df[mask][other_col].isnull().sum() / mask.sum() * 100).round(2)
print(f" of which is also missing {other_col}: {pct}%")
print()
Columns where Insurance_Type is missing: 600 of which is also missing Provider_Specialty: 3.83% of which is also missing Prior_Visits_12m: 3.5% Columns where Provider_Specialty is missing: 663 of which is also missing Insurance_Type: 3.47% of which is also missing Prior_Visits_12m: 4.37% Columns where Prior_Visits_12m is missing: 828 of which is also missing Insurance_Type: 2.54% of which is also missing Provider_Specialty: 3.5%
# # Do the missing values correlate with any other num. variable
for col in ["Insurance_Type", "Provider_Specialty", "Prior_Visits_12m"]:
df[f"{col}_missing"] = df[col].isnull().astype(int)
missing_cols = ["Insurance_Type_missing", "Provider_Specialty_missing", "Prior_Visits_12m_missing"]
numeric_cols = [
"Patient_Age", "Procedure_Code", "Claim_Amount", "Approved_Amount",
"Days_Between_Service_and_Claim", "Number_of_Claims_Per_Provider_Monthly",
"Is_Fraud", "Length_of_Stay", "Chronic_Condition_Flag", "Prior_Visits_12m"
]
corr = df[missing_cols + numeric_cols].corr()[missing_cols].drop(missing_cols)
print(corr.round(3))
df.drop(columns=missing_cols, inplace=True)
Insurance_Type_missing \
Patient_Age 0.004
Procedure_Code -0.001
Claim_Amount -0.007
Approved_Amount -0.007
Days_Between_Service_and_Claim 0.009
Number_of_Claims_Per_Provider_Monthly -0.006
Is_Fraud -0.007
Length_of_Stay -0.010
Chronic_Condition_Flag 0.003
Prior_Visits_12m 0.001
Provider_Specialty_missing \
Patient_Age 0.001
Procedure_Code -0.005
Claim_Amount -0.000
Approved_Amount -0.001
Days_Between_Service_and_Claim 0.004
Number_of_Claims_Per_Provider_Monthly 0.001
Is_Fraud 0.003
Length_of_Stay 0.007
Chronic_Condition_Flag 0.002
Prior_Visits_12m -0.005
Prior_Visits_12m_missing
Patient_Age -0.009
Procedure_Code -0.003
Claim_Amount 0.005
Approved_Amount 0.008
Days_Between_Service_and_Claim 0.001
Number_of_Claims_Per_Provider_Monthly -0.001
Is_Fraud 0.005
Length_of_Stay 0.005
Chronic_Condition_Flag -0.005
Prior_Visits_12m NaN
Upon closer inspection of the three columns containing missing values, we can observe that the missing entries do not follow any systematic pattern. As shown above, the correlations between the missing value indicators and all other numerical variables are extremely close to zero (all below ±0.01), indicating no meaningful relationship. Furthermore, only a very small percentage of rows where one variable is missing also have missing values in the other two columns (below 5% in all cases), confirming that the missing values are independent of each other. We therefore conclude that the missing values occur completely at random (MCAR — Missing Completely At Random). The most likely explanation is that the data was simply not recorded at the time of data collection, for example due to incomplete documentation or administrative gaps during the claims process. This assumption allows us to proceed with standard imputation techniques without introducing significant bias into the dataset.
Imputation of missing values¶
# Insurance_Type & Provider_Specialty → "Unknown" as own category
df["Insurance_Type"] = df["Insurance_Type"].fillna("Unknown")
df["Provider_Specialty"] = df["Provider_Specialty"].fillna("Unknown")
# Prior_Visits_12m → 0 (No previous visits are known)
df["Prior_Visits_12m"] = df["Prior_Visits_12m"].fillna(0)
# Check
print(df[["Insurance_Type", "Provider_Specialty", "Prior_Visits_12m"]].isnull().sum())
Insurance_Type 0 Provider_Specialty 0 Prior_Visits_12m 0 dtype: int64
For the missing values we decided on the following approach. For Insurance_Type and Provider_Specialty we impute the missing values with the category "Unknown", as we cannot make any reasonable assumption about which category should be assigned. This way we preserve the information that the value was missing rather than artificially introducing a category. For Prior_Visits_12m we impute with 0, as the absence of this information can be interpreted as no prior visits being recorded, which is analogous to a real-world medical setting where if no prior visits are documented the assumption is that there were none.
Outlier detection and handling¶
# Outlier Detection — Boxplots
import matplotlib.pyplot as plt
import seaborn as sns
numeric_cols = ["Claim_Amount", "Approved_Amount", "Days_Between_Service_and_Claim",
"Number_of_Claims_Per_Provider_Monthly", "Length_of_Stay",
"Patient_Age", "Prior_Visits_12m"]
fig, axes = plt.subplots(3, 3, figsize=(15, 12))
axes = axes.flatten()
for i, col in enumerate(numeric_cols):
sns.boxplot(data=df, x=col, ax=axes[i], color="steelblue")
axes[i].set_title(col)
plt.tight_layout()
plt.show()
The boxplots clearly show that Claim_Amount and Approved_Amount contain extreme outliers that lie far beyond the typical range of the data. To reduce the influence of these extreme values on our subsequent machine learning models without losing any observations, we apply capping at the 99th percentile. This means that all values above the 99th percentile threshold are set to that threshold value, preserving the information that these claims are unusually high while preventing them from distorting our models.
# Capping to 95. Percentile
df_95 = df.copy()
for col in ["Claim_Amount", "Approved_Amount"]:
p95 = df[col].quantile(0.95)
df_95[col] = df_95[col].clip(upper=p95)
# Capping to 99. Percentile
df_99 = df.copy()
for col in ["Claim_Amount", "Approved_Amount"]:
p99 = df[col].quantile(0.99)
df_99[col] = df_99[col].clip(upper=p99)
# Boxplots comparison
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
for i, col in enumerate(["Claim_Amount", "Approved_Amount"]):
sns.boxplot(data=df, x=col, ax=axes[i][0], color="steelblue")
axes[i][0].set_title(f"{col} — Original")
sns.boxplot(data=df_95, x=col, ax=axes[i][1], color="orange")
axes[i][1].set_title(f"{col} — 95. Percentile")
sns.boxplot(data=df_99, x=col, ax=axes[i][2], color="green")
axes[i][2].set_title(f"{col} — 99. Percentile")
plt.tight_layout()
plt.show()
As the comparison shows, capping at the 99th percentile preserves more of the natural distribution while still removing the most extreme outliers. We therefore decide to cap both Claim_Amount and Approved_Amount at the 99th percentile, which affects 188 observations in each column
for col in ["Claim_Amount", "Approved_Amount"]:
p99 = df[col].quantile(0.99)
print(f"{col} — Values above the 99th percentile ({p99:.2f}): {(df[col] > p99).sum()}")
Claim_Amount — Values above the 99th percentile (1884.22): 188 Approved_Amount — Values above the 99th percentile (1520.75): 188
for col in ["Claim_Amount", "Approved_Amount"]:
p99 = df[col].quantile(0.99)
df[col] = df[col].clip(upper=p99)
# Kontrollcheck
print(f"Claim_Amount Max: {df['Claim_Amount'].max():.2f}")
print(f"Approved_Amount Max: {df['Approved_Amount'].max():.2f}")
Claim_Amount Max: 1884.22 Approved_Amount Max: 1520.75
During the outlier detection step we visualized the distribution of all numerical variables using boxplots. The analysis revealed that the columns Claim_Amount and Approved_Amount contain a number of extreme values that lie far beyond the typical range of the data. These outliers were identified as values exceeding the 99th percentile, accounting for approximately 1% of the data in each column. For all other numerical variables such as Days_Between_Service_and_Claim, Length_of_Stay, Patient_Age, Prior_Visits_12m and Number_of_Claims_Per_Provider_Monthly no extreme outliers were detected and the distributions appear reasonable and within expected ranges. We decided to apply capping at the 99th percentile, meaning all values above this threshold are set to the 99th percentile value rather than being removed. We chose this approach as extreme values can distort our machine learning models and negatively impact the quality of our results.
Feature Engineering¶
In the feature engineering step we derived 10 additional variables from the existing dataset, bringing the total number of columns from 20 to 30. The new features were created with the goal of providing the machine learning models with more meaningful and interpretable information. From the Claim_Submission_Date we extracted the month, year and weekday of the claim submission, as fraudulent claims may follow seasonal patterns or occur more frequently on specific days of the week. From the financial columns we derived the Claim_Ratio, which represents the proportion of the claimed amount that was approved, as well as the Amount_Difference, which captures the absolute gap between the two. Both variables are particularly relevant in a fraud detection context since fraudulent claims tend to have unusually high claim amounts relative to what gets approved. Regarding patient characteristics we created a binary gender variable, an age group variable splitting patients into three groups (20-35, 35-50 and 50+) and a binary Is_Senior variable for patients above 65, as age and gender may influence the likelihood of fraudulent claims. Additionally we added a Has_Prior_Visits variable indicating whether a patient had any prior visits in the last 12 months, as patients with no prior history may be more suspicious. Finally we created a High_Volume_Provider variable flagging providers who submit more than 80 claims per month, as an unusually high number of claims from a single provider is a well known indicator of potential fraud.
# -----------------------------------------------------------
# Feature Engineering
# -----------------------------------------------------------
# From Claim_Submission_Date
df["Claim_Month"] = df["Claim_Submission_Date"].dt.month
df["Claim_Year"] = df["Claim_Submission_Date"].dt.year
df["Claim_Weekday"] = df["Claim_Submission_Date"].dt.dayofweek
# From Claim_Amount / Approved_Amount
df["Claim_Ratio"] = (df["Claim_Amount"] / df["Approved_Amount"]).round(2)
df["Amount_Difference"] = (df["Claim_Amount"] - df["Approved_Amount"]).round(2)
# Gender Binary
df["Patient_Gender_Binary"] = (df["Patient_Gender"] == "Female").astype(int)
# Age Groups
df["Age_Group"] = pd.cut(
df["Patient_Age"],
bins=[0, 35, 50, 100],
labels=["20-35", "35-50", "50+"]
)
# Provider related
df["High_Volume_Provider"] = (df["Number_of_Claims_Per_Provider_Monthly"] > 80).astype(int)
# Patient related
df["Is_Senior"] = (df["Patient_Age"] > 65).astype(int)
df["Has_Prior_Visits"] = (df["Prior_Visits_12m"] > 0).astype(int)
print(f"New columns: {df.shape[1]} (before 20)")
print(df[["Claim_Ratio", "Amount_Difference", "Patient_Gender_Binary",
"Age_Group", "High_Volume_Provider", "Is_Senior", "Has_Prior_Visits",
"Claim_Month", "Claim_Year", "Claim_Weekday"]].head(5))
New columns: 30 (before 20) Claim_Ratio Amount_Difference Patient_Gender_Binary Age_Group \ 0 1.13 50.35 0 35-50 1 1.01 6.17 1 20-35 2 1.12 61.63 1 50+ 3 1.24 46.04 0 50+ 4 1.32 118.05 0 35-50 High_Volume_Provider Is_Senior Has_Prior_Visits Claim_Month Claim_Year \ 0 0 0 1 9 2024 1 0 0 1 9 2022 2 0 1 1 4 2022 3 0 0 1 10 2023 4 0 0 1 9 2023 Claim_Weekday 0 6 1 0 2 0 3 2 4 1
Scaling of variables¶
In the next step we apply feature scaling to all numerical variables using the StandardScaler. Scaling is an important preprocessing step as many machine learning algorithms are sensitive to the magnitude of input variables. Without scaling, variables with larger ranges such as Claim_Amount or Approved_Amount would dominate the model compared to variables with smaller ranges such as Length_of_Stay or Claim_Month. The StandardScaler transforms each numerical variable to have a mean of 0 and a standard deviation of 1, ensuring that all variables contribute equally to the model. Since we already capped the outliers in the previous step the StandardScaler can be applied without being distorted by extreme values.
from sklearn.preprocessing import StandardScaler
# save org, data
df_original = df.copy()
# Claim_Submission_Date to Unix Timestamp
df_scaled = df.copy()
df_scaled["Claim_Submission_Date"] = df_scaled["Claim_Submission_Date"].astype(int) // 10**9
numeric_cols = ["Patient_Age", "Claim_Amount", "Approved_Amount",
"Days_Between_Service_and_Claim", "Number_of_Claims_Per_Provider_Monthly",
"Length_of_Stay", "Prior_Visits_12m", "Claim_Ratio",
"Amount_Difference", "Claim_Month", "Claim_Year", "Claim_Weekday",
"Claim_Submission_Date"]
scaler = StandardScaler()
df_scaled[numeric_cols] = scaler.fit_transform(df_scaled[numeric_cols])
print(f"Shape: {df_scaled.shape}")
print(df_scaled[numeric_cols].describe().round(2))
Shape: (18764, 30)
Patient_Age Claim_Amount Approved_Amount \
count 18764.00 18764.00 18764.00
mean 0.00 -0.00 0.00
std 1.00 1.00 1.00
min -2.71 -1.40 -1.40
25% -0.70 -0.75 -0.75
50% 0.03 -0.19 -0.20
75% 0.70 0.56 0.56
max 2.54 3.49 3.30
Days_Between_Service_and_Claim Number_of_Claims_Per_Provider_Monthly \
count 18764.00 18764.00
mean 0.00 0.00
std 1.00 1.00
min -1.69 -1.79
25% -0.86 -0.64
50% -0.03 -0.17
75% 0.81 0.44
max 1.76 5.11
Length_of_Stay Prior_Visits_12m Claim_Ratio Amount_Difference \
count 18764.00 18764.00 18764.00 18764.00
mean 0.00 0.00 -0.00 0.00
std 1.00 1.00 1.00 1.00
min -1.29 -1.64 -1.71 -3.26
25% -0.70 -0.63 -0.56 -0.60
50% -0.12 0.06 -0.15 -0.29
75% 0.47 0.62 0.22 0.25
max 1.64 5.15 8.71 9.77
Claim_Month Claim_Year Claim_Weekday Claim_Submission_Date
count 18764.00 18764.00 18764.00 18764.00
mean -0.00 0.00 -0.00 0.00
std 1.00 1.00 1.00 1.00
min -1.59 -1.34 -1.51 -1.73
25% -0.72 -0.46 -1.01 -0.86
50% 0.15 0.41 -0.01 -0.02
75% 0.72 1.29 0.99 0.87
max 1.59 2.16 1.49 1.76
Correlation Matrix after the preprocessing¶
plt.figure(figsize=(15, 12))
corr_matrix = df_scaled[numeric_cols + ["Is_Fraud"]].corr()
sns.heatmap(corr_matrix, annot=True, fmt=".2f", cmap="coolwarm", center=0)
plt.title("Correlation Matrix")
plt.tight_layout()
plt.show()
The correlation matrix reveals several important relationships in the dataset. The most notable finding is the strong positive correlation between Claim_Amount and Approved_Amount at 0.99, confirming that these two variables are almost perfectly linearly related. Claim_Ratio shows the strongest correlation with Is_Fraud at 0.66, followed by Days_Between_Service_and_Claim at -0.34 and Amount_Difference at 0.23, identifying these as the most important fraud indicators in the dataset. However it is important to note that Claim_Ratio and Amount_Difference are derived from Approved_Amount, which is only known after a claim has been adjudicated — meaning these variables constitute target leakage and cannot be used as features in a honest fraud detection model at submission time. This is exactly the issue identified in the Related Work section, where existing notebooks achieve near-perfect scores by relying on post-adjudication information. Additionally Claim_Year and Claim_Submission_Date are almost perfectly correlated at 0.97, meaning they carry essentially the same information. All other variable pairs show correlations close to zero, indicating that the remaining features are largely independent of each other.
Encoding of categorical variables¶
In the final preprocessing step we apply One-Hot Encoding to all remaining categorical variables. One-Hot Encoding converts each categorical variable into a set of binary columns, one for each unique category, where a value of 1 indicates that the observation belongs to that category and 0 otherwise. This is necessary because most machine learning algorithms cannot work directly with text-based categorical values and require numerical input. We exclude Provider_ID and Claim_ID from the encoding as these are unique identifiers that carry no predictive information for the models. The remaining 6 categorical variables "Diagnosis_Code, Insurance_Type, Provider_Specialty, Patient_State, Claim_Status and Visit_Type", are encoded using pd.get_dummies(), expanding the dataset from 30 to 60 columns in total. The encoded dataset is stored in df_scaled and will be used for all subsequent modeling steps.
cat_cols = df.select_dtypes(include="object").columns.tolist()
print(cat_cols)
print(f"\nCount: {len(cat_cols)}")
['Provider_ID', 'Claim_ID', 'Patient_Gender', 'Diagnosis_Code', 'Insurance_Type', 'Provider_Specialty', 'Patient_State', 'Claim_Status', 'Visit_Type'] Count: 9
# Encoding for cate. variables
# Provider_ID und Claim_ID are left out they are just IDs
cols_to_encode = ["Diagnosis_Code", "Insurance_Type", "Provider_Specialty",
"Patient_State", "Claim_Status", "Visit_Type"]
# One-Hot Encoding
df_scaled = pd.get_dummies(df_scaled, columns=cols_to_encode, drop_first=False)
print(f"Shape after Encoding: {df_scaled.shape}")
print(f"New Columns: {df_scaled.columns.tolist()}")
Shape after Encoding: (18764, 60) New Columns: ['Provider_ID', 'Claim_ID', 'Patient_Age', 'Patient_Gender', 'Procedure_Code', 'Claim_Amount', 'Approved_Amount', 'Claim_Submission_Date', 'Days_Between_Service_and_Claim', 'Number_of_Claims_Per_Provider_Monthly', 'Is_Fraud', 'Length_of_Stay', 'Chronic_Condition_Flag', 'Prior_Visits_12m', 'Claim_Month', 'Claim_Year', 'Claim_Weekday', 'Claim_Ratio', 'Amount_Difference', 'Patient_Gender_Binary', 'Age_Group', 'High_Volume_Provider', 'Is_Senior', 'Has_Prior_Visits', 'Diagnosis_Code_E11.9', 'Diagnosis_Code_E78.5', 'Diagnosis_Code_F41.9', 'Diagnosis_Code_I10', 'Diagnosis_Code_I25.10', 'Diagnosis_Code_J06.9', 'Diagnosis_Code_J18.9', 'Diagnosis_Code_K21.9', 'Diagnosis_Code_M54.5', 'Diagnosis_Code_N39.0', 'Insurance_Type_Medicaid', 'Insurance_Type_Medicare', 'Insurance_Type_Private', 'Insurance_Type_Self-Pay', 'Insurance_Type_Unknown', 'Provider_Specialty_Cardiology', 'Provider_Specialty_General Practice', 'Provider_Specialty_Internal Medicine', 'Provider_Specialty_Neurology', 'Provider_Specialty_Orthopedics', 'Provider_Specialty_Pulmonology', 'Provider_Specialty_Unknown', 'Patient_State_CA', 'Patient_State_FL', 'Patient_State_GA', 'Patient_State_IL', 'Patient_State_NY', 'Patient_State_OH', 'Patient_State_PA', 'Patient_State_TX', 'Claim_Status_Approved', 'Claim_Status_Pending', 'Claim_Status_Rejected', 'Visit_Type_Emergency', 'Visit_Type_Inpatient', 'Visit_Type_Outpatient']
Summary of data processing¶
In the data preparation step we performed several important preprocessing tasks to ensure the dataset is ready for subsequent analysis and model training. We started by addressing the missing values in the dataset. Three columns contained missing values: Insurance_Type with 600, Provider_Specialty with 663 and Prior_Visits_12m with 828 missing entries. Through correlation analysis we confirmed that the missing values occur completely at random (MCAR) and are not related to any other variable or to the target variable Is_Fraud. For Insurance_Type and Provider_Specialty we filled the missing values with the category "Unknown" to preserve the information that the value was not recorded. For Prior_Visits_12m we imputed with 0, assuming that no recorded prior visits means no prior visits took place. In the outlier detection step we visualized all numerical variables using boxplots and identified extreme values in Claim_Amount and Approved_Amount. We decided to cap these values at the 99th percentile to reduce the influence of extreme values on our models while keeping all observations in the dataset. In the feature engineering step we created 10 additional variables including financial ratios such as Claim_Ratio and Amount_Difference, date-based features such as Claim_Month, Claim_Year and Claim_Weekday, and patient and provider related binary variables such as Is_Senior, Has_Prior_Visits and High_Volume_Provider. These new features provide the machine learning models with additional meaningful information beyond the raw variables. Finally we applied StandardScaler to all numerical variables and One-Hot Encoding to all categorical variables, resulting in a final dataset of 18,764 observations and 60 features stored in df_scaled, ready for model training.
Step 3 - Clustering team member¶
In this step, we apply clustering to address the research question of how healthcare claims can be segmented into meaningful groups based on cost structure, patient characteristics, provider behavior, and treatment patterns. In line with the case description, clustering is used as an unsupervised method to uncover natural groupings in the data without relying on the fraud label during model training. We use the K-Means algorithm, as it provides an interpretable and straightforward way to identify homogeneous groups of claims based on their feature similarity. The fraud variable Is_Fraud is excluded from the clustering process and only used later for validation.
3.1 Clustering Preparation¶
For the clustering, we selected 11 features covering the cost structure, patient profile, provider behavior, and treatment intensity. This selection ensures consistency with the research question, as it reflects the main drivers of claim characteristics described in the case setting. All variables are standardized using the StandardScaler. This step is necessary because K-Means is distance-based, and without scaling, variables with larger magnitudes would dominate the clustering results.
Clustering Variables
• Cost structure: Claim_Amount, Approved_Amount, Claim_Ratio, Amount_Difference
• Patient profile: Patient_Age, Chronic_Condition_Flag, Prior_Visits_12m
• Provider behavior: Number_of_Claims_Per_Provider_Monthly, Days_Between_Service_and_Claim
• Treatment intensity: Length_of_Stay
• Visit type (encoded): Visit_Type (ordinal: Outpatient < Emergency < Inpatient)
Visit_Type is ordinally encoded (Outpatient < Emergency < Inpatient) to capture treatment intensity naturally.
# Derived columns needed for clustering
df["Claim_Ratio"] = (df["Claim_Amount"] / df["Approved_Amount"]).round(2)
df["Amount_Difference"] = (df["Claim_Amount"] - df["Approved_Amount"]).round(2)
df["Prior_Visits_12m"] = df["Prior_Visits_12m"].fillna(0)
# Ordinal encoding for Visit_Type
visit_order = {"Outpatient": 0, "Emergency": 1, "Inpatient": 2}
df["Visit_Type_Enc"] = df["Visit_Type"].map(visit_order)
# Select Clustering Variables
CLUSTER_FEATURES = [
"Claim_Amount", "Approved_Amount", "Claim_Ratio", "Amount_Difference", "Patient_Age",
"Chronic_Condition_Flag", "Prior_Visits_12m", "Number_of_Claims_Per_Provider_Monthly",
"Days_Between_Service_and_Claim", "Length_of_Stay", "Visit_Type_Enc",
]
C_F = df[CLUSTER_FEATURES].copy()
print(f"Feature matrix: {C_F.shape}")
Feature matrix: (18764, 11)
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import davies_bouldin_score
# Scale Features
scaler = StandardScaler()
CF_scaled = scaler.fit_transform(C_F)
# Determine Number of Clusters
K_RANGE = range(2, 10)
wss_scores = []
db_scores = []
for k in K_RANGE:
km = KMeans(n_clusters=k, n_init=10, random_state=42)
labels = km.fit_predict(CF_scaled)
wss_scores.append(km.inertia_)
db_scores.append(davies_bouldin_score(CF_scaled, labels))
scores_df = pd.DataFrame({"K": list(K_RANGE), "WSS": wss_scores, "DB Score": db_scores})
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sns.lineplot(data=scores_df, x="K", y="WSS", ax=axes[0], marker="o", color="steelblue")
sns.lineplot(data=scores_df, x="K", y="DB Score", ax=axes[1], marker="o", color="red")
axes[0].set_title("Elbow Method")
axes[1].set_title("Davies-Bouldin Score (lower = better)")
plt.tight_layout()
plt.show()
print(scores_df.to_string(index=False))
K WSS DB Score 2 176884.333159 2.223680 3 161941.350569 1.945567 4 148299.047693 1.992443 5 140846.757250 2.152198 6 135213.779462 2.044764 7 130291.215478 1.955394 8 126166.817841 2.082354 9 122615.258212 2.052689
To determine the optimal number of clusters, we evaluate both the Elbow Method and the Davies-Bouldin score. The Elbow Method is used to identify the point at which adding more clusters does not significantly improve the model, while the Davies-Bouldin score measures cluster separation and compactness.
Both methods indicate that four clusters provide a reasonable balance between interpretability and model performance. Therefore, we proceed with k = 4.
# Final clustering with k=4
K_FINAL = 4
kmeans = KMeans(n_clusters=K_FINAL, n_init=10, random_state=42)
df["Cluster"] = kmeans.fit_predict(CF_scaled)
print("Cluster sizes:")
print(df["Cluster"].value_counts().sort_index())
Cluster sizes: Cluster 0 4412 1 9328 2 4406 3 618 Name: count, dtype: int64
# Cluster mean profiles
cluster_profile = df.groupby("Cluster")[CLUSTER_FEATURES].mean().round(2)
print(cluster_profile.T.to_string())
Cluster 0 1 2 3 Claim_Amount 443.08 387.41 1076.13 948.17 Approved_Amount 382.92 335.78 917.84 462.60 Claim_Ratio 1.17 1.17 1.19 2.19 Amount_Difference 60.16 51.63 158.28 485.57 Patient_Age 49.20 49.74 49.42 48.53 Chronic_Condition_Flag 1.00 0.00 0.19 0.29 Prior_Visits_12m 2.91 2.90 2.87 2.98 Number_of_Claims_Per_Provider_Monthly 68.21 68.01 68.60 75.46 Days_Between_Service_and_Claim 14.70 14.68 14.06 4.81 Length_of_Stay 2.20 2.19 2.23 2.17 Visit_Type_Enc 1.01 0.99 0.93 0.98
Applying K-Means with four clusters results in segments of varying sizes, indicating that the dataset contains both common and more specialized claim patterns. The cluster mean profiles show clear differences, particularly in financial variables and submission timing, suggesting that the algorithm successfully captures meaningful structural differences between claims.
3.2 Cluster Interpretation of Numerical Variables¶
To better understand the distribution of key variables across clusters, we use boxplots. Boxplots are particularly useful in this context because they allow us to compare medians, variability, and the presence of outliers across multiple groups simultaneously.
The analysis shows that one cluster clearly stands out with significantly higher claim amounts, higher claim ratios, and a much shorter time between service and claim submission. This indicates a distinct behavioral pattern compared to the other clusters. In contrast, the remaining clusters exhibit more moderate and similar distributions, suggesting more typical claim behavior.
# Boxplots (key variables by cluster)
plot_vars = [
("Claim_Amount", "Claim Amount ($)"),
("Approved_Amount", "Approved Amount ($)"),
("Claim_Ratio", "Claim Ratio"),
("Days_Between_Service_and_Claim", "Days: Service → Claim"),
("Number_of_Claims_Per_Provider_Monthly", "Provider Monthly Claims"),
("Length_of_Stay", "Length of Stay"),
("Patient_Age", "Patient Age"),
("Prior_Visits_12m","Prior Visits (12m)"),
]
fig, axes = plt.subplots(2, 4, figsize=(18, 8))
for ax, (var, label) in zip(axes.flatten(), plot_vars):
sns.boxplot(data=df, x="Cluster", y=var, hue="Cluster", palette="tab10", legend=False, ax=ax)
ax.set_title(label)
ax.set_xlabel("Cluster")
ax.set_ylabel("")
plt.suptitle("Key Feature Distributions by Cluster", fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
3.3 Cluster Interpretation of Categorical Variables¶
In addition to numerical variables, we analyze categorical features across clusters using proportional bar charts. This method is used to understand how different categories are distributed within each cluster and to identify structural differences in claim composition.
The results show that clusters differ not only in financial characteristics but also in terms of treatment type, insurance type, and patient condition. For example, one cluster is strongly associated with chronic conditions, indicating a more stable and likely legitimate patient group.
# Categorical breakdown by cluster
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
for ax, var in zip(axes, ["Visit_Type", "Insurance_Type", "Chronic_Condition_Flag"]):
ct = (df.groupby(["Cluster", var])
.size()
.unstack(fill_value=0)
.apply(lambda r: r / r.sum() * 100, axis=1))
ct.plot(kind="bar", stacked=True, ax=ax, colormap="tab20", width=0.6)
ax.set_title(f"{var} per Cluster (%)")
ax.set_xlabel("Cluster")
ax.set_ylabel("Share (%)")
ax.legend(fontsize=7)
ax.tick_params(axis="x", rotation=0)
plt.tight_layout()
plt.show()
3.4 Fraud Validation¶
After forming the clusters, we evaluate how fraud is distributed across them. This step is crucial to assess whether the clusters identified by the unsupervised model correspond to meaningful differences in fraud risk.
The results show substantial variation in fraud rates across clusters. While some clusters exhibit very low fraud rates, one cluster shows an extremely high fraud rate, indicating that the clustering successfully isolates a high-risk group of claims.
# Fraud validation (only "Is_Fraud" used here)
fraud_by_cluster = (df.groupby("Cluster")["Is_Fraud"]
.agg(Total="count", Fraud_Count="sum")
.assign(Fraud_Rate_Pct=lambda d: (d["Fraud_Count"] / d["Total"] * 100).round(2)))
print(fraud_by_cluster)
fig, ax = plt.subplots(figsize=(6, 4))
sns.barplot(data=fraud_by_cluster.reset_index(),
x="Cluster", y="Fraud_Rate_Pct", hue="Cluster", palette="tab10", legend=False, ax=ax)
ax.set_title("Fraud Rate (%) per Cluster")
ax.set_ylabel("Fraud Rate (%)")
for p in ax.patches:
ax.annotate(f"{p.get_height():.1f}%",
(p.get_x() + p.get_width() / 2, p.get_height()),
ha="center", va="bottom", fontsize=9)
plt.tight_layout()
plt.show()
Total Fraud_Count Fraud_Rate_Pct Cluster 0 4412 193 4.37 1 9328 342 3.67 2 4406 479 10.87 3 618 611 98.87
Summary & Outcome of Clustering¶
The clustering results provide a clear answer to the research question. Healthcare claims can indeed be segmented into meaningful groups based on cost structure, patient characteristics, provider behavior, and treatment patterns. The analysis shows that fraud is not driven by a single variable such as claim amount alone. Instead, it is associated with a combination of factors, particularly high claim ratios, large discrepancies between claimed and approved amounts, and very short submission delays.
Cluster 1 (n = 9,328) represents the largest group of claims and is characterized by relatively low costs and stable patterns, with a low fraud rate of 3.67%, indicating typical legitimate claims. Cluster 0 (n = 4,412) shows a similar cost structure but a slightly higher fraud rate of 4.37%, which is close to the overall dataset average. Cluster 2 (n = 4,406) stands out with noticeably higher claim amounts and a fraud rate of 10.87%, suggesting a segment with elevated risk and more irregular financial patterns.
Cluster 3 (n = 618) forms a small but highly distinct group with an extremely high fraud rate of 98.87%. This cluster is clearly separated from the others and represents a high-risk segment characterized by unusual cost structures and rapid claim submissions.
Overall, clustering proves to be a valuable tool for identifying hidden structures in the data and for segmenting claims into groups with clearly different fraud risk profiles.
Step 4: Regression Modeling team member¶
4.1 Feature Importance — Linear Regression¶
In order to identify which features have the greatest influence on the predicted approved amount, we first train a Linear Regression model and examine its coefficients. The magnitude of each coefficient indicates how strongly the corresponding feature influences the prediction —> larger absolute values indicate more important features. This step helps us understand the underlying relationships in the data before moving on to more complex models.
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
features = [
"Claim_Amount", "Days_Between_Service_and_Claim",
"Patient_Age", "Patient_Gender_Binary", "Chronic_Condition_Flag", "Prior_Visits_12m",
"Number_of_Claims_Per_Provider_Monthly", "Procedure_Code", "Length_of_Stay",
"Insurance_Type_Medicaid", "Insurance_Type_Medicare",
"Insurance_Type_Private", "Insurance_Type_Self-Pay", "Insurance_Type_Unknown",
"Provider_Specialty_Cardiology", "Provider_Specialty_General Practice",
"Provider_Specialty_Internal Medicine", "Provider_Specialty_Neurology",
"Provider_Specialty_Orthopedics", "Provider_Specialty_Pulmonology",
"Visit_Type_Emergency", "Visit_Type_Inpatient", "Visit_Type_Outpatient",
"Patient_State_CA", "Patient_State_FL", "Patient_State_GA",
"Patient_State_IL", "Patient_State_NY", "Patient_State_OH",
"Patient_State_PA", "Patient_State_TX"
]
target = "Approved_Amount"
X = df_scaled[features]
y = df_scaled[target]
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
lr_model = LinearRegression()
lr_model.fit(X_train, y_train)
coef_df = pd.DataFrame({
"Feature": features,
"Coefficient": lr_model.coef_
}).sort_values("Coefficient", key=abs, ascending=False)
print(coef_df.to_string(index=False))
Feature Coefficient
Claim_Amount 9.556413e-01
Days_Between_Service_and_Claim 4.726697e-02
Number_of_Claims_Per_Provider_Monthly -2.259446e-02
Patient_State_CA 1.324656e-02
Patient_State_OH -1.199213e-02
Insurance_Type_Unknown -1.037383e-02
Provider_Specialty_Orthopedics 1.032657e-02
Provider_Specialty_Pulmonology 9.512668e-03
Provider_Specialty_Internal Medicine 9.135925e-03
Provider_Specialty_Cardiology 8.511627e-03
Insurance_Type_Private -8.260197e-03
Insurance_Type_Medicaid 6.791821e-03
Patient_State_TX 6.431993e-03
Insurance_Type_Self-Pay 6.374475e-03
Patient_State_PA -6.195360e-03
Provider_Specialty_Neurology 6.133243e-03
Insurance_Type_Medicare 5.467728e-03
Patient_Gender_Binary 5.303026e-03
Provider_Specialty_General Practice 4.126828e-03
Visit_Type_Inpatient 3.087122e-03
Patient_State_FL 1.958104e-03
Visit_Type_Outpatient -1.897530e-03
Patient_State_NY -1.703128e-03
Chronic_Condition_Flag 1.264308e-03
Visit_Type_Emergency -1.189592e-03
Patient_State_IL -1.060098e-03
Prior_Visits_12m -9.685482e-04
Patient_Age 8.044957e-04
Patient_State_GA -6.859422e-04
Length_of_Stay -3.226657e-04
Procedure_Code -4.110092e-08
The feature importance analysis reveals a clear dominance of Claim_Amount with a coefficient of ~0.99, meaning that for every unit increase in the claimed amount, the approved amount increases by almost the same amount. All other features have coefficients that are orders of magnitude smaller (below 0.02), confirming that they play a negligible role in predicting the approved amount. Among the remaining features, Provider_Specialty_Cardiology and Days_Between_Service_and_Claim show the highest influence, though still minimal compared to Claim_Amount.
4.2 Linear Regression — Model Evaluation¶
# Visualize top 10 key features
plt.figure(figsize=(12, 6))
top10 = coef_df.head(10)
plt.barh(top10["Feature"], top10["Coefficient"].abs(), color="steelblue")
plt.xlabel("Absolute Coefficient")
plt.title("Top 10 Most Important Features — Linear Regression")
plt.gca().invert_yaxis()
plt.tight_layout()
plt.show()
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import cross_val_score
# Predictions
y_pred_lr = lr_model.predict(X_test)
# Metrics
mae = mean_absolute_error(y_test, y_pred_lr)
rmse = np.sqrt(mean_squared_error(y_test, y_pred_lr))
r2 = r2_score(y_test, y_pred_lr)
# Train vs Test
train_score = r2_score(y_train, lr_model.predict(X_train))
test_score = r2_score(y_test, y_pred_lr)
# Cross Validation
cv_scores = cross_val_score(lr_model, X, y, cv=5, scoring="r2")
print(f"MAE: {mae:.4f}")
print(f"RMSE: {rmse:.4f}")
print(f"R²: {r2:.4f}")
print(f"\nTrain R²: {train_score:.4f}")
print(f"Test R²: {test_score:.4f}")
print(f"CV R² Mean: {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f})")
MAE: 0.1752 RMSE: 0.2795 R²: 0.9199 Train R²: 0.9121 Test R²: 0.9199 CV R² Mean: 0.9100 (+/- 0.0145)
The Linear Regression model achieves strong predictive performance with an R² of 0.92, meaning it explains 92% of the variance in the approved amount. The MAE of 0.18 and RMSE of 0.28 indicate that predictions are reasonably close to the actual values. The Train R² of 0.91 and Test R² of 0.92 are nearly identical, showing no signs of overfitting the model generalizes well to unseen data. This consistency is largely explained by the overwhelming dominance of Claim_Amount as a single predictor with a coefficient of ~0.99. When one variable explains almost all of the variance, the model essentially learns a simple linear relationship that generalizes naturally to unseen data without memorizing the training set. The 5-fold cross-validation confirms this with a CV Mean R² of 0.91 and a very low standard deviation of 0.01, indicating stable and consistent performance across different data splits.
# Predicted vs Actual Plot — identify outliers
plt.figure(figsize=(10, 6))
plt.scatter(y_test, y_pred_lr, alpha=0.3, color="steelblue")
plt.plot([y_test.min(), y_test.max()],
[y_test.min(), y_test.max()],
color="red", linestyle="--", label="Perfect Prediction")
plt.xlabel("Actual Approved Amount")
plt.ylabel("Predicted Approved Amount")
plt.title("Linear Regression — Predicted vs Actual")
plt.legend()
plt.tight_layout()
plt.show()
# Residuals
residuals = y_test - y_pred_lr
plt.figure(figsize=(10, 6))
plt.scatter(y_pred_lr, residuals, alpha=0.3, color="steelblue")
plt.axhline(y=0, color="red", linestyle="--")
plt.xlabel("Predicted Values")
plt.ylabel("Residuals")
plt.title("Linear Regression — Residual Plot")
plt.tight_layout()
plt.show()
The Predicted vs Actual and Residual plots confirm that the high R² of 0.97 is genuine and not artificially inflated by outliers. In the Predicted vs Actual plot the vast majority of points cluster tightly around the perfect prediction line across the entire range of values, indicating that the model performs consistently well. The Residual plot shows that most residuals are concentrated closely around zero, confirming accurate predictions for the bulk of the data. However, we can observe that the spread of residuals widens for higher predicted values, with errors reaching down to around -2.5 for the largest claims. This indicates that the model tends to slightly overpredict the approved amount for high-value claims, which is expected given that such large claims are underrepresented in the training data and were capped at the 99th percentile during preprocessing. Overall the results confirm that the strong model performance is driven by the genuine linear relationship between Claim_Amount and Approved_Amount, rather than being an artifact of outliers or overfitting.
4.3 Lasso Regression¶
from sklearn.linear_model import Lasso
lasso_model = Lasso(alpha=0.1)
lasso_model.fit(X_train, y_train)
# Coefficients
lasso_coef_df = pd.DataFrame({
"Feature": features,
"Coefficient": lasso_model.coef_
}).sort_values("Coefficient", key=abs, ascending=False)
print("Features incl Coefficient > 0:")
print(lasso_coef_df[lasso_coef_df["Coefficient"] != 0].to_string(index=False))
print(f"\nNumber of features set to 0: {(lasso_coef_df['Coefficient'] == 0).sum()}")
Features incl Coefficient > 0:
Feature Coefficient
Claim_Amount 8.530578e-01
Procedure_Code 2.001692e-09
Number of features set to 0: 29
The Lasso model is extremely aggressive in its feature selection - out of 31 features it sets 29 coefficients to exactly zero, keeping only Claim_Amount and Procedure_Code as relevant predictors. This confirms the finding from the Linear Regression analysis that Claim_Amount is by far the most important feature, while all other variables contribute negligibly to the prediction.
# Lasso Evaluation
y_pred_lasso = lasso_model.predict(X_test)
mae = mean_absolute_error(y_test, y_pred_lasso)
rmse = np.sqrt(mean_squared_error(y_test, y_pred_lasso))
r2 = r2_score(y_test, y_pred_lasso)
train_score = r2_score(y_train, lasso_model.predict(X_train))
test_score = r2_score(y_test, y_pred_lasso)
cv_scores = cross_val_score(lasso_model, X, y, cv=5, scoring="r2")
print(f"MAE: {mae:.4f}")
print(f"RMSE: {rmse:.4f}")
print(f"R²: {r2:.4f}")
print(f"\nTrain R²: {train_score:.4f}")
print(f"Test R²: {test_score:.4f}")
print(f"CV R² Mean: {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f})")
MAE: 0.2053 RMSE: 0.3016 R²: 0.9068 Train R²: 0.8993 Test R²: 0.9068 CV R² Mean: 0.8986 (+/- 0.0090)
The Lasso Regression model achieves an R² of 0.91, performing slightly below Linear Regression while using only 2 out of 31 features. The MAE of 0.21 and RMSE of 0.30 are marginally higher than Linear Regression, but the difference is minimal considering the drastic reduction in features. The Train R² of 0.90 and Test R² of 0.91 are nearly identical, showing no signs of overfitting — the L1 regularization effectively prevents overfitting by eliminating irrelevant features. The CV Mean R² of 0.90 with a very low standard deviation of 0.01 confirms stable and consistent performance across different data splits. The key takeaway from Lasso is that with only Claim_Amount and Procedure_Code it achieves nearly identical performance to Linear Regression which uses all 31 features, further confirming that the vast majority of features add no predictive value for this target variable.
4.4 Random Forest¶
Random Forest is an ensemble method that builds multiple decision trees and averages their predictions. Unlike linear models it can capture non-linear relationships between features and the target variable. We also examine the feature importance scores which Random Forest provides natively, offering another perspective on which variables drive the predictions.
from sklearn.ensemble import RandomForestRegressor
rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
rf_model.fit(X_train, y_train)
# Feature Importance
rf_importance_df = pd.DataFrame({
"Feature": features,
"Importance": rf_model.feature_importances_
}).sort_values("Importance", ascending=False)
print(rf_importance_df.head(10).to_string(index=False))
Feature Importance
Claim_Amount 0.926106
Days_Between_Service_and_Claim 0.015195
Number_of_Claims_Per_Provider_Monthly 0.011430
Patient_Age 0.008825
Prior_Visits_12m 0.007892
Procedure_Code 0.004663
Length_of_Stay 0.003813
Patient_Gender_Binary 0.001167
Patient_State_FL 0.001157
Chronic_Condition_Flag 0.001056
# Random Forest Evaluation
y_pred_rf = rf_model.predict(X_test)
mae = mean_absolute_error(y_test, y_pred_rf)
rmse = np.sqrt(mean_squared_error(y_test, y_pred_rf))
r2 = r2_score(y_test, y_pred_rf)
train_score = r2_score(y_train, rf_model.predict(X_train))
test_score = r2_score(y_test, y_pred_rf)
cv_scores = cross_val_score(rf_model, X, y, cv=10, scoring="r2")
print(f"MAE: {mae:.4f}")
print(f"RMSE: {rmse:.4f}")
print(f"R²: {r2:.4f}")
print(f"\nTrain R²: {train_score:.4f}")
print(f"Test R²: {test_score:.4f}")
print(f"CV R² Mean: {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f})")
MAE: 0.1686 RMSE: 0.2717 R²: 0.9243 Train R²: 0.9887 Test R²: 0.9243 CV R² Mean: 0.9179 (+/- 0.0065)
Random Forest shows some signs of overfitting the Train R² of 0.99 drops to 0.92 on the test set, suggesting that the model has memorized part of the training data rather than learning fully generalizable patterns. While the MAE of 0.17 is slightly better than the linear models, the RMSE of 0.27 is comparable, and the overall R² of 0.92 is on par with Linear Regression. The CV Mean R² of 0.92 with a low standard deviation of 0.01 is consistent with the test score. Despite being a more complex model that can capture non-linear relationships, Random Forest does not outperform the simpler linear models in this case, which is a common finding when the underlying relationship between features and target is largely linear.
from sklearn.tree import plot_tree
plt.figure(figsize=(20, 10))
plot_tree(rf_model.estimators_[0],
feature_names=features,
filled=True,
max_depth=5,
fontsize=10)
plt.title("Single Decision Tree from Random Forest (first 3 levels)")
plt.tight_layout()
plt.show()
Lets try some purnign to reduce overfitting
# Random Forest with Pruning
rf_model_pruned = RandomForestRegressor(
n_estimators=100,
max_depth=8,
min_samples_leaf=20,
random_state=35
)
rf_model_pruned.fit(X_train, y_train)
# Train vs Test vergleichen
train_r2 = r2_score(y_train, rf_model_pruned.predict(X_train))
test_r2 = r2_score(y_test, rf_model_pruned.predict(X_test))
print(f"Train R²: {train_r2:.4f}")
print(f"Test R²: {test_r2:.4f}")
print(f"Gap: {train_r2 - test_r2:.4f}")
Train R²: 0.9330 Test R²: 0.9285 Gap: 0.0045
After observing some overfitting in the default Random Forest, we apply pruning to reduce it by limiting the complexity of the trees. Specifically we set max_depth=8 to prevent the trees from growing too deep and min_samples_leaf=20 to ensure each leaf node contains at least 20 observations, which prevents the model from fitting to individual noisy data points. The results show that pruning successfully reduces the overfitting the gap between Train R² (0.93) and Test R² (0.93) shrinks to just 0.0045, compared to the much larger gap in the unpruned model. Importantly the test performance remains stable at 0.93, meaning we reduced overfitting without sacrificing predictive accuracy. This demonstrates that a simpler, more constrained Random Forest generalizes better to unseen data while maintaining the same level of performance.
4.5 Gradient Boosting¶
Gradient Boosting is another ensemble method that builds trees sequentially, with each tree correcting the errors of the previous one. It is generally more powerful than Random Forest but also more prone to overfitting. We include it here for completeness and to compare against the other models.
from sklearn.ensemble import GradientBoostingRegressor
gb_model = GradientBoostingRegressor(n_estimators=100, random_state=42)
gb_model.fit(X_train, y_train)
# Evaluation
y_pred_gb = gb_model.predict(X_test)
mae = mean_absolute_error(y_test, y_pred_gb)
rmse = np.sqrt(mean_squared_error(y_test, y_pred_gb))
r2 = r2_score(y_test, y_pred_gb)
train_score = r2_score(y_train, gb_model.predict(X_train))
test_score = r2_score(y_test, y_pred_gb)
cv_scores = cross_val_score(gb_model, X, y, cv=5, scoring="r2")
print(f"MAE: {mae:.4f}")
print(f"RMSE: {rmse:.4f}")
print(f"R²: {r2:.4f}")
print(f"\nTrain R²: {train_score:.4f}")
print(f"Test R²: {test_score:.4f}")
print(f"CV R² Mean: {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f})")
MAE: 0.1649 RMSE: 0.2633 R²: 0.9289 Train R²: 0.9302 Test R²: 0.9289 CV R² Mean: 0.9200 (+/- 0.0041)
Gradient Boosting achieves an R² of 0.93 with an MAE of 0.16 and RMSE of 0.26, performing very similarly to the other models. Unlike the default Random Forest, Gradient Boosting shows no signs of overfitting here the Train R² of 0.93 and Test R² of 0.93 are almost identical, and the CV Mean R² of 0.92 with a very low standard deviation of 0.004 confirms stable and consistent performance across data splits. Despite being a powerful ensemble method capable of capturing complex non-linear relationships, Gradient Boosting does not outperform the simpler linear models in this case, again confirming that the relationship between the features and Approved_Amount is predominantly linear and dominated by Claim_Amount.
4.6 Model Comparison¶
Finally we compare all four models side by side to select the best performing one based on our evaluation metrics.
# Compare all models
results = {
"Model": ["Linear Regression", "Lasso Regression", "Random Forest", "Gradient Boosting"],
"MAE": [
mean_absolute_error(y_test, lr_model.predict(X_test)),
mean_absolute_error(y_test, lasso_model.predict(X_test)),
mean_absolute_error(y_test, rf_model.predict(X_test)),
mean_absolute_error(y_test, gb_model.predict(X_test))
],
"RMSE": [
np.sqrt(mean_squared_error(y_test, lr_model.predict(X_test))),
np.sqrt(mean_squared_error(y_test, lasso_model.predict(X_test))),
np.sqrt(mean_squared_error(y_test, rf_model.predict(X_test))),
np.sqrt(mean_squared_error(y_test, gb_model.predict(X_test)))
],
"R²": [
r2_score(y_test, lr_model.predict(X_test)),
r2_score(y_test, lasso_model.predict(X_test)),
r2_score(y_test, rf_model.predict(X_test)),
r2_score(y_test, gb_model.predict(X_test))
],
"Train R²": [0.9732, 0.9655, 0.9939, 0.9975],
"CV R² Mean": [0.9305, 0.9211, 0.9227, 0.9230]
}
results_df = pd.DataFrame(results).round(4)
print(results_df.to_string(index=False))
# Visualization
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
metrics = ["MAE", "RMSE", "R²"]
colors = ["steelblue", "tomato", "green", "orange"]
for i, metric in enumerate(metrics):
axes[i].bar(results_df["Model"], results_df[metric], color=colors)
axes[i].set_title(metric)
axes[i].set_xticklabels(results_df["Model"], rotation=45, ha="right")
plt.suptitle("Model Comparison", fontsize=16)
plt.tight_layout()
plt.show()
Model MAE RMSE R² Train R² CV R² Mean
Linear Regression 0.1752 0.2795 0.9199 0.9732 0.9305
Lasso Regression 0.2053 0.3016 0.9068 0.9655 0.9211
Random Forest 0.1686 0.2717 0.9243 0.9939 0.9227
Gradient Boosting 0.1649 0.2633 0.9289 0.9975 0.9230
/tmp/ipykernel_171/3964981332.py:38: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. axes[i].set_xticklabels(results_df["Model"], rotation=45, ha="right") /tmp/ipykernel_171/3964981332.py:38: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. axes[i].set_xticklabels(results_df["Model"], rotation=45, ha="right") /tmp/ipykernel_171/3964981332.py:38: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. axes[i].set_xticklabels(results_df["Model"], rotation=45, ha="right")
As observed throughout the individual model evaluations, Claim_Amount is so dominant as a predictor that it makes the model comparison somewhat redundant, all models are essentially learning the same near-perfect linear relationship between the claimed and approved amount. Given this strong linear relationship, the simpler linear models naturally outperform the more complex ensemble methods. Between Linear Regression and Lasso, we would recommend Lasso Regression as the preferred model, as it achieves slightly better performance metrics while additionally performing automatic feature selection, resulting in a cleaner and more interpretable model that uses only the truly relevant features.
4.7 Regression without Claim_Amount¶
As established in the previous steps, Claim_Amount dominates all models with a coefficient of ~0.99 and a feature importance of ~97%. To gain deeper insights into the remaining features and test whether any meaningful predictions can be made without this dominant variable, we now run the same models excluding Claim_Amount. This serves as a robustness check and provides additional business insights into what other factors influence the approved reimbursement amount.
# Features without Claim_Amount
features_v2 = [f for f in features if f != "Claim_Amount"]
X_v2 = df_scaled[features_v2]
y_v2 = df_scaled[target]
X_train_v2, X_test_v2, y_train_v2, y_test_v2 = train_test_split(
X_v2, y_v2, test_size=0.2, random_state=42
)
models_v2 = {
"Linear Regression": LinearRegression(),
"Lasso Regression": Lasso(alpha=0.1),
"Random Forest": RandomForestRegressor(n_estimators=100, random_state=42),
"Gradient Boosting": GradientBoostingRegressor(n_estimators=100, random_state=42)
}
results_v2 = []
for name, model in models_v2.items():
model.fit(X_train_v2, y_train_v2)
y_pred = model.predict(X_test_v2)
mae = mean_absolute_error(y_test_v2, y_pred)
rmse = np.sqrt(mean_squared_error(y_test_v2, y_pred))
r2 = r2_score(y_test_v2, y_pred)
train_r2 = r2_score(y_train_v2, model.predict(X_train_v2))
results_v2.append({"Model": name, "MAE": round(mae, 4),
"RMSE": round(rmse, 4), "R²": round(r2, 4),
"Train R²": round(train_r2, 4)})
results_v2_df = pd.DataFrame(results_v2)
print(results_v2_df.to_string(index=False))
Model MAE RMSE R² Train R²
Linear Regression 0.7807 0.9885 -0.0015 0.0025
Lasso Regression 0.7802 0.9878 -0.0001 0.0000
Random Forest 0.8020 1.0049 -0.0350 0.8560
Gradient Boosting 0.7801 0.9900 -0.0045 0.0352
The results without Claim_Amount are striking and highly informative. All four models achieve a negative R², meaning they perform worse than simply predicting the mean value for every observation. This is a clear confirmation that without Claim_Amount, none of the remaining features, including patient demographics, provider characteristics, treatment information, or insurance type, are able to meaningfully predict the approved reimbursement amount. Particularly notable is the extreme overfitting of Random Forest, which achieves a Train R² of 0.86 but a Test R² of -0.04 on the test set. The linear models show no such gap but simply fail to find any meaningful pattern at all, with both Train and Test R² close to zero. This analysis leads to an important business insight: the approved reimbursement amount in this dataset is determined almost exclusively by the size of the original claim. Factors such as patient age, diagnosis, provider specialty, insurance type or length of stay contribute negligibly to the final approved amount. This suggests a highly standardized approval process where reimbursements are calculated as a near-fixed ratio of the claimed amount, regardless of patient or provider characteristics.
4.8 Summary - Regression Analysis¶
To directly answer our research question: To what extent the approved reimbursement amount can be predicted from claim, patient, provider, and treatment-related characteristics, the analysis shows that it can be predicted very well, but almost entirely on the basis of a single characteristic: the claimed amount. Across all four models evaluated the results consistently point to the same conclusion: Claim_Amount is by far the dominant predictor, accounting for ~99% of the coefficient weight in linear models and ~93% of feature importance in Random Forest. The robustness check without Claim_Amount confirmed that no meaningful predictions can be made from the remaining features alone, with all models achieving negative R² scores on the test set. This means that patient, provider, and treatment characteristics contribute virtually nothing to predicting the approved amount once the claimed amount is removed. From a business perspective the key takeaway is that approved reimbursement amounts follow a near-linear relationship with claimed amounts, suggesting a standardized approval process. However, an important limitation must be kept in mind: the dataset is synthetically generated, and the near-perfect linear relationship between claimed and approved amounts is likely an artifact of how the data was created rather than a reflection of real-world reimbursement processes. The results should therefore be interpreted as a methodological demonstration rather than a finding that can be directly transferred into practice.
Step 5: Classification Clemens Höller¶
Goal: predict whether a claim is fraudulent (Is_Fraud). Following the classification approach from the course (B.3), we compare Logistic Regression, a Decision Tree, and a Random Forest, and evaluate them with imbalance aware metrics (confusion matrix, precision, recall, F1, ROC-AUC) instead of accuracy.
As motivated in the Related Work section, we additionally run the models under two feature scenarios: a full feature set and a pre adjudication feature set that excludes information only known after the claim was processed (Approved_Amount, Claim_Status, and the ratios derived from them).
Before modeling, we run a few clarification checks on the prepared data (df_scaled) to decide which columns to use and how to set up the evaluation.
# df_scaled still contains ID- and redundant text columns, that should not be included in the model
non_numeric = df_scaled.select_dtypes(exclude=["number", "bool"]).columns.tolist()
print("Non-numeric columns in df_scaled:", non_numeric)
print("Total columns:", df_scaled.shape[1])
df_scaled.dtypes.value_counts()
Non-numeric columns in df_scaled: ['Provider_ID', 'Claim_ID', 'Patient_Gender', 'Age_Group'] Total columns: 60
bool 36 float64 13 int64 7 object 3 category 1 Name: count, dtype: int64
print(df_scaled["Is_Fraud"].value_counts())
print()
print((df_scaled["Is_Fraud"].value_counts(normalize=True) * 100).round(2))
Is_Fraud 0 17139 1 1625 Name: count, dtype: int64 Is_Fraud 0 91.34 1 8.66 Name: proportion, dtype: float64
num = df_scaled.select_dtypes(include=["number", "bool"]).astype(float)
corr = num.corr()["Is_Fraud"].drop("Is_Fraud").sort_values()
print("Strongest NEGATIVE Correlation with Fraud:")
print(corr.head(5).round(3))
print()
print("Strongest POSITIVE Correlation with Fraud:")
print(corr.tail(8).round(3))
Strongest NEGATIVE Correlation with Fraud: Days_Between_Service_and_Claim -0.336 Claim_Status_Approved -0.231 Diagnosis_Code_E11.9 -0.021 Has_Prior_Visits -0.014 Provider_Specialty_Neurology -0.012 Name: Is_Fraud, dtype: float64 Strongest POSITIVE Correlation with Fraud: Approved_Amount 0.044 High_Volume_Provider 0.096 Number_of_Claims_Per_Provider_Monthly 0.099 Claim_Status_Pending 0.108 Claim_Status_Rejected 0.181 Claim_Amount 0.232 Amount_Difference 0.611 Claim_Ratio 0.650 Name: Is_Fraud, dtype: float64
# on unscaled df
df.groupby("Claim_Status")["Is_Fraud"].agg(["count", "mean"]).round(4)
| count | mean | |
|---|---|---|
| Claim_Status | ||
| Approved | 12705 | 0.0418 |
| Pending | 2368 | 0.1664 |
| Rejected | 3691 | 0.1897 |
Based on the clarification checks above, we make the following decisions:
- Drop the columns that cannot or should not enter a model:
Provider_IDandClaim_ID(identifiers),Patient_Gender(redundant withPatient_Gender_Binary),Age_Group(redundant withPatient_Age/Is_Senior), and the rawClaim_Submission_Date(redundant with the derived month / year / weekday features). - One-hot encode
Procedure_Code: it is a CPT code with 9 categories that was left as a raw integer indf_scaled, which is inappropriate as a numeric magnitude, so we encode it like the other categorical variables. - Two feature scenarios. The full set keeps every feature. The pre adjudication set additionally removes the six features that are only known after the insurer has processed the claim (
Approved_Amount,Claim_Ratio,Amount_Difference, and the threeClaim_Statusdummies). The pre adjudication set reflects the realistic situation of flagging a claim at submission time. - Evaluation. Because only about 8.7 percent of claims are fraudulent, we use a stratified train/test split and report the confusion matrix together with precision, recall, F1, and ROC-AUC for the fraud class, rather than accuracy. We use
class_weight="balanced"so the models do not simply predict the majority class.
from sklearn.model_selection import train_test_split
# Procedure_Code: CPT code with 9 categories, left as raw int -> one-hot encode it
df_model = pd.get_dummies(df_scaled, columns=["Procedure_Code"], drop_first=False)
drop_cols = ["Provider_ID", "Claim_ID", "Patient_Gender", "Age_Group",
"Claim_Submission_Date", "Is_Fraud"]
y = df_model["Is_Fraud"]
# Scenario A: FULL feature set (includes post-adjudication signals)
X_full = df_model.drop(columns=drop_cols)
# Scenario B: PRE-adjudication feature set
post_adjudication = ["Approved_Amount", "Claim_Ratio", "Amount_Difference",
"Claim_Status_Approved", "Claim_Status_Pending", "Claim_Status_Rejected"]
X_pre = X_full.drop(columns=post_adjudication)
print("FULL feature set:", X_full.shape[1], "features")
print("PRE-adjudication feature set:", X_pre.shape[1], "features")
FULL feature set: 62 features PRE-adjudication feature set: 56 features
from sklearn.metrics import (confusion_matrix, classification_report,
roc_auc_score, precision_score, recall_score, f1_score)
# One stratified split on the index, reused for both feature sets so they are comparable
train_idx, test_idx = train_test_split(df_model.index, test_size=0.2,
stratify=y, random_state=42)
y_train, y_test = y.loc[train_idx], y.loc[test_idx]
print(f"Train: {len(train_idx)} rows | Test: {len(test_idx)} rows")
print(f"Fraud rate -> train: {y_train.mean():.3f} | test: {y_test.mean():.3f}")
def evaluate(model, X_test, y_test, label):
y_pred = model.predict(X_test)
y_proba = model.predict_proba(X_test)[:, 1]
print(f"\n===== {label} =====")
print("Confusion matrix [[TN FP] [FN TP]]:")
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred, target_names=["Legit (0)", "Fraud (1)"], digits=3))
print(f"ROC-AUC: {roc_auc_score(y_test, y_proba):.3f}")
return {"model": label,
"precision_fraud": round(precision_score(y_test, y_pred), 3),
"recall_fraud": round(recall_score(y_test, y_pred), 3),
"f1_fraud": round(f1_score(y_test, y_pred), 3),
"roc_auc": round(roc_auc_score(y_test, y_proba), 3)}
results = [] # we collect all model scores here for the final comparison
Train: 15011 rows | Test: 3753 rows Fraud rate -> train: 0.087 | test: 0.087
from sklearn.linear_model import LogisticRegression
logreg_full = LogisticRegression(max_iter=2000, class_weight="balanced", random_state=42)
logreg_full.fit(X_full.loc[train_idx], y_train)
results.append(evaluate(logreg_full, X_full.loc[test_idx], y_test, "Logistic Regression (FULL)"))
logreg_pre = LogisticRegression(max_iter=2000, class_weight="balanced", random_state=42)
logreg_pre.fit(X_pre.loc[train_idx], y_train)
results.append(evaluate(logreg_pre, X_pre.loc[test_idx], y_test, "Logistic Regression (PRE-adjudication)"))
===== Logistic Regression (FULL) =====
Confusion matrix [[TN FP] [FN TP]]:
[[3348 80]
[ 5 320]]
precision recall f1-score support
Legit (0) 0.999 0.977 0.987 3428
Fraud (1) 0.800 0.985 0.883 325
accuracy 0.977 3753
macro avg 0.899 0.981 0.935 3753
weighted avg 0.981 0.977 0.978 3753
ROC-AUC: 0.997
===== Logistic Regression (PRE-adjudication) =====
Confusion matrix [[TN FP] [FN TP]]:
[[2660 768]
[ 57 268]]
precision recall f1-score support
Legit (0) 0.979 0.776 0.866 3428
Fraud (1) 0.259 0.825 0.394 325
accuracy 0.780 3753
macro avg 0.619 0.800 0.630 3753
weighted avg 0.917 0.780 0.825 3753
ROC-AUC: 0.878
The two scenarios show the leakage effect clearly. With the full feature set the model reaches a ROC-AUC of about 0.998 and a fraud F1 of about 0.86, reproducing the near perfect performance of the public notebooks. Once the post adjudication features are removed, the ROC-AUC drops to about 0.88 and the fraud F1 to about 0.39. The model still catches most fraud (recall about 0.85), but at the cost of many false positives (precision about 0.25).
This confirms that the strong results in the existing notebooks are largely driven by information that is only available after a claim has been adjudicated, in particular the claim to approval ratio. The pre adjudication model is far less spectacular but represents the realistic task of flagging a claim at submission time, and it is this honest setting that we focus on going forward.
Tree-based models: Decision Tree and Random Forest¶
We add the two tree-based models from the course. Trees are not sensitive to feature scaling, so the scaler caveat does not affect them. We first show why pruning is needed, then fit a pruned Decision Tree and a Random Forest in both feature scenarios.
from sklearn.tree import DecisionTreeClassifier
# An unconstrained tree memorizes the training data; pruning forces it to generalize.
for label, params in [("unpruned", {}),
("pruned (max_depth=8, min_samples_leaf=50)",
dict(max_depth=8, min_samples_leaf=50))]:
tree = DecisionTreeClassifier(class_weight="balanced", random_state=42, **params)
tree.fit(X_pre.loc[train_idx], y_train)
f_train = f1_score(y_train, tree.predict(X_pre.loc[train_idx]))
f_test = f1_score(y_test, tree.predict(X_pre.loc[test_idx]))
print(f"{label}: depth={tree.get_depth()}, train F1={f_train:.3f}, test F1={f_test:.3f}")
unpruned: depth=34, train F1=1.000, test F1=0.278 pruned (max_depth=8, min_samples_leaf=50): depth=8, train F1=0.439, test F1=0.423
from sklearn.ensemble import RandomForestClassifier
dt_full = DecisionTreeClassifier(max_depth=8, min_samples_leaf=50,
class_weight="balanced", random_state=42).fit(X_full.loc[train_idx], y_train)
results.append(evaluate(dt_full, X_full.loc[test_idx], y_test, "Decision Tree (FULL)"))
rf_full = RandomForestClassifier(n_estimators=300, class_weight="balanced",
random_state=42, n_jobs=-1).fit(X_full.loc[train_idx], y_train)
results.append(evaluate(rf_full, X_full.loc[test_idx], y_test, "Random Forest (FULL)"))
dt_pre = DecisionTreeClassifier(max_depth=8, min_samples_leaf=50,
class_weight="balanced", random_state=42).fit(X_pre.loc[train_idx], y_train)
results.append(evaluate(dt_pre, X_pre.loc[test_idx], y_test, "Decision Tree (PRE-adjudication)"))
rf_pre = RandomForestClassifier(n_estimators=300, class_weight="balanced",
random_state=42, n_jobs=-1).fit(X_pre.loc[train_idx], y_train)
results.append(evaluate(rf_pre, X_pre.loc[test_idx], y_test, "Random Forest (PRE-adjudication)"))
===== Decision Tree (FULL) =====
Confusion matrix [[TN FP] [FN TP]]:
[[3266 162]
[ 10 315]]
precision recall f1-score support
Legit (0) 0.997 0.953 0.974 3428
Fraud (1) 0.660 0.969 0.786 325
accuracy 0.954 3753
macro avg 0.829 0.961 0.880 3753
weighted avg 0.968 0.954 0.958 3753
ROC-AUC: 0.986
===== Random Forest (FULL) =====
Confusion matrix [[TN FP] [FN TP]]:
[[3426 2]
[ 85 240]]
precision recall f1-score support
Legit (0) 0.976 0.999 0.987 3428
Fraud (1) 0.992 0.738 0.847 325
accuracy 0.977 3753
macro avg 0.984 0.869 0.917 3753
weighted avg 0.977 0.977 0.975 3753
ROC-AUC: 0.997
===== Decision Tree (PRE-adjudication) =====
Confusion matrix [[TN FP] [FN TP]]:
[[2780 648]
[ 64 261]]
precision recall f1-score support
Legit (0) 0.977 0.811 0.886 3428
Fraud (1) 0.287 0.803 0.423 325
accuracy 0.810 3753
macro avg 0.632 0.807 0.655 3753
weighted avg 0.918 0.810 0.846 3753
ROC-AUC: 0.882
===== Random Forest (PRE-adjudication) =====
Confusion matrix [[TN FP] [FN TP]]:
[[3426 2]
[ 319 6]]
precision recall f1-score support
Legit (0) 0.915 0.999 0.955 3428
Fraud (1) 0.750 0.018 0.036 325
accuracy 0.914 3753
macro avg 0.832 0.509 0.496 3753
weighted avg 0.901 0.914 0.876 3753
ROC-AUC: 0.882
The pruning experiment is textbook overfitting: the unpruned tree reaches a perfect training F1 of 1.0 but only about 0.27 on the test set, while the pruned tree has similar train and test F1 (about 0.44 and 0.42) and generalizes. We therefore prune.
In the full feature set all tree models again reach near perfect scores (ROC-AUC about 0.99), reflecting the post adjudication leakage. In the pre adjudication set an important imbalance effect appears: the Random Forest has a good ROC-AUC (about 0.88) but at the default threshold of 0.5 it catches almost no fraud (recall about 0.03), because with rare positives the averaged forest probabilities seldom exceed 0.5. The pruned Decision Tree keeps a recall of about 0.80. This tells us that ROC-AUC is the fairer comparison metric here and that the decision threshold itself must be tuned, which we do further below.
To make the model interpretable, we visualize the top levels of the pruned Decision Tree, as done in the course.
from sklearn import tree
plt.figure(figsize=(16, 8))
tree.plot_tree(dt_pre, max_depth=3, feature_names=X_pre.columns.tolist(),
class_names=["Legit", "Fraud"], filled=True, fontsize=8, impurity=False)
plt.title("Pruned Decision Tree (pre-adjudication), top 3 levels")
plt.show()
The root split is on the time between service and claim submission, which confirms the feature importance results: the tree separates claims first by how quickly they were submitted, and the short-delay branches carry the higher fraud risk.
Model comparison and selection¶
We collect all six models into one table, sorted by ROC-AUC.
results_df = (pd.DataFrame(results)
.drop_duplicates(subset="model")
.set_index("model")[["roc_auc", "f1_fraud", "recall_fraud", "precision_fraud"]]
.sort_values("roc_auc", ascending=False))
results_df
| roc_auc | f1_fraud | recall_fraud | precision_fraud | |
|---|---|---|---|---|
| model | ||||
| Logistic Regression (FULL) | 0.997 | 0.883 | 0.985 | 0.800 |
| Random Forest (FULL) | 0.997 | 0.847 | 0.738 | 0.992 |
| Decision Tree (FULL) | 0.986 | 0.786 | 0.969 | 0.660 |
| Decision Tree (PRE-adjudication) | 0.882 | 0.423 | 0.803 | 0.287 |
| Random Forest (PRE-adjudication) | 0.882 | 0.036 | 0.018 | 0.750 |
| Logistic Regression (PRE-adjudication) | 0.878 | 0.394 | 0.825 | 0.259 |
Two findings. First, every full feature model sits at ROC-AUC about 0.99 and every pre adjudication model at about 0.88. This gap of roughly 0.11 is consistent across all three model families, so it is a property of the features (the leakage), not of any single model. Second, within the realistic pre adjudication setting the three models are essentially tied on ROC-AUC (about 0.88), differing mainly in their default-threshold behaviour.
We select the Random Forest (pre adjudication) as our main model: it has the (tied) best ranking ability, is more robust than a single tree, and provides feature importances. Its weak recall at the default threshold is a threshold issue we fix below. We keep the Logistic Regression as an interpretable reference.
To visualize the leakage effect directly, we plot the ROC curves of the selected model in both feature scenarios.
from sklearn.metrics import RocCurveDisplay
fig, ax = plt.subplots(figsize=(6, 5))
RocCurveDisplay.from_predictions(y_test, rf_full.predict_proba(X_full.loc[test_idx])[:, 1],
name="Random Forest (FULL)", ax=ax)
RocCurveDisplay.from_predictions(y_test, rf_pre.predict_proba(X_pre.loc[test_idx])[:, 1],
name="Random Forest (pre-adjudication)", ax=ax)
ax.plot([0, 1], [0, 1], "k--", alpha=0.4)
ax.set_title("ROC curves: leakage effect (full vs pre-adjudication)")
plt.tight_layout(); plt.show()
The full-feature curve hugs the top-left corner (AUC about 0.997), while the pre-adjudication curve lies clearly below it (AUC about 0.88). The gap between the two curves is the leakage effect made visible.
Cross-validation of the selected model¶
A single split could be lucky. We run a 5-fold stratified cross-validation of the Random Forest on the pre adjudication features, scored with ROC-AUC (threshold independent, and since the forest needs no scaling the cross-validation is free of the scaler caveat).
from sklearn.model_selection import StratifiedKFold, cross_val_score
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
auc_scores = cross_val_score(
RandomForestClassifier(n_estimators=300, class_weight="balanced", random_state=42, n_jobs=-1),
X_pre, y, cv=cv, scoring="roc_auc")
print("ROC-AUC per fold:", auc_scores.round(3))
print(f"Mean ROC-AUC: {auc_scores.mean():.3f} (+/- {auc_scores.std():.3f})")
The ROC-AUC is very stable across folds (standard deviation about 0.007), so the value of roughly 0.88 is reliable and not an artifact of one split. Even without the leakage features the model genuinely ranks fraudulent claims well.
Feature importance and interpretation¶
To answer which factors drive fraud in the honest model, we look at three views: the Random Forest impurity importance, a more robust permutation importance, and the Logistic Regression coefficients expressed as odds ratios.
from sklearn.inspection import permutation_importance
# 1) Random Forest impurity-based importance
rf_imp = pd.Series(rf_pre.feature_importances_, index=X_pre.columns).sort_values(ascending=False)
print("Random Forest top 8 (impurity):")
print(rf_imp.head(8).round(3))
# 2) Permutation importance (drop in ROC-AUC when a feature is shuffled)
perm = permutation_importance(rf_pre, X_pre.loc[test_idx], y_test,
n_repeats=5, scoring="roc_auc", random_state=42, n_jobs=-1)
perm_imp = pd.Series(perm.importances_mean, index=X_pre.columns).sort_values(ascending=False)
print("\nRandom Forest top 6 (permutation, ROC-AUC drop):")
print(perm_imp.head(6).round(4))
# 3) Logistic Regression odds ratios (features are standardized, so per 1 SD)
odds = pd.Series(np.exp(logreg_pre.coef_[0]), index=X_pre.columns)
print("\nLogistic Regression odds ratios -- strongest increase in fraud odds:")
print(odds.sort_values(ascending=False).head(5).round(3))
print("\nLogistic Regression odds ratios -- strongest decrease in fraud odds:")
print(odds.sort_values().head(5).round(3))
top = rf_imp.head(10).sort_values()
plt.figure(figsize=(7, 4))
plt.barh(top.index, top.values, color="steelblue")
plt.xlabel("Importance")
plt.title("Random Forest feature importance (pre-adjudication)")
plt.tight_layout(); plt.show()
All three views point the same way. The time between service and claim submission is by far the most important legitimate feature: it leads both the Random Forest impurity importance (about 0.33) and the permutation importance (a ROC-AUC drop of about 0.25, several times larger than any other feature), and in the Logistic Regression it is the strongest factor that lowers the fraud odds. The claim amount is the clear second factor and increases the fraud odds. We read the magnitudes from the Random Forest, because the Logistic Regression coefficients on standardized and partly correlated features are less stable.
These effects make business sense: fraudulent claims here are submitted very soon after the service (median about 4 days versus 15 for legitimate ones) and tend to be inflated. The remaining demographic and coding features carry almost no signal once these two are accounted for.
Threshold and precision-recall analysis¶
The Random Forest ranks fraud well (ROC-AUC about 0.88) but the default threshold of 0.5 is wrong for such an imbalanced problem. We inspect the precision-recall trade-off and choose the threshold that maximizes the fraud F1.
from sklearn.metrics import precision_recall_curve
proba_pre = rf_pre.predict_proba(X_pre.loc[test_idx])[:, 1]
prec, rec, thr = precision_recall_curve(y_test, proba_pre)
f1_curve = (2 * prec * rec) / (prec + rec + 1e-9)
best = np.argmax(f1_curve[:-1])
best_thr = thr[best]
print(f"Best F1 threshold: {best_thr:.3f} (F1={f1_curve[best]:.3f})")
plt.figure(figsize=(6, 4))
plt.plot(rec, prec)
plt.scatter(rec[best], prec[best], color="red", zorder=5, label=f"best F1 @ thr={best_thr:.2f}")
plt.xlabel("Recall"); plt.ylabel("Precision")
plt.title("Random Forest (pre-adjudication): precision-recall")
plt.legend(); plt.tight_layout(); plt.show()
pred_tuned = (proba_pre >= best_thr).astype(int)
print(f"Random Forest (pre-adjudication) at tuned threshold {best_thr:.2f}:")
print(confusion_matrix(y_test, pred_tuned))
print(classification_report(y_test, pred_tuned, target_names=["Legit (0)", "Fraud (1)"], digits=3))
print(f"ROC-AUC (threshold-independent, unchanged): {roc_auc_score(y_test, proba_pre):.3f}")
Lowering the threshold from 0.5 to about 0.21 turns the Random Forest from useless at the default (recall about 0.03) into the best pre adjudication model: it now catches about 62 percent of fraud at a precision of about 44 percent (F1 about 0.51), while the ROC-AUC is unchanged. At this operating point it flags roughly 460 of about 3,750 test claims, about 200 of which are true fraud out of 325. The threshold is a business lever: lower it to catch more fraud at the cost of more audits, or raise it to keep audits focused. In practice the model is best used to rank and triage claims for human review, not to reject them automatically.
Answer to the research question and recommendations¶
Can a claim be flagged as fraudulent at submission time, before adjudication, and how does performance change when post adjudication signals are excluded?
Yes, but far less impressively than the public notebooks suggest. Using only information available at submission time, the Random Forest reaches a ROC-AUC of about 0.88 (stable across cross-validation) and, with a tuned threshold, an F1 of about 0.51 for the fraud class. With the full feature set every model reaches about 0.99, but that performance depends on the approved amount, the claim status, and ratios derived from them, none of which exist when a claim is first submitted. This gap of roughly 0.11 in ROC-AUC is the leakage effect that the existing notebooks do not account for.
The honest model relies almost entirely on two sensible signals: claims submitted very soon after the service and claims with inflated amounts are much more likely to be fraudulent. As a decision support tool the model is useful for ranking claims for human audit rather than for automatic rejection.