Explainable Stacking Ensemble for Software Defect Prediction¶

1. Setup and configuration¶

A single fixed random seed is used everywhere (splits, SMOTE, every stochastic model) so the whole notebook is reproducible. Change DATASETS here if your CSV files live elsewhere.

In [2]:
import os, warnings, json, time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import wilcoxon

from sklearn.model_selection import StratifiedKFold, cross_validate
from sklearn.preprocessing import MinMaxScaler
from sklearn.impute import SimpleImputer
from sklearn.feature_selection import RFE
from sklearn.ensemble import RandomForestClassifier, StackingClassifier
from sklearn.svm import SVC
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (accuracy_score, precision_score, recall_score,
                             f1_score, roc_auc_score, matthews_corrcoef)

from imblearn.over_sampling import SMOTE
from imblearn.pipeline import Pipeline as ImbPipeline   # NOTE: imblearn Pipeline, not sklearn's

from xgboost import XGBClassifier

warnings.filterwarnings("ignore")
sns.set_style("whitegrid")
plt.rcParams["figure.dpi"] = 120
plt.rcParams["savefig.dpi"] = 200
plt.rcParams["font.size"] = 11

RANDOM_STATE = 42
N_FOLDS = 10
N_FEATURES = 15          # number of features RFE keeps; tune in Section 8 if desired
OUTDIR = "outputs"
os.makedirs(OUTDIR, exist_ok=True)

# dataset file -> label column name. The four core datasets from the proposal, plus KC2.
DATASETS = {
    "CM1": ("D:/thesis/nabil/report/code/data/cm1.csv", "defects"),
    "KC1": ("D:/thesis/nabil/report/code/data/kc1.csv", "defects"),
    "PC1": ("D:/thesis/nabil/report/code/data/pc1.csv", "defects"),
    "JM1": ("D:/thesis/nabil/report/code/data/jm1.csv", "defects"),
}
CORE = ["CM1", "KC1", "PC1", "JM1"]   # headline results reported on these four
np.random.seed(RANDOM_STATE)
print("Setup complete. Output directory:", os.path.abspath(OUTDIR))

Setup complete. Output directory: d:\thesis\nabil\report\code\outputs

2. Data loading, label normalisation and validation¶

In [3]:
def load_dataset(path, label_col):
    df = pd.read_csv(path)
    y_raw = df[label_col]
    # normalise a variety of encodings to 1/0
    def to_binary(v):
        s = str(v).strip().lower()
        return 1 if s in ("true", "yes", "y", "1") else 0
    y = y_raw.map(to_binary).astype(int)
    X = df.drop(columns=[label_col])
    # coerce features to numeric (any stray '?' becomes NaN and is imputed later)
    X = X.apply(pd.to_numeric, errors="coerce")
    return X, y

summ = []
data = {}
for name, (path, label) in DATASETS.items():
    X, y = load_dataset(path, label)
    data[name] = (X, y)
    pos = int(y.sum()); n = len(y)
    summ.append({
        "Dataset": name, "Modules": n, "Features": X.shape[1],
        "Defective": pos, "Clean": n - pos,
        "Defect ratio (%)": round(pos / n * 100, 1),
        "Missing cells": int(X.isnull().sum().sum()),
        "Duplicate rows": int(X.duplicated().sum()),
    })

summary_df = pd.DataFrame(summ).set_index("Dataset")
summary_df.to_csv(f"{OUTDIR}/table_dataset_summary.csv")
summary_df
Out[3]:
Modules Features Defective Clean Defect ratio (%) Missing cells Duplicate rows
Dataset
CM1 498 21 49 449 9.8 0 57
KC1 2109 21 326 1783 15.5 0 917
PC1 1109 21 77 1032 6.9 0 161
JM1 13204 21 2103 11101 15.9 0 4384
In [4]:
DROP_DUPLICATES = False   # set True to run a duplicate-free sensitivity check

if DROP_DUPLICATES:
    for name in list(data.keys()):
        X, y = data[name]
        mask = ~X.duplicated().values
        data[name] = (X[mask].reset_index(drop=True), y[mask].reset_index(drop=True))
    print("Duplicates dropped.")
else:
    print("Duplicates kept (default).")

Duplicates kept (default).

3. Exploratory data analysis¶

Three EDA outputs, each saved to outputs/:

  • class-balance bar chart across datasets (motivates SMOTE),
  • a feature-correlation heatmap on the largest dataset (motivates RFE),
  • metric distribution summary.
