# -*- coding: utf-8 -*-
"""VisaApplication_vFinal.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1Uc8ptXpI8EA_wO9jokKJIKB-BbXAMfzW

<a href="https://colab.research.google.com/github/EvagAIML/Case-Study-V3/blob/main/FEasyVisa_Full_Code_vFinal_ErikVagshenian.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

### Executive Summary

In order to improve the fairness and efficiency of the U.S. visa certification process, this analysis uses data from 25,480 historical visa applications to predict whether a visa will be certified or denied. Each record includes both employee and employer details, such as education level, job experience, prevailing wage, company size, and region of employment. The target variable, case_status, represents visa outcomes — Certified (1) or Denied (0) — making this a supervised binary classification problem.

The objective is to support the Office of Foreign Labor Certification (OFLC) by identifying which applications are most likely to be approved or denied, enabling faster and more consistent decision-making.
The final Gradient Boosting model, optimized through class weighting, resampling, and threshold tuning, achieved 70.3% accuracy for certified cases and 69.7% accuracy for denied cases, with an overall accuracy of 70.1%, F1 score of 0.76, and ROC-AUC of 0.77.

The F1 score of 0.76 is acceptable because it reflects a balanced tradeoff between precision (avoiding false approvals) and recall (catching true certifications), which is critical in this context where both incorrect approvals and incorrect denials carry significant operational and compliance risks.

Outcome: A balanced machine learning model was developed that correctly identifies about 7 in 10 visa certifications and denials, offering a fair and data-driven decision-support system that improves processing efficiency while maintaining compliance integrity.
"""



"""### Problem Statement

**Business Context**

Business communities in the United States are facing high demand for human resources, but one of the constant challenges is identifying and attracting the right talent, which is perhaps the most important element in remaining competitive. Companies in the United States look for hard-working, talented, and qualified individuals both locally as well as abroad.

The Immigration and Nationality Act (INA) of the US permits foreign workers to come to the United States to work on either a temporary or permanent basis. The act also protects US workers against adverse impacts on their wages or working conditions by ensuring US employers’ compliance with statutory requirements when they hire foreign workers to fill workforce shortages. The immigration programs are administered by the Office of Foreign Labor Certification (OFLC).

OFLC processes job certification applications for employers seeking to bring foreign workers into the United States and grants certifications in those cases where employers can demonstrate that there are not sufficient US workers available to perform the work at wages that meet or exceed the wage paid for the occupation in the area of intended employment.

**Objective**

In FY 2016, the OFLC processed 775,979 employer applications for 1,699,957 positions for temporary and permanent labor certifications. This was a nine percent increase in the overall number of processed applications from the previous year. The process of reviewing every case is becoming a tedious task as the number of applicants is increasing every year.

The increasing number of applicants every year calls for a Machine Learning based solution that can help in shortlisting the candidates having higher chances of VISA approval. OFLC has hired the firm EasyVisa for data-driven solutions. You as a data scientist at EasyVisa have to analyze the data provided and, with the help of a classification model:

Facilitate the process of visa approvals.
Recommend a suitable profile for the applicants for whom the visa should be certified or denied based on the drivers that significantly influence the case status.

**Data Description**

The data contains the different attributes of employee and the employer. The detailed data dictionary is given below.

case_id: ID of each visa application
continent: Information of continent the employee
education_of_employee: Information of education of the employee
has_job_experience: Does the employee has any job experience? Y= Yes; N = No
requires_job_training: Does the employee require any job training? Y = Yes; N = No
no_of_employees: Number of employees in the employer’s company
yr_of_estab: Year in which the employer’s company was established
region_of_employment: Information of foreign worker’s intended region of employment in the US.
prevailing_wage: Average wage paid to similarly employed workers in a specific occupation in the area of intended employment. The purpose of the prevailing wage is to ensure that the foreign worker is not underpaid compared to other workers offering the same or similar service in the same area of employment.
unit_of_wage: Unit of prevailing wage. Values include Hourly, Weekly, Monthly, and Yearly.
full_time_position: Is the position of work full-time? Y = Full Time Position; N = Part Time Position
case_status: Flag indicating if the Visa was certified or denied
"""



"""### 0.	Library Installation

**Outcome**: All libraries were successfully installed, creating a reliable technical environment for end-to-end data analysis and model development.

In order to enable comprehensive analysis and modeling, all necessary Python libraries are installed, including those for data manipulation, visualization, and machine learning such as pandas, NumPy, scikit-learn, XGBoost, and imbalanced-learn. This ensures that every dependency required for numerical computation, data cleaning, and classification modeling is properly configured for execution.
"""

# ============================================================
# 0) Library Installation
# ============================================================
!pip -q install --user numpy pandas scikit-learn matplotlib seaborn xgboost imbalanced-learn



"""### 2.	Load Dataset

**Outcome**: The dataset is successfully loaded with 25,480 records and 12 attributes, establishing the foundation for classification modeling.

In order to begin exploration, the visa dataset is loaded, containing 25,480 observations across 12 variables representing employee and employer attributes. These include numeric fields such as company size, year of establishment, and prevailing wage, and categorical fields such as continent, education, and region of employment. The target variable, case_status, indicates whether each visa was Certified or Denied, defining this as a binary classification problem.
"""

# ============================================================
# 1) Imports and Configuration
# ============================================================
import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import train_test_split, StratifiedKFold, RandomizedSearchCV
from sklearn.metrics import (
    accuracy_score, precision_score, recall_score, f1_score, roc_auc_score,
    confusion_matrix, classification_report
)
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.utils.class_weight import compute_class_weight

from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import BaggingClassifier, RandomForestClassifier, AdaBoostClassifier, GradientBoostingClassifier

from imblearn.pipeline import Pipeline as ImbPipeline
from imblearn.over_sampling import SMOTE
from imblearn.under_sampling import RandomUnderSampler

try:
    from xgboost import XGBClassifier
    XGB_OK = True
except Exception:
    XGB_OK = False

RANDOM_STATE = 42
np.random.seed(RANDOM_STATE)
sns.set_theme(style="whitegrid", context="notebook")



"""### 2.	Load Dataset

**Outcome**: The dataset is successfully loaded with 25,480 records and 12 attributes, establishing the foundation for classification modeling.

In order to begin exploration, the visa dataset is loaded, containing 25,480 observations across 12 variables representing employee and employer attributes. These include numeric fields such as company size, year of establishment, and prevailing wage, and categorical fields such as continent, education, and region of employment. The target variable, case_status, indicates whether each visa was Certified or Denied, defining this as a binary classification problem.
"""

# ============================================================
# 2) Data Loading
# ============================================================
URL = "https://raw.githubusercontent.com/EvagAIML/Case-Study-V3/refs/heads/main/EasyVisa%20(13).csv"
df = pd.read_csv(URL)

print("Shape:", df.shape)
display(df.head())
print("\nFirst 5 rows:")
display(df.head())
print("\nLast 5 rows:")
display(df.tail())



"""### 3.	Data Understanding (Types, Missingness, Target Balance)

**Outcome**: The dataset is complete and balanced enough for modeling, with Certified ≈ 67% and Denied ≈ 33%.

In order to assess data integrity and understand class proportions, the dataset’s variable types, missing values, and target distribution are reviewed. No missing values are detected. The class balance shows that approximately 67% of visa applications are certified and 33% are denied, confirming moderate imbalance. Numeric fields display wide but valid ranges, indicating the need for normalization during preprocessing.
"""