In [5]:
# 3a. Class balance across datasets
fig, ax = plt.subplots(figsize=(8, 4.5))
ratios = summary_df.loc[list(DATASETS.keys()), "Defect ratio (%)"]
bars = ax.bar(ratios.index, ratios.values, color="#c2185b", edgecolor="black")
ax.set_ylabel("Defective modules (%)")
ax.set_title("Class imbalance across NASA PROMISE datasets")
for b, v in zip(bars, ratios.values):
    ax.text(b.get_x() + b.get_width()/2, v + 0.3, f"{v}%", ha="center", fontsize=10)
plt.tight_layout()
plt.savefig(f"{OUTDIR}/fig_class_balance.png", bbox_inches="tight")
plt.show()
No description has been provided for this image
In [6]:
# 3b. Correlation heatmap on JM1 (largest dataset) — shows redundancy RFE will trim
Xj, yj = data["JM1"]
corr = Xj.corr()
fig, ax = plt.subplots(figsize=(11, 9))
sns.heatmap(corr, cmap="coolwarm", center=0, square=True,
            cbar_kws={"shrink": .7}, linewidths=.3, ax=ax)
ax.set_title("Feature correlation heatmap (JM1)")
plt.tight_layout()
plt.savefig(f"{OUTDIR}/fig_correlation_jm1.png", bbox_inches="tight")
plt.show()
No description has been provided for this image
In [7]:
# 3c. Descriptive statistics table (JM1) for the appendix
desc = Xj.describe().T[["mean", "std", "min", "50%", "max"]].round(2)
desc.to_csv(f"{OUTDIR}/table_jm1_describe.csv")
desc.head(22)
Out[7]:
mean std min 50% max
loc 39.68 71.90 1.0 22.00 3442.00
v(g) 6.02 12.19 1.0 3.00 470.00
ev(g) 3.21 6.47 1.0 1.00 165.00
iv(g) 3.82 8.48 1.0 2.00 402.00
n 108.86 234.81 0.0 47.00 8441.00
v 633.89 1804.95 0.0 206.58 80843.08
l 0.14 0.17 0.0 0.08 1.30
d 13.86 18.45 0.0 8.90 418.20
i 28.76 32.50 0.0 21.52 569.78
e 33656.99 401839.48 0.0 1877.33 31079782.27
b 0.21 0.60 0.0 0.07 26.95
t 1869.83 22324.42 0.0 104.30 1726654.57
lOCode 25.10 55.69 0.0 12.00 2824.00
lOComment 2.60 8.69 0.0 0.00 344.00
lOBlank 4.43 9.70 0.0 2.00 447.00
locCodeAndComment 0.35 1.81 0.0 0.00 108.00
uniq_Op 11.02 9.54 0.0 10.00 411.00
uniq_Opnd 16.15 24.92 0.0 11.00 1026.00
total_Op 64.78 142.73 0.0 28.00 5420.00
total_Opnd 44.20 94.31 0.0 19.00 3021.00
branchCount 10.67 21.29 1.0 5.00 826.00

4. Building the leak-free pipeline¶

In [8]:
def make_base_learners():
    rf = RandomForestClassifier(n_estimators=200, random_state=RANDOM_STATE, n_jobs=-1)
    xgb = XGBClassifier(n_estimators=200, learning_rate=0.1, max_depth=4,
                        subsample=0.9, colsample_bytree=0.9, eval_metric="logloss",
                        random_state=RANDOM_STATE, n_jobs=-1)
    svm = SVC(kernel="rbf", probability=True, random_state=RANDOM_STATE)
    return [("rf", rf), ("xgb", xgb), ("svm", svm)]

def make_stacking():
    return StackingClassifier(
        estimators=make_base_learners(),
        final_estimator=LogisticRegression(max_iter=1000, random_state=RANDOM_STATE),
        stack_method="predict_proba", cv=5, n_jobs=-1,
    )

def make_pipeline(model):
    """Wrap a classifier in the full leak-free preprocessing pipeline."""
    rfe = RFE(estimator=RandomForestClassifier(n_estimators=100, random_state=RANDOM_STATE),
              n_features_to_select=N_FEATURES)
    return ImbPipeline([
        ("impute", SimpleImputer(strategy="median")),
        ("scale",  MinMaxScaler()),
        ("smote",  SMOTE(random_state=RANDOM_STATE)),
        ("rfe",    rfe),
        ("model",  model),
    ])

# the proposed artefact plus the single-classifier baselines, all through the SAME pipeline
MODELS = {
    "Proposed (Stacking)": make_stacking,
    "Random Forest":       lambda: RandomForestClassifier(n_estimators=200, random_state=RANDOM_STATE, n_jobs=-1),
    "XGBoost":             lambda: XGBClassifier(n_estimators=200, learning_rate=0.1, max_depth=4,
                                                 eval_metric="logloss", random_state=RANDOM_STATE, n_jobs=-1),
    "SVM":                 lambda: SVC(kernel="rbf", probability=True, random_state=RANDOM_STATE),
    "Logistic Regression": lambda: LogisticRegression(max_iter=1000, random_state=RANDOM_STATE),
}
print("Pipeline factory ready. Models:", list(MODELS.keys()))

Pipeline factory ready. Models: ['Proposed (Stacking)', 'Random Forest', 'XGBoost', 'SVM', 'Logistic Regression']

5. Evaluation protocol (10-fold stratified CV)¶

In [9]:
SCORERS = {
    "Accuracy":  "accuracy",
    "Precision": "precision",
    "Recall":    "recall",
    "F1":        "f1",
    "AUC-ROC":   "roc_auc",
    "MCC":       "matthews_corrcoef",
}

def evaluate_model(model_factory, X, y):
    """Return a dict metric -> array of per-fold scores."""
    pipe = make_pipeline(model_factory())
    cv = StratifiedKFold(n_splits=N_FOLDS, shuffle=True, random_state=RANDOM_STATE)
    res = cross_validate(pipe, X, y, cv=cv, scoring=SCORERS, n_jobs=-1, error_score="raise")
    return {m: res[f"test_{m}"] for m in SCORERS}

print("Evaluation helpers ready.")

Evaluation helpers ready.

6. Run all models on all datasets¶

In [10]:
per_fold = {}   # per_fold[dataset][model][metric] = array of fold scores
mean_rows = []  # tidy table of mean±std

t0 = time.time()
for ds in DATASETS:
    X, y = data[ds]
    per_fold[ds] = {}
    for mname, factory in MODELS.items():
        scores = evaluate_model(factory, X, y)
        per_fold[ds][mname] = scores
        row = {"Dataset": ds, "Model": mname}
        for m in SCORERS:
            row[m] = round(float(np.mean(scores[m])), 4)
            row[m + "_std"] = round(float(np.std(scores[m])), 4)
        mean_rows.append(row)
        print(f"{ds:5} | {mname:20} | F1={row['F1']:.3f}  AUC={row['AUC-ROC']:.3f}  MCC={row['MCC']:.3f}")
print(f"\nDone in {time.time()-t0:.0f}s")

results_df = pd.DataFrame(mean_rows)
results_df.to_csv(f"{OUTDIR}/table_all_results.csv", index=False)
results_df

CM1   | Proposed (Stacking)  | F1=0.169  AUC=0.724  MCC=0.094
CM1   | Random Forest        | F1=0.172  AUC=0.743  MCC=0.090
CM1   | XGBoost              | F1=0.228  AUC=0.712  MCC=0.148
CM1   | SVM                  | F1=0.350  AUC=0.756  MCC=0.284
CM1   | Logistic Regression  | F1=0.358  AUC=0.790  MCC=0.308
KC1   | Proposed (Stacking)  | F1=0.445  AUC=0.804  MCC=0.367
KC1   | Random Forest        | F1=0.466  AUC=0.819  MCC=0.367
KC1   | XGBoost              | F1=0.437  AUC=0.791  MCC=0.340
KC1   | SVM                  | F1=0.437  AUC=0.788  MCC=0.332
KC1   | Logistic Regression  | F1=0.434  AUC=0.793  MCC=0.320
PC1   | Proposed (Stacking)  | F1=0.383  AUC=0.880  MCC=0.369
PC1   | Random Forest        | F1=0.429  AUC=0.869  MCC=0.397
PC1   | XGBoost              | F1=0.417  AUC=0.868  MCC=0.384
PC1   | SVM                  | F1=0.303  AUC=0.838  MCC=0.283
PC1   | Logistic Regression  | F1=0.261  AUC=0.814  MCC=0.209
JM1   | Proposed (Stacking)  | F1=0.408  AUC=0.806  MCC=0.334
JM1   | Random Forest        | F1=0.444  AUC=0.811  MCC=0.358
JM1   | XGBoost              | F1=0.359  AUC=0.743  MCC=0.271
JM1   | SVM                  | F1=0.403  AUC=0.726  MCC=0.268
JM1   | Logistic Regression  | F1=0.386  AUC=0.707  MCC=0.248