# ============================================================
# 3) Data Overview
# ============================================================
print("Data Types:\n", df.dtypes)
print("\nMissing Values per Column:")
display(df.isna().sum().sort_values(ascending=False))

print("\nTarget Distribution:")
print(df["case_status"].value_counts())

print("\nNumeric Summary:")
display(df.select_dtypes(include="number").describe().T)

print("\nSample Categorical Values:")
for col in df.select_dtypes(include="object").columns:
    print(f"\n{col}:\n", df[col].value_counts().head(5))



"""### 4.	Primary Cleaning (Essential, Non-destructive)

**Outcome**: The dataset is cleaned and standardized, with consistent categorical and numeric features across all 25,480 records.

In order to ensure reliable modeling, categorical inconsistencies are standardized, invalid numeric values (such as negative employee counts) are corrected using medians, and implausible establishment years are removed. Prevailing wages are normalized to an annualized rate for comparability across pay units. The target variable case_status is encoded numerically (Denied = 0, Certified = 1).
"""

# ============================================================
# 4) Data Cleaning and Feature Preparation
# ============================================================
df_clean = df.copy()

def normalize_yn(x):
    s = str(x).strip().upper()
    if s in {"Y","YES","TRUE","T","1"}: return "YES"
    if s in {"N","NO","FALSE","F","0"}: return "NO"
    return s

for col in ["has_job_experience", "requires_job_training", "full_time_position"]:
    df_clean[col] = df_clean[col].map(normalize_yn)

for col in df_clean.select_dtypes(include="object").columns:
    df_clean[col] = df_clean[col].astype(str).str.strip().str.upper()

if "no_of_employees" in df_clean.columns:
    emp = pd.to_numeric(df_clean["no_of_employees"], errors="coerce")
    neg = emp < 0
    if neg.sum() > 0:
        df_clean.loc[neg, "no_of_employees"] = emp.median()
        print(f"Negative values corrected in 'no_of_employees': {neg.sum()}")

if {"prevailing_wage", "unit_of_wage"}.issubset(df_clean.columns):
    factors = {"HOUR":2080, "WEEK":52, "MONTH":12, "YEAR":1}
    df_clean["unit_of_wage"] = df_clean["unit_of_wage"].astype(str).str.upper().str.strip()
    df_clean["prevailing_wage_yearly"] = df_clean.apply(
        lambda r: pd.to_numeric(r["prevailing_wage"], errors="coerce") * factors.get(r["unit_of_wage"], 1), axis=1
    )
    print("Annualized prevailing wages created.")

if "yr_of_estab" in df_clean.columns:
    years = pd.to_numeric(df_clean["yr_of_estab"], errors="coerce")
    valid = years.between(1800, 2025)
    removed = (~valid).sum()
    df_clean = df_clean.loc[valid]
    print(f"Invalid establishment years removed: {removed}")

df_clean["case_status"] = df_clean["case_status"].map({"DENIED":0, "CERTIFIED":1})
df_clean.dropna(subset=["case_status"], inplace=True)
df_clean["case_status"] = df_clean["case_status"].astype(int)

if "case_id" in df_clean.columns:
    df_clean.drop(columns=["case_id"], inplace=True)

print("Final shape after cleaning:", df_clean.shape)
display(df_clean.head())



"""### 5.	Exploratory Data Analysis (Univariate)

**Outcome**: Certified cases form the majority, and key applicant traits include high education levels, experience, and full-time employment in major U.S. regions.

In order to understand the characteristics of each variable individually, distributions of numeric and categorical features are analyzed. Most applicants hold bachelor’s or master’s degrees, and the majority apply for full-time roles. Applications are concentrated in the Northeast, South, and West regions. Companies range widely in size and wage levels, with annual wages spanning from $2 to over $300,000. Most applicants have prior job experience, while few require training. Certified cases dominate, indicating potential class imbalance.
"""

# ============================================================
# 5) Exploratory Analysis (Univariate)
# ============================================================
num_cols = df_clean.select_dtypes(include="number").columns.tolist()
cat_cols = df_clean.select_dtypes(exclude="number").columns.tolist()

for c in [x for x in num_cols if x != "case_status"]:
    fig, axs = plt.subplots(2, 1, figsize=(10, 6), gridspec_kw={'height_ratios':[1,2]}, sharex=True)
    sns.boxplot(x=df_clean[c], ax=axs[0])
    sns.histplot(df_clean[c], kde=True, ax=axs[1])
    axs[0].set_title(f"Boxplot — {c}")
    axs[1].set_title(f"Histogram — {c}")
    plt.tight_layout(); plt.show()

for c in [x for x in cat_cols if x.upper() != "CASE_STATUS"]:
    plt.figure(figsize=(10,4))
    order = df_clean[c].value_counts().index
    sns.countplot(data=df_clean, x=c, order=order)
    plt.title(f"Counts — {c}")
    plt.xticks(rotation=45, ha="right")
    plt.tight_layout(); plt.show()



"""### 6.	Exploratory Data Analysis (Bivariate)

**Outcome**: Education, experience, region, and wage emerge as key determinants of visa certification outcomes.

In order to identify factors influencing approval outcomes, relationships between independent variables and case_status are examined. Applicants with master’s or doctorate degrees have notably higher certification rates than those with high school education. Job experience significantly increases the likelihood of certification, while training requirements decrease it. Regionally, the Northeast and West show higher certification rates, corresponding with higher prevailing wages.
"""

# ============================================================
# 6) Exploratory Analysis (Bivariate)
# ============================================================
target = "case_status"

for c in [x for x in num_cols if x != target]:
    fig, axs = plt.subplots(1, 2, figsize=(12, 4))
    sns.boxplot(data=df_clean, x=target, y=c, ax=axs[0])
    sns.kdeplot(data=df_clean, x=c, hue=target, common_norm=False, ax=axs[1])
    axs[0].set_title(f"{c} by Case Outcome")
    axs[1].set_title(f"Density of {c} by Case Outcome")
    plt.tight_layout(); plt.show()

def stacked_bar(data, feature):
    ctab = pd.crosstab(data[feature], data[target], normalize="index")*100
    ctab.plot(kind="bar", stacked=True, figsize=(10,4))
    plt.title(f"{feature}: % Certified vs Denied")
    plt.ylabel("Percentage")
    plt.xticks(rotation=45, ha="right")
    plt.tight_layout()
    plt.show()

for c in [x for x in cat_cols if x.upper() != "CASE_STATUS"]:
    stacked_bar(df_clean, c)

plt.figure(figsize=(8,6))
sns.heatmap(df_clean.select_dtypes(include="number").corr(), cmap="coolwarm", annot=False)
plt.title("Correlation Map (Numeric Features)")
plt.tight_layout(); plt.show()



"""### 7.	Analytical Insights

**Outcome**: Certification rates are highest among applicants with advanced education and high prevailing wages, confirming strong positive influence from wage and education.

In order to identify deeper relationships before modeling, the code defines a helper function to compute certification rates by feature category and explores how certification varies across wage levels and education. Applicants earning above the 75th percentile of prevailing wage are analyzed by education level, region, and experience. The analysis reveals that higher wages, advanced education (master’s and doctorate), and prior job experience strongly increase the likelihood of certification.
"""

# ============================================================
# 7) Analytical Insights
# ============================================================
def cert_rate_by(df, col):
    table = df.groupby(col)[target].mean().rename("cert_rate").to_frame()
    table["count"] = df.groupby(col)[target].size()
    return table.sort_values("cert_rate", ascending=False)

df_clean["wage_bin"] = pd.qcut(df_clean["prevailing_wage_yearly"], q=4, duplicates="drop")
high_cut = df_clean["prevailing_wage_yearly"].quantile(0.75)
df_highpay = df_clean[df_clean["prevailing_wage_yearly"] >= high_cut]

edu_highpay = cert_rate_by(df_highpay, "education_of_employee")
continent_cert = cert_rate_by(df_clean, "continent")
exp_cert = cert_rate_by(df_clean, "has_job_experience")
region_wage = df_clean.groupby("region_of_employment")["prevailing_wage_yearly"].agg(["count","median","mean"]).sort_values("median", ascending=False)
wage_cert = cert_rate_by(df_clean, "wage_bin")
unit_cert = cert_rate_by(df_clean, "unit_of_wage")

print("Education among high-paid roles:")
display(edu_highpay)
print("\nCertification rate by continent:")
display(continent_cert)
print("\nCertification rate by prior experience:")
display(exp_cert)
print("\nPrevailing wage by region:")
display(region_wage)
print("\nCertification by wage quartile:")
display(wage_cert)
print("\nCertification by wage unit:")
display(unit_cert)



"""### 8.	Data Preparation for Modeling

**Outcome**: Dataset split into Train (15,288), Validation (5,096), and Test (5,096) subsets; preprocessing pipeline successfully built with numeric scaling and categorical encoding.

In order to prepare data for unbiased modeling, the dataset is split into 60% training, 20% validation, and 20% test subsets while preserving class proportions (Certified ≈ 67%, Denied ≈ 33%). A preprocessing pipeline is built using scikit-learn’s ColumnTransformer to handle numerical and categorical features—numerical columns are imputed with median values and scaled, while categorical variables are imputed with the most frequent value and one-hot encoded. This ensures consistent transformations across all models and prevents data leakage.
"""

# ============================================================
# 8) Data Preparation for Modeling
# ============================================================
X = df_clean.drop(columns=[target])
y = df_clean[target].astype(int)

X_train, X_temp, y_train, y_temp = train_test_split(X, y, test_size=0.4, stratify=y, random_state=RANDOM_STATE)
X_valid, X_test, y_valid, y_test = train_test_split(X_temp, y_temp, test_size=0.5, stratify=y_temp, random_state=RANDOM_STATE)

num_cols = X_train.select_dtypes(include="number").columns.tolist()
cat_cols = X_train.select_dtypes(exclude="number").columns.tolist()

numeric_tf = Pipeline([("impute", SimpleImputer(strategy="median")), ("scale", StandardScaler())])
categorical_tf = Pipeline([("impute", SimpleImputer(strategy="most_frequent")), ("encode", OneHotEncoder(handle_unknown="ignore"))])

preprocessor = ColumnTransformer([("num", numeric_tf, num_cols), ("cat", categorical_tf, cat_cols)], remainder="drop")
print("Data preprocessed successfully.")



"""### 9.	Metric and Visualization Utilities

**Outcome**: Utility functions were created to calculate model performance metrics including accuracy, precision, recall, F1, ROC-AUC, and specificity, and to visualize confusion matrices for result interpretation.

In order to streamline model evaluation, reusable functions were built to compute key classification metrics and plot confusion matrices. These tools standardize performance tracking across all models and ensure that both certified (approvals) and denied (rejections) outcomes are assessed fairly. The inclusion of specificity allows explicit measurement of how accurately the model identifies denied cases, which is central to the business goal of reducing approval errors.
"""

# ============================================================
# 9) Metric and Visualization Utilities
# ============================================================
def evaluate(model, X, y):
    y_pred = model.predict(X)
    y_prob = model.predict_proba(X)[:,1] if hasattr(model.named_steps["clf"], "predict_proba") else None
    acc = accuracy_score(y, y_pred)
    prec = precision_score(y, y_pred)
    rec = recall_score(y, y_pred)
    f1 = f1_score(y, y_pred)
    auc = roc_auc_score(y, y_prob) if y_prob is not None else np.nan
    cm = confusion_matrix(y, y_pred)
    tn, fp, fn, tp = cm.ravel()
    spec = tn/(tn+fp)
    return acc, prec, rec, f1, auc, spec, cm

def plot_cm(cm, title="Confusion Matrix"):
    plt.figure(figsize=(5,4))
    sns.heatmap(cm, annot=True, fmt="d", cmap="Greens",
                xticklabels=["Pred Denied","Pred Certified"],
                yticklabels=["Actual Denied","Actual Certified"])
    plt.title(title)
    plt.tight_layout(); plt.show()



"""### 10.	Baseline Models on ORIGINAL TRAIN (No Resampling)

**Outcome**: Gradient Boosting achieved the best baseline with approximately 82% F1, 86% recall for certified cases, and 49% denied specificity, establishing strong approval prediction but weaker denial detection.

In order to establish baseline performance, multiple machine learning models were trained on the original, imbalanced data using class weights to account for certification dominance. Models included Logistic Regression, Decision Tree, Random Forest, Bagging, AdaBoost, Gradient Boosting, and XGBoost. Gradient Boosting and XGBoost outperformed others, achieving F1 scores near 0.82, high recall for certified applications, and moderate specificity for denials. This demonstrates strong capability in predicting approvals but the need for better recognition of denials.
"""

# ============================================================
# 10) Baseline Models on ORIGINAL TRAIN (No Resampling)
#     - Class weights emphasize DENIED (0)
#     - Validate on untouched VALID split
# ============================================================
def make_plain(model):
    return Pipeline([("pre", preprocessor), ("clf", model)])

# Compute weights from the train split (robust to class ratio)
classes = np.array([0,1])
weights = compute_class_weight(class_weight="balanced", classes=classes, y=y_train)
class_weights = {0: float(weights[0]), 1: float(weights[1])}
print("Class weights (train-derived):", class_weights)

orig_models = {
    "Logistic (weighted)":   make_plain(LogisticRegression(
        solver="liblinear", class_weight=class_weights, max_iter=3000, random_state=RANDOM_STATE)),
    "Decision Tree (weighted)": make_plain(DecisionTreeClassifier(
        class_weight=class_weights, random_state=RANDOM_STATE)),
    "Random Forest (weighted)": make_plain(RandomForestClassifier(
        n_estimators=400, class_weight=class_weights, n_jobs=-1, random_state=RANDOM_STATE)),
    "Bagging (DT base)":     make_plain(BaggingClassifier(
        n_estimators=300, random_state=RANDOM_STATE)),
    "AdaBoost":              make_plain(AdaBoostClassifier(
        n_estimators=400, random_state=RANDOM_STATE)),
    "Gradient Boosting":     make_plain(GradientBoostingClassifier(
        random_state=RANDOM_STATE)),
}

# Add XGBoost if available
if XGB_OK:
    orig_models["XGBoost"] = make_plain(XGBClassifier(
        random_state=RANDOM_STATE,
        n_estimators=500, learning_rate=0.1, max_depth=4,
        subsample=0.9, colsample_bytree=0.9, eval_metric="logloss",
        tree_method="hist", n_jobs=-1))

rows, fitted_orig = [], {}
for name, pipe in orig_models.items():
    pipe.fit(X_train, y_train)
    acc, prec, rec, f1, auc, spec, cm = evaluate(pipe, X_valid, y_valid)
    rows.append([name, acc, prec, rec, f1, auc, spec])
    fitted_orig[name] = pipe

orig_perf = (pd.DataFrame(rows, columns=["Model","Accuracy","Precision(1)","Recall(1)",
                                        "F1(1)","ROC_AUC","Spec(0)"])
             .sort_values(by=["F1(1)","ROC_AUC"], ascending=False).reset_index(drop=True))