Done in 946s
Out[10]:
Dataset Model Accuracy Accuracy_std Precision Precision_std Recall Recall_std F1 F1_std AUC-ROC AUC-ROC_std MCC MCC_std
0 CM1 Proposed (Stacking) 0.8555 0.0354 0.1786 0.1969 0.1650 0.1761 0.1694 0.1819 0.7243 0.1069 0.0944 0.1939
1 CM1 Random Forest 0.8414 0.0351 0.1521 0.1726 0.2050 0.2534 0.1724 0.1997 0.7432 0.0980 0.0905 0.2168
2 CM1 XGBoost 0.8515 0.0201 0.2114 0.1335 0.2500 0.1857 0.2278 0.1549 0.7120 0.1383 0.1485 0.1624
3 CM1 SVM 0.7692 0.0444 0.2436 0.0771 0.6300 0.2002 0.3498 0.1074 0.7560 0.1077 0.2840 0.1451
4 CM1 Logistic Regression 0.7492 0.0476 0.2427 0.0666 0.7150 0.2025 0.3583 0.0900 0.7902 0.1014 0.3082 0.1252
5 KC1 Proposed (Stacking) 0.8483 0.0113 0.5189 0.0518 0.3982 0.0711 0.4448 0.0476 0.8036 0.0360 0.3668 0.0422
6 KC1 Random Forest 0.8260 0.0259 0.4553 0.0748 0.4903 0.0772 0.4657 0.0491 0.8189 0.0362 0.3673 0.0611
7 KC1 XGBoost 0.8269 0.0207 0.4507 0.0614 0.4384 0.0802 0.4371 0.0446 0.7909 0.0324 0.3403 0.0486
8 KC1 SVM 0.6932 0.0324 0.3063 0.0327 0.7665 0.0522 0.4371 0.0381 0.7876 0.0345 0.3321 0.0525
9 KC1 Logistic Regression 0.7198 0.0286 0.3167 0.0327 0.6927 0.0546 0.4341 0.0383 0.7930 0.0333 0.3196 0.0522
10 PC1 Proposed (Stacking) 0.9252 0.0189 0.5195 0.2477 0.3500 0.1631 0.3831 0.1494 0.8804 0.0472 0.3689 0.1513
11 PC1 Random Forest 0.9189 0.0184 0.4424 0.1572 0.4554 0.1765 0.4293 0.1330 0.8694 0.0657 0.3972 0.1411
12 PC1 XGBoost 0.9180 0.0218 0.4426 0.1525 0.4286 0.1559 0.4171 0.1242 0.8684 0.0580 0.3839 0.1319
13 PC1 SVM 0.7719 0.0407 0.1953 0.0416 0.7036 0.1127 0.3031 0.0525 0.8382 0.0451 0.2829 0.0663
14 PC1 Logistic Regression 0.7989 0.0382 0.1775 0.0701 0.5054 0.1953 0.2608 0.0997 0.8140 0.0693 0.2089 0.1308
15 JM1 Proposed (Stacking) 0.8450 0.0099 0.5204 0.0504 0.3362 0.0371 0.4082 0.0417 0.8060 0.0123 0.3342 0.0472
16 JM1 Random Forest 0.8432 0.0078 0.5101 0.0331 0.3932 0.0300 0.4438 0.0295 0.8112 0.0129 0.3585 0.0335
17 JM1 XGBoost 0.8282 0.0088 0.4433 0.0392 0.3024 0.0314 0.3590 0.0325 0.7434 0.0220 0.2707 0.0363
18 JM1 SVM 0.7253 0.0096 0.3084 0.0148 0.5839 0.0453 0.4033 0.0220 0.7262 0.0242 0.2679 0.0296
19 JM1 Logistic Regression 0.7495 0.0062 0.3168 0.0132 0.4959 0.0352 0.3864 0.0195 0.7067 0.0219 0.2482 0.0240

6.1 Headline results table (core four datasets, mean over folds)¶

In [11]:
core_df = results_df[results_df["Dataset"].isin(CORE)]
headline = core_df.pivot_table(index=["Dataset", "Model"],
                               values=list(SCORERS.keys()))