print("=== ORIGINAL TRAIN — Validation Results ===")
display(orig_perf.round(4))



"""### 11.	Model Building — OVERSAMPLED TRAIN (SMOTE)

**Outcome**: Oversampling increased denied detection to around 61% specificity while maintaining certified recall above 80%, showing progress in balancing predictions.

In order to address class imbalance, the Synthetic Minority Oversampling Technique (SMOTE) was applied to the training set, ensuring models received equal exposure to denied cases. The same suite of models was retrained using this balanced dataset. Oversampling improved denied specificity and stabilized F1 performance while slightly reducing overall accuracy. These results confirmed that balancing the dataset helps the model better recognize denials without sacrificing too much approval precision.
"""

# ============================================================
# 11) Model Building — OVERSAMPLED TRAIN (SMOTE)
#     - Resample only inside training fold via Pipeline
#     - Evaluate on untouched VALID set
# ============================================================

from imblearn.pipeline import Pipeline as ImbPipeline
from imblearn.over_sampling import SMOTE

def make_smote(model):
    return ImbPipeline([
        ("pre", preprocessor),
        ("smote", SMOTE(random_state=RANDOM_STATE)),
        ("clf", model)
    ])

smote_models = {
    "Logistic + SMOTE (weighted)": make_smote(LogisticRegression(
        solver="liblinear", class_weight=class_weights, max_iter=3000, random_state=RANDOM_STATE)),
    "Decision Tree + SMOTE (weighted)": make_smote(DecisionTreeClassifier(
        class_weight=class_weights, random_state=RANDOM_STATE)),
    "Random Forest + SMOTE (weighted)": make_smote(RandomForestClassifier(
        n_estimators=400, class_weight=class_weights, n_jobs=-1, random_state=RANDOM_STATE)),
    "Gradient Boosting + SMOTE": make_smote(GradientBoostingClassifier(random_state=RANDOM_STATE)),
    "AdaBoost + SMOTE": make_smote(AdaBoostClassifier(n_estimators=400, random_state=RANDOM_STATE))
}

rows, fitted_smote = [], {}
for name, pipe in smote_models.items():
    pipe.fit(X_train, y_train)
    acc, prec, rec, f1, auc, spec, cm = evaluate(pipe, X_valid, y_valid)
    rows.append([name, acc, prec, rec, f1, auc, spec])
    fitted_smote[name] = pipe

smote_perf = pd.DataFrame(rows, columns=["Model","Accuracy","Precision(1)","Recall(1)",
                                         "F1(1)","ROC_AUC","Spec(0)"])\
                            .sort_values(by=["F1(1)","Spec(0)"], ascending=False)\
                            .reset_index(drop=True)

print("=== SMOTE TRAIN — Validation Results ===")
display(smote_perf.round(4))



"""### 12.	Model Building — UNDERSAMPLED TRAIN

**Outcome**: Undersampling improved denied recognition to approximately 66% specificity but reduced certified recall to around 75%, illustrating the tradeoff between the two outcomes.

In order to further evaluate imbalance handling, the majority class (certified) was reduced to match the number of denied cases, creating a balanced but smaller training dataset. Models trained on this data showed stronger denial prediction but weaker performance for approvals. The reduced dataset also slightly decreased generalization performance, confirming that undersampling improves fairness but requires optimization to avoid overfitting to smaller data.
"""

# ============================================================
# 12) Model Building — UNDERSAMPLED TRAIN
#     - Undersample only inside training fold via Pipeline
#     - Evaluate on untouched VALID set
# ============================================================

from imblearn.under_sampling import RandomUnderSampler

def make_under(model):
    return ImbPipeline([
        ("pre", preprocessor),
        ("under", RandomUnderSampler(random_state=RANDOM_STATE)),
        ("clf", model)
    ])

under_models = {
    "Logistic + Under (weighted)": make_under(LogisticRegression(
        solver="liblinear", class_weight=class_weights, max_iter=3000, random_state=RANDOM_STATE)),
    "Decision Tree + Under (weighted)": make_under(DecisionTreeClassifier(
        class_weight=class_weights, random_state=RANDOM_STATE)),
    "Random Forest + Under (weighted)": make_under(RandomForestClassifier(
        n_estimators=400, class_weight=class_weights, n_jobs=-1, random_state=RANDOM_STATE)),
    "Gradient Boosting + Under": make_under(GradientBoostingClassifier(random_state=RANDOM_STATE)),
    "AdaBoost + Under": make_under(AdaBoostClassifier(n_estimators=400, random_state=RANDOM_STATE))
}

rows, fitted_under = [], {}
for name, pipe in under_models.items():
    pipe.fit(X_train, y_train)
    acc, prec, rec, f1, auc, spec, cm = evaluate(pipe, X_valid, y_valid)
    rows.append([name, acc, prec, rec, f1, auc, spec])
    fitted_under[name] = pipe

under_perf = pd.DataFrame(rows, columns=["Model","Accuracy","Precision(1)","Recall(1)",
                                         "F1(1)","ROC_AUC","Spec(0)"])\
                            .sort_values(by=["F1(1)","Spec(0)"], ascending=False)\
                            .reset_index(drop=True)

print("=== UNDERSAMPLED TRAIN — Validation Results ===")
display(under_perf.round(4))



"""### 13.	Select Top Candidates Across All Methods

**Outcome**: Gradient Boosting, Random Forest (weighted), and XGBoost were selected as top performers for tuning, based on the highest F1 and balanced recall-specificity profiles.

In order to optimize computational effort and focus on the most promising approaches, model results from the original, oversampled, and undersampled sets were compared. The top three models were selected for hyperparameter tuning based on their F1 scores (balance of precision and recall) and their ability to detect both certified and denied outcomes effectively.
"""

# ============================================================
# 13) Select Top Candidates Across All Methods
# ============================================================

def pick_top(df, k=3):
    return df.sort_values(by=["F1(1)","Spec(0)"], ascending=False).head(k)["Model"].tolist()

candidates = list(dict.fromkeys(
    pick_top(orig_perf, 3) + pick_top(smote_perf, 3) + pick_top(under_perf, 3)
))

print("=== Top Candidates for Tuning ===")
for c in candidates:
    print(" -", c)



"""### 14.	Hyperparameter Tuning (CV + Randomized Search)

**Outcome**: Optimized Gradient Boosting achieved an F1 score of about 0.83 with balanced recall and specificity, confirming it as the final model.

In order to maximize model performance, RandomizedSearchCV with stratified cross-validation was applied to fine-tune the leading models. Key parameters like learning rate, tree depth, and estimator count were adjusted. Gradient Boosting consistently outperformed others after tuning, achieving the best mix of recall, F1, and denial specificity. This step confirmed its suitability for operational deployment.
"""

# ============================================================
# 14) Hyperparameter Tuning (CV + Randomized Search)
#     - Focus: F1 optimization while tracking DENIED specificity
# ============================================================

from sklearn.model_selection import RandomizedSearchCV, StratifiedKFold
from scipy.stats import loguniform, randint

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)

def build_model_by_name(name):
    if "Gradient Boosting" in name:
        base = GradientBoostingClassifier(random_state=RANDOM_STATE)
        pipe = make_plain(base)
        params = {
            "clf__n_estimators": [200, 300, 500],
            "clf__learning_rate": [0.05, 0.1, 0.2],
            "clf__max_depth": [2, 3, 4],
            "clf__subsample": [0.8, 1.0]
        }
        return pipe, params

    if "Random Forest" in name:
        base = RandomForestClassifier(
            n_estimators=400, class_weight=class_weights, n_jobs=-1, random_state=RANDOM_STATE)
        pipe = make_plain(base)
        params = {
            "clf__max_depth": [None, 10, 20],
            "clf__min_samples_split": [2, 5, 10],
            "clf__min_samples_leaf": [1, 2, 4],
            "clf__max_features": ["sqrt", "log2"]
        }
        return pipe, params

    if "Logistic" in name:
        base = LogisticRegression(
            solver="liblinear", class_weight=class_weights, max_iter=4000, random_state=RANDOM_STATE)
        pipe = make_plain(base)
        params = {"clf__C": loguniform(1e-3, 1e2), "clf__penalty": ["l1", "l2"]}
        return pipe, params

    return None

tuned_models, results = {}, []
for name in candidates[:3]:
    built = build_model_by_name(name)
    if built is None:
        continue
    pipe, grid = built
    tuner = RandomizedSearchCV(
        estimator=pipe,
        param_distributions=grid,
        n_iter=20,
        scoring="f1",
        cv=cv,
        n_jobs=-1,
        random_state=RANDOM_STATE,
        refit=True,
        verbose=0
    )
    tuner.fit(X_train, y_train)
    tuned_models[name] = tuner.best_estimator_

    acc, prec, rec, f1, auc, spec, cm = evaluate(tuner.best_estimator_, X_valid, y_valid)
    results.append([name, tuner.best_params_, acc, prec, rec, f1, auc, spec])

tuned_df = pd.DataFrame(results, columns=[
    "Model","Best Params","Accuracy","Precision(1)","Recall(1)","F1(1)","ROC_AUC","Spec(0)"
]).sort_values(by=["F1(1)","Spec(0)"], ascending=False).reset_index(drop=True)

print("=== Tuned Models — Validation ===")
display(tuned_df)

best_name = tuned_df.iloc[0]["Model"]
best_model = tuned_models[best_name]
print("Selected tuned model:", best_name)



"""### 15.	Threshold Optimization

**Outcome**: Optimal threshold tuning achieved approximately 74% accuracy, 86% recall for certified applications, and 49% denied specificity under the default setting, later improved through calibration.

In order to better balance the trade-off between identifying approvals and denials, probability thresholds were adjusted. This allowed the model to shift sensitivity toward detecting denials while retaining strong approval prediction. The chosen threshold offered a better operational balance, addressing business goals of fairness and reliability.
"""

# ============================================================
# 15) Threshold Optimization
#     - Balance Certified Recall vs Denied Specificity
# ============================================================

def proba_positive(pipe, X):
    clf = pipe.named_steps["clf"]
    if hasattr(clf, "predict_proba"):
        return pipe.predict_proba(X)[:, 1]
    raise RuntimeError("This model does not support predict_proba().")

val_prob = proba_positive(best_model, X_valid)

thresholds = np.linspace(0.2, 0.8, 121)
records = []
for t in thresholds:
    y_pred = (val_prob >= t).astype(int)
    acc = accuracy_score(y_valid, y_pred)
    prec = precision_score(y_valid, y_pred, zero_division=0)
    rec = recall_score(y_valid, y_pred, zero_division=0)
    f1 = f1_score(y_valid, y_pred, zero_division=0)
    tn, fp, fn, tp = confusion_matrix(y_valid, y_pred, labels=[0,1]).ravel()
    spec = tn / (tn + fp) if (tn + fp) else 0
    score = 0.6*f1 + 0.4*spec
    records.append([t, acc, prec, rec, f1, spec, score])

thresh_df = pd.DataFrame(records, columns=["Threshold","Accuracy","Precision(1)","Recall(1)",
                                           "F1(1)","Spec(0)","Score"])\
                        .sort_values(by="Score", ascending=False).reset_index(drop=True)

best_thr = float(thresh_df.iloc[0]["Threshold"])
print("Best operating threshold:", round(best_thr, 3))
display(thresh_df.head(10))



"""### 16.	Validation Threshold Sweep — Multi-Threshold Confusion Metrics

**Outcome**: At the balanced threshold, certified and denied predictions were both correctly identified around 70% of the time, showing equal performance across classes.

In order to understand performance trade-offs at different probability thresholds, a series of confusion matrices and metric calculations were generated across a range of cutoff values. This analysis highlighted how increasing the threshold improved denied specificity while reducing certified recall, and vice versa. The balanced threshold was identified as the point where both outcomes achieved similar accuracy, supporting fairer model decisions.
"""

# ============================================================
# 16) Validation Threshold Sweep — Multi-Threshold Confusion Metrics
#     - Computes metrics at multiple thresholds (incl. default 0.50)
#     - Finds: best-F1 threshold, "even" (recall vs specificity) threshold
#     - Shows confusion stats for a curated set of thresholds
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import (
    accuracy_score, precision_score, recall_score, f1_score, roc_auc_score,
    confusion_matrix
)

# In case not already defined earlier
def proba_positive(est, X):
    clf = est.named_steps["clf"]
    if hasattr(clf, "predict_proba"):
        return est.predict_proba(X)[:, 1]
    raise RuntimeError("Estimator lacks predict_proba for threshold tuning.")

def eval_threshold(y_true, y_prob, thr):
    y_pred = (y_prob >= thr).astype(int)
    acc  = accuracy_score(y_true, y_pred)
    prec = precision_score(y_true, y_pred, pos_label=1, zero_division=0)
    rec  = recall_score(y_true, y_pred, pos_label=1, zero_division=0)   # Certified recall
    f1c  = f1_score(y_true, y_pred, pos_label=1, zero_division=0)
    cm   = confusion_matrix(y_true, y_pred, labels=[0, 1])
    tn, fp, fn, tp = cm.ravel()
    spec = tn/(tn+fp) if (tn+fp) else 0.0                               # Denied specificity
    cert_ok = tp/(tp+fn) if (tp+fn) else 0.0
    den_ok  = tn/(tn+fp) if (tn+fp) else 0.0
    return {
        "threshold": thr, "accuracy": acc, "precision(1)": prec, "recall(1)": rec,
        "f1(1)": f1c, "spec(0)": spec, "cert_ok%": cert_ok, "denied_ok%": den_ok,
        "tn": tn, "fp": fp, "fn": fn, "tp": tp, "cm": cm
    }

def summarize_thresholds(y_true, y_prob, thr_list, title_prefix="Validation", show_heatmaps=False):
    rows = []
    for t in thr_list:
        rows.append(eval_threshold(y_true, y_prob, t))
    df = pd.DataFrame(rows).drop(columns=["cm"])
    df_sorted = df.sort_values(by=["f1(1)", "spec(0)"], ascending=False).reset_index(drop=True)
    display(df_sorted.round(4))

    if show_heatmaps:
        for r in rows:
            cm = r["cm"]
            cm_df = pd.DataFrame(cm,
                                 index=["Actual DENIED (0)", "Actual CERTIFIED (1)"],
                                 columns=["Pred DENIED (0)", "Pred CERTIFIED (1)"])
            plt.figure(figsize=(5.2,4.2))
            sns.heatmap(cm_df, annot=True, fmt="d", cmap="Purples", cbar=True, linewidths=.5)
            plt.title(f"{title_prefix} — Confusion Matrix @ Thr {r['threshold']:.3f}")
            plt.tight_layout(); plt.show()
    return df_sorted

# ---- Run on VALIDATION ----
val_prob = proba_positive(best_model, X_valid)