headline = headline.reindex(columns=list(SCORERS.keys()))
headline.to_csv(f"{OUTDIR}/table_headline_core.csv")
headline.round(3)
Out[11]:
Accuracy Precision Recall F1 AUC-ROC MCC
Dataset Model
CM1 Logistic Regression 0.749 0.243 0.715 0.358 0.790 0.308
Proposed (Stacking) 0.856 0.179 0.165 0.169 0.724 0.094
Random Forest 0.841 0.152 0.205 0.172 0.743 0.090
SVM 0.769 0.244 0.630 0.350 0.756 0.284
XGBoost 0.852 0.211 0.250 0.228 0.712 0.148
JM1 Logistic Regression 0.750 0.317 0.496 0.386 0.707 0.248
Proposed (Stacking) 0.845 0.520 0.336 0.408 0.806 0.334
Random Forest 0.843 0.510 0.393 0.444 0.811 0.358
SVM 0.725 0.308 0.584 0.403 0.726 0.268
XGBoost 0.828 0.443 0.302 0.359 0.743 0.271
KC1 Logistic Regression 0.720 0.317 0.693 0.434 0.793 0.320
Proposed (Stacking) 0.848 0.519 0.398 0.445 0.804 0.367
Random Forest 0.826 0.455 0.490 0.466 0.819 0.367
SVM 0.693 0.306 0.766 0.437 0.788 0.332
XGBoost 0.827 0.451 0.438 0.437 0.791 0.340
PC1 Logistic Regression 0.799 0.178 0.505 0.261 0.814 0.209
Proposed (Stacking) 0.925 0.520 0.350 0.383 0.880 0.369
Random Forest 0.919 0.442 0.455 0.429 0.869 0.397
SVM 0.772 0.195 0.704 0.303 0.838 0.283
XGBoost 0.918 0.443 0.429 0.417 0.868 0.384

6.2 Proposed vs best baseline grouped bar chart per dataset¶

In [12]:
metrics_to_plot = ["F1", "AUC-ROC", "MCC"]
for ds in CORE:
    sub = results_df[results_df["Dataset"] == ds].set_index("Model")
    fig, ax = plt.subplots(figsize=(8, 4.5))
    sub[metrics_to_plot].plot(kind="bar", ax=ax, edgecolor="black")
    ax.set_title(f"{ds}: proposed ensemble vs baselines")
    ax.set_ylabel("Score"); ax.set_ylim(0, 1)
    ax.legend(title="Metric", loc="lower right")
    plt.xticks(rotation=20, ha="right")
    plt.tight_layout()
    plt.savefig(f"{OUTDIR}/fig_bars_{ds}.png", bbox_inches="tight")
    plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

7. Statistical significance (Wilcoxon signed-rank)¶

In [13]:
sig_rows = []
PROP = "Proposed (Stacking)"
for ds in DATASETS:
    for mname in MODELS:
        if mname == PROP:
            continue
        for metric in ["F1", "AUC-ROC", "MCC"]:
            a = per_fold[ds][PROP][metric]
            b = per_fold[ds][mname][metric]
            try:
                stat, p = wilcoxon(a, b)
            except ValueError:
                p = np.nan   # happens if all differences are zero
            sig_rows.append({
                "Dataset": ds, "Baseline": mname, "Metric": metric,
                "Proposed_mean": round(float(np.mean(a)), 4),
                "Baseline_mean": round(float(np.mean(b)), 4),
                "p_value": round(float(p), 4) if p == p else np.nan,
                "Significant(<0.05)": ("yes" if (p == p and p < 0.05) else "no"),
            })