# Fine sweep to find best-F1 and best "even" (harmonic mean of recall & specificity)
grid = np.linspace(0.20, 0.80, 121)
scan = [eval_threshold(y_valid, val_prob, t) for t in grid]
scan_df = pd.DataFrame(scan).drop(columns=["cm"])
scan_df["even_hmean"] = (2*scan_df["recall(1)"]*scan_df["spec(0)"]) / (
    (scan_df["recall(1)"] + scan_df["spec(0)"]).replace(0, np.nan)
)

VAL_best_f1_threshold   = float(scan_df.loc[scan_df["f1(1)"].idxmax(), "threshold"])
VAL_even_threshold      = float(scan_df.loc[scan_df["even_hmean"].idxmax(), "threshold"])
VAL_default_threshold   = 0.50

print(f"Validation best-F1 threshold   : {VAL_best_f1_threshold:.3f}")
print(f"Validation even-balance thr    : {VAL_even_threshold:.3f}")
print(f"Validation default threshold   : {VAL_default_threshold:.3f}")

# Curated set of thresholds to inspect on VALIDATION
val_thresholds_to_view = sorted({
    VAL_default_threshold,
    round(VAL_best_f1_threshold, 3),
    round(VAL_even_threshold, 3),
    0.60, 0.65, 0.68, 0.70
})

print("\n=== Validation — Metrics at Selected Thresholds ===")
val_table = summarize_thresholds(y_valid, val_prob, val_thresholds_to_view,
                                 title_prefix="Validation", show_heatmaps=False)



"""### 17.	Threshold Sweep with Key Operating Points

**Outcome**: On the test data, the balanced threshold maintained 70.3% certified and 69.7% denied correctness, confirming stable generalization.

In order to validate real-world reliability, the same threshold comparisons were applied to the untouched test set. The model preserved balanced accuracy across both classes, demonstrating consistency under new data. This confirmed the tuned model’s robustness and its potential for real-world deployment in the visa decision process.
"""

# ============================================================
# 17) Threshold Sweep with Key Operating Points
#     - Best F1 (Certified)
#     - Balanced predicted counts (≈ equal Certified vs Denied predictions)
#     - Parity of correctness (Certified recall ≈ Denied specificity)
#     - Always include default 0.50
#     Prints counts, % correct/incorrect, and confusion matrices
# ============================================================

import numpy as np
import pandas as pd
from sklearn.metrics import (
    accuracy_score, precision_score, recall_score, f1_score, roc_auc_score, confusion_matrix
)

# --- Choose dataset to inspect ---
DATASET = "test"   # "valid" or "test"
if DATASET == "valid":
    X_eval, y_eval = X_valid, y_valid
else:
    X_eval, y_eval = X_test, y_test

# --- Probabilities from the selected model ---
if not hasattr(best_model.named_steps["clf"], "predict_proba"):
    raise RuntimeError("The selected model does not support predict_proba; threshold tuning requires probabilities.")

y_prob = best_model.predict_proba(X_eval)[:, 1]  # P(CERTIFIED)
y_true = y_eval.values if hasattr(y_eval, "values") else y_eval

# --- Helpers ---
def eval_at_threshold(y_true, y_prob, thr):
    y_pred = (y_prob >= thr).astype(int)
    cm = confusion_matrix(y_true, y_pred, labels=[0, 1])
    tn, fp, fn, tp = cm.ravel()
    total = tn + fp + fn + tp

    # Core metrics
    acc  = accuracy_score(y_true, y_pred)
    prec = precision_score(y_true, y_pred, pos_label=1, zero_division=0)
    rec1 = recall_score(y_true, y_pred, pos_label=1, zero_division=0)        # Certified recall (TPR for class 1)
    f1c  = f1_score(y_true, y_pred, pos_label=1, zero_division=0)
    spec0 = tn / (tn + fp) if (tn + fp) else 0.0                              # Denied specificity (TNR for class 0)

    # Predicted counts & shares
    pred_1 = (y_pred == 1).sum()
    pred_0 = (y_pred == 0).sum()
    share_1 = pred_1 / total
    share_0 = pred_0 / total

    # Class-wise correctness
    tot_denied = tn + fp
    tot_cert   = tp + fn
    denied_correct_pct   = (tn / tot_denied) * 100 if tot_denied else 0.0
    denied_incorrect_pct = (fp / tot_denied) * 100 if tot_denied else 0.0
    cert_correct_pct     = (tp / tot_cert) * 100 if tot_cert else 0.0
    cert_incorrect_pct   = (fn / tot_cert) * 100 if tot_cert else 0.0

    return {
        "thr": thr,
        "acc": acc,
        "prec1": prec,
        "rec1": rec1,
        "spec0": spec0,
        "f1c": f1c,
        "pred_0": int(pred_0), "pred_1": int(pred_1),
        "share_0": share_0, "share_1": share_1,
        "tn": int(tn), "fp": int(fp), "fn": int(fn), "tp": int(tp),
        "denied_correct_pct": denied_correct_pct,
        "denied_incorrect_pct": denied_incorrect_pct,
        "cert_correct_pct": cert_correct_pct,
        "cert_incorrect_pct": cert_incorrect_pct,
        "cm": cm
    }

# --- Sweep thresholds ---
grid = np.linspace(0.20, 0.80, 121)
rows = [eval_at_threshold(y_true, y_prob, t) for t in grid]
df = pd.DataFrame(rows)

# --- Pick operating points ---
# 1) Default 0.50
pt_default = df.iloc[(df['thr']-0.50).abs().argmin()]

# 2) Best F1 (Certified) — strongest single-score balance for the positive class
pt_best_f1 = df.iloc[df['f1c'].values.argmax()]

# 3) Balanced predicted counts — roughly equal number of predicted 0 and 1
pt_balanced_counts = df.iloc[(df['pred_0'] - df['pred_1']).abs().argmin()]

# 4) Parity of correctness — Certified recall ≈ Denied specificity
pt_parity = df.iloc[(df['rec1'] - df['spec0']).abs().argmin()]

# --- Collate for display ---
def short_row(pt, label):
    return {
        "Point": label,
        "Threshold": round(float(pt['thr']), 3),
        "Accuracy": round(float(pt['acc']), 4),
        "Precision(1)": round(float(pt['prec1']), 4),
        "Recall(1) (Certified)": round(float(pt['rec1']), 4),
        "Specificity(0) (Denied)": round(float(pt['spec0']), 4),
        "F1(1)": round(float(pt['f1c']), 4),
        "Pred Certified (count)": int(pt['pred_1']),
        "Pred Denied (count)": int(pt['pred_0']),
        "Certified correct %": f"{pt['cert_correct_pct']:.2f}%",
        "Denied correct %": f"{pt['denied_correct_pct']:.2f}%"
    }

summary = pd.DataFrame([
    short_row(pt_default, "Default (0.50)"),
    short_row(pt_best_f1, "Best F1 (Certified)"),
    short_row(pt_balanced_counts, "Balanced Predicted Counts"),
    short_row(pt_parity, "Parity of Correctness (Recall≈Specificity)")
])

print(f"=== Threshold Summary on {DATASET.upper()} ===")
display(summary)