sig_df = pd.DataFrame(sig_rows)
sig_df.to_csv(f"{OUTDIR}/table_wilcoxon.csv", index=False)
sig_df[sig_df["Dataset"].isin(CORE)]
Out[13]:
Dataset Baseline Metric Proposed_mean Baseline_mean p_value Significant(<0.05)
0 CM1 Random Forest F1 0.1694 0.1724 1.0000 no
1 CM1 Random Forest AUC-ROC 0.7243 0.7432 0.0840 no
2 CM1 Random Forest MCC 0.0944 0.0905 0.6406 no
3 CM1 XGBoost F1 0.1694 0.2278 0.3906 no
4 CM1 XGBoost AUC-ROC 0.7243 0.7120 0.9219 no
5 CM1 XGBoost MCC 0.0944 0.1485 0.4258 no
6 CM1 SVM F1 0.1694 0.3498 0.0410 yes
7 CM1 SVM AUC-ROC 0.7243 0.7560 0.6250 no
8 CM1 SVM MCC 0.0944 0.2840 0.0645 no
9 CM1 Logistic Regression F1 0.1694 0.3583 0.0137 yes
10 CM1 Logistic Regression AUC-ROC 0.7243 0.7902 0.0273 yes
11 CM1 Logistic Regression MCC 0.0944 0.3082 0.0098 yes
12 KC1 Random Forest F1 0.4448 0.4657 0.2754 no
13 KC1 Random Forest AUC-ROC 0.8036 0.8189 0.0020 yes
14 KC1 Random Forest MCC 0.3668 0.3673 1.0000 no
15 KC1 XGBoost F1 0.4448 0.4371 0.7695 no
16 KC1 XGBoost AUC-ROC 0.8036 0.7909 0.0371 yes
17 KC1 XGBoost MCC 0.3668 0.3403 0.1309 no
18 KC1 SVM F1 0.4448 0.4371 0.3223 no
19 KC1 SVM AUC-ROC 0.8036 0.7876 0.1934 no
20 KC1 SVM MCC 0.3668 0.3321 0.1055 no
21 KC1 Logistic Regression F1 0.4448 0.4341 0.4922 no
22 KC1 Logistic Regression AUC-ROC 0.8036 0.7930 0.0840 no
23 KC1 Logistic Regression MCC 0.3668 0.3196 0.0645 no
24 PC1 Random Forest F1 0.3831 0.4293 0.0781 no
25 PC1 Random Forest AUC-ROC 0.8804 0.8694 0.2324 no
26 PC1 Random Forest MCC 0.3689 0.3972 0.2031 no
27 PC1 XGBoost F1 0.3831 0.4171 0.3223 no
28 PC1 XGBoost AUC-ROC 0.8804 0.8684 0.2422 no
29 PC1 XGBoost MCC 0.3689 0.3839 0.6250 no
30 PC1 SVM F1 0.3831 0.3031 0.1602 no
31 PC1 SVM AUC-ROC 0.8804 0.8382 0.0840 no
32 PC1 SVM MCC 0.3689 0.2829 0.1055 no
33 PC1 Logistic Regression F1 0.3831 0.2608 0.0273 yes
34 PC1 Logistic Regression AUC-ROC 0.8804 0.8140 0.0273 yes
35 PC1 Logistic Regression MCC 0.3689 0.2089 0.0195 yes
36 JM1 Random Forest F1 0.4082 0.4438 0.0020 yes
37 JM1 Random Forest AUC-ROC 0.8060 0.8112 0.0273 yes
38 JM1 Random Forest MCC 0.3342 0.3585 0.0098 yes
39 JM1 XGBoost F1 0.4082 0.3590 0.0098 yes
40 JM1 XGBoost AUC-ROC 0.8060 0.7434 0.0020 yes
41 JM1 XGBoost MCC 0.3342 0.2707 0.0059 yes
42 JM1 SVM F1 0.4082 0.4033 0.6250 no
43 JM1 SVM AUC-ROC 0.8060 0.7262 0.0020 yes
44 JM1 SVM MCC 0.3342 0.2679 0.0059 yes
45 JM1 Logistic Regression F1 0.4082 0.3864 0.0645 no
46 JM1 Logistic Regression AUC-ROC 0.8060 0.7067 0.0020 yes
47 JM1 Logistic Regression MCC 0.3342 0.2482 0.0039 yes

8. Explainability (SHAP global + LIME local)¶

The pipeline is fitted once on a chosen dataset, then SHAP gives the global ranking of the most influential defect-predictive metrics and LIME explains individual module predictions. To keep SHAP tractable, it is applied to the fitted Random Forest base learner on the RFE-selected, resampled feature space (the tree explainer is fast and exact for RF). Change XAI_DATASET to run the analysis on a different dataset.

In [14]:
import shap
from lime.lime_tabular import LimeTabularExplainer
from sklearn.model_selection import train_test_split

XAI_DATASET = "KC1"
Xx, yx = data[XAI_DATASET]

# Fit the full leak-free pipeline on a train split, explain on it
X_tr, X_te, y_tr, y_te = train_test_split(Xx, yx, test_size=0.3,
                                          stratify=yx, random_state=RANDOM_STATE)
xai_pipe = make_pipeline(RandomForestClassifier(n_estimators=200, random_state=RANDOM_STATE, n_jobs=-1))
xai_pipe.fit(X_tr, y_tr)

# recover the transformed feature space and selected feature names
Xt = X_tr.copy()
Xt = pd.DataFrame(xai_pipe.named_steps["impute"].transform(Xt), columns=Xx.columns)
Xt = pd.DataFrame(xai_pipe.named_steps["scale"].transform(Xt), columns=Xx.columns)
# SMOTE only changes rows, not columns; RFE selects columns:
rfe = xai_pipe.named_steps["rfe"]
selected = Xx.columns[rfe.support_].tolist()
print(f"RFE selected {len(selected)} features on {XAI_DATASET}:")
print(selected)