# --- Pretty print confusion matrices & correctness breakdown ---
def print_block(pt, label):
    cm = pt['cm']
    tn, fp, fn, tp = pt['tn'], pt['fp'], pt['fn'], pt['tp']
    print(f"\n--- {label} @ threshold={pt['thr']:.3f} ---")
    print(pd.DataFrame(cm,
                       index=["Actual DENIED (0)", "Actual CERTIFIED (1)"],
                       columns=["Pred DENIED (0)", "Pred CERTIFIED (1)"]))
    print(f"\nCounts:")
    print(f"  Predicted CERTIFIED (1): {pt['pred_1']}")
    print(f"  Predicted DENIED (0)   : {pt['pred_0']}")

    print(f"\nCorrectness (% of each actual class):")
    print(f"  Certified correctly predicted: {pt['cert_correct_pct']:.2f}%   (Incorrect: {pt['cert_incorrect_pct']:.2f}%)")
    print(f"  Denied correctly predicted   : {pt['denied_correct_pct']:.2f}%   (Incorrect: {pt['denied_incorrect_pct']:.2f}%)")

    print(f"\nMetrics:")
    print(f"  Accuracy: {pt['acc']:.4f} | Precision(1): {pt['prec1']:.4f} | Recall(1): {pt['rec1']:.4f} | Specificity(0): {pt['spec0']:.4f} | F1(1): {pt['f1c']:.4f}")

print_block(pt_default, "Default (0.50)")
print_block(pt_best_f1, "Best F1 (Certified)")
print_block(pt_balanced_counts, "Balanced Predicted Counts")
print_block(pt_parity, "Parity of Correctness (Recall≈Specificity)")



"""### 18.	Business Alignment: Outcomes vs. Problem Statement

**Outcome**: The final model accurately predicts about 70% of both certified and denied applications, identifies key drivers of approval, and supports faster, fairer visa reviews.

In order to align findings with the original business problem, results were mapped to operational goals. The model achieves balanced predictive accuracy, identifying education, experience, prevailing wage, and employment region as major influences on certification outcomes. These insights allow EasyVisa to streamline the visa evaluation process, ensuring decisions are more consistent, data-driven, and efficient while reducing manual workload and potential bias.
"""

# ============================================================
# 18) Business Alignment: Outcomes vs. Problem Statement
#     - Quantify how we facilitate approvals AND catch denials
#     - Report results at: default 0.50, Best-F1(thr on VALID), and Even(Recall≈Specificity)
#     - Summarize in business language + technical metrics
#     - Show top drivers (features) used for recommendations
# ============================================================

import numpy as np
import pandas as pd
from sklearn.metrics import (
    accuracy_score, precision_score, recall_score, f1_score, roc_auc_score, confusion_matrix
)

# -------- Helpers --------
def proba_positive(est, X):
    """P(CERTIFIED=1)."""
    clf = est.named_steps["clf"]
    if hasattr(clf, "predict_proba"):
        return est.predict_proba(X)[:, 1]
    raise RuntimeError("Selected model does not provide predict_proba; required for threshold analysis.")

def evaluate_at_threshold(y_true, y_prob, thr):
    y_pred = (y_prob >= thr).astype(int)
    cm = confusion_matrix(y_true, y_pred, labels=[0, 1])
    tn, fp, fn, tp = cm.ravel()
    acc  = accuracy_score(y_true, y_pred)
    prec = precision_score(y_true, y_pred, pos_label=1, zero_division=0)
    rec1 = recall_score(y_true, y_pred, pos_label=1, zero_division=0)      # Certified recall
    f1c  = f1_score(y_true, y_pred, pos_label=1, zero_division=0)
    spec0 = tn / (tn + fp) if (tn + fp) else 0.0                            # Denied specificity
    auc  = roc_auc_score(y_true, y_prob)

    total_denied = tn + fp
    total_cert   = tp + fn
    out = {
        "Threshold": thr,
        "Accuracy": acc,
        "Precision(1)": prec,
        "Recall(1)_Certified": rec1,
        "Specificity(0)_Denied": spec0,
        "F1(1)": f1c,
        "ROC_AUC": auc,
        "TP": int(tp), "FP": int(fp), "TN": int(tn), "FN": int(fn),
        "Total_Denied": int(total_denied), "Total_Certified": int(total_cert),
        "Cert_Correct_%": (tp / total_cert * 100) if total_cert else 0.0,
        "Denied_Correct_%": (tn / total_denied * 100) if total_denied else 0.0
    }
    return out, cm

def best_f1_threshold_from_valid(model, X_valid, y_valid):
    """Pick threshold that maximizes F1 on VALID; fallback to 0.50 if anything goes wrong."""
    try:
        vprob = proba_positive(model, X_valid)
        grid = np.linspace(0.20, 0.80, 121)
        best_t, best_f1 = 0.50, -1
        for t in grid:
            ypv = (vprob >= t).astype(int)
            f1c = f1_score(y_valid, ypv, pos_label=1, zero_division=0)
            if f1c > best_f1:
                best_f1, best_t = f1c, t
        return float(best_t)
    except Exception:
        return 0.50

def even_threshold_from_valid(model, X_valid, y_valid):
    """Pick threshold that balances Certified recall and Denied specificity on VALID."""
    try:
        vprob = proba_positive(model, X_valid)
        grid = np.linspace(0.20, 0.80, 121)
        best_t, best_hm = 0.50, -1
        for t in grid:
            ypv = (vprob >= t).astype(int)
            cm = confusion_matrix(y_valid, ypv, labels=[0,1])
            tn, fp, fn, tp = cm.ravel()
            rec1 = tp / (tp + fn) if (tp + fn) else 0.0
            spec0 = tn / (tn + fp) if (tn + fp) else 0.0
            hm = (2 * rec1 * spec0) / (rec1 + spec0) if (rec1 + spec0) else 0.0
            if hm > best_hm:
                best_hm, best_t = hm, t
        return float(best_t)
    except Exception:
        return 0.50

# -------- Determine thresholds (use existing ones if defined earlier) --------
try:
    THR_BEST_F1 = VAL_best_f1_threshold  # from earlier cell if present
except NameError:
    THR_BEST_F1 = best_f1_threshold_from_valid(best_model, X_valid, y_valid)

try:
    THR_EVEN = EVEN_thr  # from earlier cell if present
except NameError:
    THR_EVEN = even_threshold_from_valid(best_model, X_valid, y_valid)

THR_DEFAULT = 0.50

# -------- Evaluate on TEST (business acceptance requires held-out test) --------
test_prob = proba_positive(best_model, X_test)

res_default, cm_default = evaluate_at_threshold(y_test, test_prob, THR_DEFAULT)
res_bestf1, cm_bestf1   = evaluate_at_threshold(y_test, test_prob, THR_BEST_F1)
res_even, cm_even       = evaluate_at_threshold(y_test, test_prob, THR_EVEN)

summary = pd.DataFrame([
    {"Point": "Default (0.50)", **{k: v for k, v in res_default.items() if k not in ["TP","FP","TN","FN"]}},
    {"Point": "Best-F1 (Valid-chosen)", **{k: v for k, v in res_bestf1.items() if k not in ["TP","FP","TN","FN"]}},
    {"Point": "Even Recall/Specificity", **{k: v for k, v in res_even.items() if k not in ["TP","FP","TN","FN"]}},
]).round(4)

print("=== Business Alignment Summary (TEST) ===")
display(summary)

def print_confusion_block(name, res, cm):
    print(f"\n--- {name} @ threshold={res['Threshold']:.3f} ---")
    print(pd.DataFrame(cm,
                       index=["Actual DENIED (0)", "Actual CERTIFIED (1)"],
                       columns=["Pred DENIED (0)", "Pred CERTIFIED (1)"]))
    print(f"Certified correctly predicted: {res['Cert_Correct_%']:.2f}% "
          f"({res['TP']} / {res['Total_Certified']})")
    print(f"Denied    correctly predicted: {res['Denied_Correct_%']:.2f}% "
          f"({res['TN']} / {res['Total_Denied']})")
    print(f"Accuracy={res['Accuracy']:.4f} | Precision(1)={res['Precision(1)']:.4f} | "
          f"Recall(1)={res['Recall(1)_Certified']:.4f} | Specificity(0)={res['Specificity(0)_Denied']:.4f} | "
          f"F1(1)={res['F1(1)']:.4f} | ROC_AUC={res['ROC_AUC']:.4f}")

print_confusion_block("Default (0.50)", res_default, cm_default)
print_confusion_block("Best-F1 (Valid-chosen)", res_bestf1, cm_bestf1)
print_confusion_block("Even Recall/Specificity", res_even, cm_even)

# -------- Drivers of recommendation (top features) --------
# Try to extract model-native importances; fallback to permutation on VALID if needed.
def get_feature_names_from_ct(ct):
    """Get transformed feature names from ColumnTransformer."""
    feat_names = []
    for name, trans, cols in ct.transformers_:
        if name == 'remainder' and trans == 'drop':
            continue
        if hasattr(trans, 'get_feature_names_out'):
            # Pipelines: last step may be OHE
            try:
                fn = trans.get_feature_names_out(cols)
            except:
                fn = cols
        else:
            fn = cols
        feat_names.extend(list(fn))
    return np.array(feat_names, dtype=object)

X_names = get_feature_names_from_ct(best_model.named_steps["pre"])

importances = None
clf = best_model.named_steps["clf"]
if hasattr(clf, "feature_importances_"):
    importances = clf.feature_importances_
elif hasattr(clf, "coef_"):
    # Use absolute coefficients (multi-dim safety)
    coef = np.ravel(clf.coef_)
    importances = np.abs(coef)

top_feat_df = None
if importances is not None and len(importances) == len(X_names):
    idx = np.argsort(importances)[::-1][:15]
    top_feat_df = pd.DataFrame({
        "Feature": X_names[idx],
        "Importance": importances[idx]
    })
    print("\nTop drivers (model-native importance/coefficient):")
    display(top_feat_df)
else:
    # Lightweight permutation importance on VALID to avoid leakage
    try:
        from sklearn.inspection import permutation_importance
        val_imp = permutation_importance(best_model, X_valid, y_valid, n_repeats=5, random_state=42, scoring="f1")
        idx = np.argsort(val_imp.importances_mean)[::-1][:15]
        top_feat_df = pd.DataFrame({
            "Feature": X_names[idx],
            "Importance": val_imp.importances_mean[idx]
        })
        print("\nTop drivers (permutation importance on VALID):")
        display(top_feat_df)
    except Exception as e:
        print("\nSkipping feature importance (not supported for this model).", e)

# -------- Executive-readable mapping to the objective --------
def pct(x):
    return f"{100*x:.1f}%"

# Pick the most “balanced” operating point for narrative (Even or Best-F1 depending on context)
narr = res_even if abs(res_even["Recall(1)_Certified"] - res_even["Specificity(0)_Denied"]) <= \
                  abs(res_bestf1["Recall(1)_Certified"] - res_bestf1["Specificity(0)_Denied"]) else res_bestf1
narr_label = "Even Recall/Specificity" if narr is res_even else "Best-F1"

print("\n=== Business Narrative (TEST) ===")
print(
    f"At the '{narr_label}' operating point (threshold={narr['Threshold']:.3f}), the model correctly identifies "
    f"{narr['Cert_Correct_%']:.1f}% of certified cases and {narr['Denied_Correct_%']:.1f}% of denied cases, "
    f"with overall accuracy {pct(narr['Accuracy'])}, F1 for certified {pct(narr['F1(1)'])}, and ROC-AUC {pct(narr['ROC_AUC'])}."
)
print(
    "This directly supports the objective to facilitate approvals (high true-positive rate on certifications) "
    "while protecting U.S. workers by screening out weak applications (strong true-negative rate on denials)."
)

if top_feat_df is not None:
    top_txt = ", ".join(top_feat_df['Feature'].head(5).tolist())
    print(f"Key decision drivers include: {top_txt}. These features inform which profiles are more likely to be certified or denied.")

print("\nInterpretation for stakeholders:")
print(
    "- For operations: the chosen threshold sets how aggressive the shortlisting is. Moving it higher increases the share of predicted denials (more conservative), "
    "moving it lower increases predicted certifications (more throughput)."
)
print(
    "- For compliance and risk: the Even setting balances catching denials with approving qualified applicants, which is easier to justify in audits."
)



"""### Comprehensive Summary

In order to address the U.S. Office of Foreign Labor Certification’s challenge of processing an increasing volume of visa applications efficiently and fairly, this analysis applies machine learning to predict visa approval outcomes using detailed applicant and employer data. The dataset includes 25,480 visa applications, covering variables such as education, experience, prevailing wage, region of employment, and company characteristics. The target variable, case_status, is binary — Certified (1) or Denied (0) — defining this as a supervised classification problem.

After thorough data cleaning and preparation, all categorical values were standardized, inconsistent entries corrected, and wages annualized to ensure comparability. Exploratory data analysis revealed strong relationships between education, experience, and wage with visa outcomes: applicants with higher education levels and prior job experience were more likely to be certified, while those with only high school education or requiring job training faced higher denial rates. Regionally, Northeast and West regions had higher average wages and slightly higher certification rates.

To ensure unbiased model evaluation, the data was split into training (60%), validation (20%), and test (20%) sets. Multiple classification models were tested, including Logistic Regression, Decision Tree, Random Forest, Bagging, AdaBoost, Gradient Boosting, and XGBoost. Class imbalance between certified (67%) and denied (33%) cases was addressed through SMOTE oversampling, undersampling, and class weighting.

Among all models, Gradient Boosting delivered the best results. Initial validation showed strong performance for certified predictions but weaker recognition of denied applications. Through hyperparameter tuning and threshold optimization, the model was refined to balance both classes effectively.

At the chosen balanced operating threshold (0.68), the final model achieved:
	•	Overall accuracy: 70.1%
	•	Certified correctly predicted: 70.35%
	•	Denied correctly predicted: 69.70%
	•	F1 score (Certified): 0.76
	•	ROC-AUC: 0.77

The F1 score of 0.76 indicates a strong overall balance between precision (how many predicted approvals were correct) and recall (how many actual approvals were detected). In this context, that tradeoff is valuable because the model’s purpose is not simply to maximize accuracy, but to maintain fairness—avoiding excessive approvals of weak applications while still identifying qualified candidates accurately.

Feature importance analysis identified the most influential predictors:
	1.	Education level (High School, Bachelor’s, Master’s, Doctorate) — higher education correlates with certification.
	2.	Job experience — applicants with prior experience are more likely to be approved.
	3.	Prevailing wage — higher wages increase the likelihood of certification.
	4.	Continent and region of employment — applications from Europe, North America, and the West region tend to have stronger approval rates.

From an operational standpoint, this balanced model allows OFLC to screen applications fairly, catching weak or noncompliant cases while maintaining high accuracy for strong candidates. The Even Recall/Specificity threshold (0.68) ensures that both certified and denied cases are identified at roughly 70% correctness, a result that aligns directly with business objectives of fairness, transparency, and regulatory compliance.

Outcome: The final Gradient Boosting model provides a fair and data-backed solution that correctly identifies 70% of both visa certifications and denials, with an F1 score of 0.76 and ROC-AUC of 0.77. The F1 score is considered strong in this application, as it reflects the model’s ability to balance false approvals and false denials—supporting consistent, compliant, and efficient visa decision-making.
"""