RFE selected 15 features on KC1:
['loc', 'iv(g)', 'n', 'v', 'l', 'd', 'i', 'e', 'b', 't', 'lOCode', 'lOBlank', 'uniq_Opnd', 'total_Op', 'total_Opnd']
In [15]:
# SHAP global on the fitted RF base learner over the selected feature space
rf_final = xai_pipe.named_steps["model"]
X_tr_sel = Xt[selected]

explainer = shap.TreeExplainer(rf_final)
Xsample = X_tr_sel.sample(min(800, len(X_tr_sel)), random_state=RANDOM_STATE)
sv = explainer.shap_values(Xsample)

# take the positive (defective) class; handle both shap output formats
if isinstance(sv, list):
    sv_pos = sv[1]
elif sv.ndim == 3:
    sv_pos = sv[:, :, 1]
else:
    sv_pos = sv

# bar plot (global importance)
plt.figure(figsize=(8, 6))
shap.summary_plot(sv_pos, Xsample, plot_type="bar", show=False, max_display=15)
plt.title(f"SHAP global feature importance ({XAI_DATASET})")
plt.tight_layout()
plt.savefig(f"{OUTDIR}/fig_shap_bar_{XAI_DATASET}.png", dpi=200, bbox_inches="tight")
plt.show()

# beeswarm (direction of effect)
plt.figure(figsize=(8, 6))
shap.summary_plot(sv_pos, Xsample, show=False, max_display=15)
plt.title(f"SHAP summary beeswarm ({XAI_DATASET})")
plt.tight_layout()
plt.savefig(f"{OUTDIR}/fig_shap_beeswarm_{XAI_DATASET}.png", dpi=200, bbox_inches="tight")
plt.show()
No description has been provided for this image
No description has been provided for this image
In [ ]:
In [16]:
# LIME local explanation for one defective and one clean module
lime_exp = LimeTabularExplainer(
    training_data=X_tr_sel.values,
    feature_names=selected,
    class_names=["clean", "defective"],
    mode="classification", random_state=RANDOM_STATE,
)
# transform the test set the same way, restrict to selected features
X_te_t = pd.DataFrame(xai_pipe.named_steps["impute"].transform(X_te), columns=Xx.columns)
X_te_t = pd.DataFrame(xai_pipe.named_steps["scale"].transform(X_te_t), columns=Xx.columns)[selected]

def explain_instance(idx, tag):
    exp = lime_exp.explain_instance(X_te_t.iloc[idx].values, rf_final.predict_proba,
                                    num_features=10)
    fig = exp.as_pyplot_figure()
    fig.set_size_inches(7, 4)
    plt.title(f"LIME local explanation ({XAI_DATASET}, {tag} module)")
    plt.tight_layout()
    plt.savefig(f"{OUTDIR}/fig_lime_{XAI_DATASET}_{tag}.png", bbox_inches="tight")
    plt.show()

# pick a predicted-defective and a predicted-clean instance
preds = rf_final.predict(X_te_t.values)
def_idx = int(np.where(preds == 1)[0][0]) if (preds == 1).any() else 0
cln_idx = int(np.where(preds == 0)[0][0]) if (preds == 0).any() else 1
explain_instance(def_idx, "defective")
explain_instance(cln_idx, "clean")
No description has been provided for this image
No description has been provided for this image

9. Leak-free vs leaky preprocessing the study's own contribution¶

In [17]:
def evaluate_leaky(model_factory, X, y):
    """WRONG-ON-PURPOSE: preprocess the whole dataset, THEN cross-validate the model alone."""
    Xi = SimpleImputer(strategy="median").fit_transform(X)
    Xs = MinMaxScaler().fit_transform(Xi)
    # SMOTE on everything (leak), then RFE on everything (leak)
    Xr, yr = SMOTE(random_state=RANDOM_STATE).fit_resample(Xs, y)
    rfe = RFE(RandomForestClassifier(n_estimators=100, random_state=RANDOM_STATE),
              n_features_to_select=N_FEATURES).fit(Xr, yr)
    Xsel = rfe.transform(Xr)
    cv = StratifiedKFold(n_splits=N_FOLDS, shuffle=True, random_state=RANDOM_STATE)
    res = cross_validate(model_factory(), Xsel, yr, cv=cv, scoring=SCORERS, n_jobs=-1)
    return {m: res[f"test_{m}"] for m in SCORERS}

leak_rows = []
for ds in CORE:
    X, y = data[ds]
    leaky = evaluate_leaky(make_stacking, X, y)
    correct = per_fold[ds][PROP]
    for m in SCORERS:
        leak_rows.append({
            "Dataset": ds, "Metric": m,
            "Leak-free (correct)": round(float(np.mean(correct[m])), 4),
            "Leaky (inflated)":    round(float(np.mean(leaky[m])), 4),
            "Inflation":           round(float(np.mean(leaky[m]) - np.mean(correct[m])), 4),
        })
leak_df = pd.DataFrame(leak_rows)
leak_df.to_csv(f"{OUTDIR}/table_leakage_ablation.csv", index=False)
leak_df
Out[17]:
Dataset Metric Leak-free (correct) Leaky (inflated) Inflation
0 CM1 Accuracy 0.8555 0.9354 0.0800
1 CM1 Precision 0.1786 0.9171 0.7385
2 CM1 Recall 0.1650 0.9578 0.7928
3 CM1 F1 0.1694 0.9366 0.7672
4 CM1 AUC-ROC 0.7243 0.9880 0.2636
5 CM1 MCC 0.0944 0.8725 0.7781
6 KC1 Accuracy 0.8483 0.9136 0.0654
7 KC1 Precision 0.5189 0.9206 0.4018
8 KC1 Recall 0.3982 0.9058 0.5076
9 KC1 F1 0.4448 0.9129 0.4681
10 KC1 AUC-ROC 0.8036 0.9674 0.1638
11 KC1 MCC 0.3668 0.8277 0.4609
12 PC1 Accuracy 0.9252 0.9700 0.0448
13 PC1 Precision 0.5195 0.9696 0.4501
14 PC1 Recall 0.3500 0.9710 0.6210
15 PC1 F1 0.3831 0.9700 0.5869
16 PC1 AUC-ROC 0.8804 0.9954 0.1150
17 PC1 MCC 0.3689 0.9404 0.5715
18 JM1 Accuracy 0.8450 0.9199 0.0749
19 JM1 Precision 0.5204 0.9332 0.4128
20 JM1 Recall 0.3362 0.9045 0.5683
21 JM1 F1 0.4082 0.9186 0.5104
22 JM1 AUC-ROC 0.8060 0.9710 0.1650
23 JM1 MCC 0.3342 0.8402 0.5060
In [18]:
# Visualise the inflation for AUC-ROC and MCC
for metric in ["AUC-ROC", "MCC", "F1"]:
    sub = leak_df[leak_df["Metric"] == metric].set_index("Dataset")
    fig, ax = plt.subplots(figsize=(8, 4.5))
    sub[["Leak-free (correct)", "Leaky (inflated)"]].plot(kind="bar", ax=ax,
        color=["#2e7d32", "#c62828"], edgecolor="black")
    ax.set_title(f"Effect of data leakage on {metric} (proposed ensemble)")
    ax.set_ylabel(metric); ax.set_ylim(0, 1)
    plt.xticks(rotation=0)
    plt.tight_layout()
    plt.savefig(f"{OUTDIR}/fig_leakage_{metric.replace('-', '')}.png", bbox_inches="tight")
    plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

10. Export summary¶

In [19]:
produced = sorted(os.listdir(OUTDIR))
print(f"{len(produced)} files written to ./{OUTDIR}/:\n")
for f in produced:
    print("  ", f)

# also dump a machine-readable roll-up of all per-fold scores for the appendix / reproducibility
roll = {ds: {m: {k: list(map(float, v)) for k, v in per_fold[ds][m].items()}
             for m in per_fold[ds]} for ds in per_fold}
with open(f"{OUTDIR}/per_fold_scores.json", "w") as fh:
    json.dump(roll, fh, indent=2)
print("\nAll per-fold scores saved to per_fold_scores.json")

20 files written to ./outputs/:

   fig_bars_CM1.png
   fig_bars_JM1.png
   fig_bars_KC1.png
   fig_bars_PC1.png
   fig_class_balance.png
   fig_correlation_jm1.png
   fig_leakage_AUCROC.png
   fig_leakage_F1.png
   fig_leakage_MCC.png
   fig_lime_KC1_clean.png
   fig_lime_KC1_defective.png
   fig_shap_bar_KC1.png
   fig_shap_beeswarm_KC1.png
   per_fold_scores.json
   table_all_results.csv
   table_dataset_summary.csv
   table_headline_core.csv
   table_jm1_describe.csv
   table_leakage_ablation.csv
   table_wilcoxon.csv

All per-fold scores saved to per_fold_scores.json