Non-linear Real Estate Price Prediction (Ames Housing):

# standard python libs
import os
import random
import time
import warnings
from typing import Dict, List, Tuple # this one helps me for typing.

warnings.filterwarnings("ignore")

# importing machine learning stuff
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

from scipy import stats

# SKLearn stuff
from sklearnex import patch_sklearn
patch_sklearn()

from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.inspection import permutation_importance
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import (
    KFold,
    RandomizedSearchCV,
    GridSearchCV,
    cross_validate,
    learning_curve,
    train_test_split,
)
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import FunctionTransformer, OneHotEncoder, StandardScaler
from sklearn.svm import SVR

sns.set_theme(style="whitegrid")

RANDOM_STATE = 67
np.random.seed(RANDOM_STATE)
random.seed(RANDOM_STATE)

print(f"[INFO] Random state fixed at {RANDOM_STATE}")

# 
SALEPRICE_PLOT_COLOR = "steelblue"
LOGSALEPRICE_PLOT_COLOR = "darkorange"

DATASET_PATH = "/kaggle/input/competitions/house-prices-advanced-regression-techniques/train.csv"
df = pd.read_csv(DATASET_PATH)
df.shape

display(df.head())
print(f"Shape: {df.shape}")
print(f"Columns: {len(df.columns)}\n")
print("Data types:")
display(df.dtypes.value_counts())

missing = df.isna().mean().sort_values(ascending=False)
missing_top = missing[missing > 0].head(20)

fig, axes = plt.subplots(1, 2, figsize=(16, 5))
sns.histplot(df['SalePrice'], kde=True, ax=axes[0], color='teal')
axes[0].set_title('SalePrice Distribution (Raw)')

if len(missing_top) > 0:
    sns.barplot(x=missing_top.values, y=missing_top.index, ax=axes[1], palette='viridis', hue=missing_top.index, legend=False)
    axes[1].set_title('Top Missingness Ratios')
    axes[1].set_xlabel('Missing Ratio')
else:
    axes[1].text(0.5, 0.5, 'No Missing Values', ha='center', va='center')
    axes[1].set_axis_off()

plt.tight_layout()
plt.show()

if "Id" in df.columns:
    df = df.drop(columns = ["Id"])
df["LogSalePrice"] = np.log1p(df["SalePrice"])

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sns.histplot(df["SalePrice"], kde=True, ax=axes[0], color=SALEPRICE_PLOT_COLOR)
axes[0].set_title("Raw SalePrice")
sns.histplot(df["LogSalePrice"], kde=True, ax=axes[1], color=LOGSALEPRICE_PLOT_COLOR)
axes[1].set_title("Log(1 + SalePrice)")
plt.tight_layout()
plt.show()

def engineerFeatures(inputDf: pd.DataFrame) -> pd.DataFrame:
    data = inputDf.copy()

    data["TotalSF"] = 0
    for col in ["TotalBsmtSF", "1stFlrSF", "2ndFlrSF"]:
        if col not in data.columns:
            data[col] = 0
        data["TotalSF"] += data[col].fillna(0)

    if {'YrSold', 'YearBuilt'}.issubset(data.columns):
        data['HouseAgeAtSale'] = data['YrSold'].fillna(data['YrSold'].median()) - data['YearBuilt'].fillna(data['YearBuilt'].median())

    if {'YrSold', 'YearRemodAdd'}.issubset(data.columns):
        data['RemodAgeAtSale'] = data['YrSold'].fillna(data['YrSold'].median()) - data['YearRemodAdd'].fillna(data['YearRemodAdd'].median())

    full_bath_cols = [c for c in ['FullBath', 'BsmtFullBath'] if c in data.columns]
    half_bath_cols = [c for c in ['HalfBath', 'BsmtHalfBath'] if c in data.columns]
    data['TotalBaths'] = data[full_bath_cols].fillna(0).sum(axis=1) + 0.5 * data[half_bath_cols].fillna(0).sum(axis=1)

    if 'GarageArea' in data.columns:
        data['HasGarage'] = (data['GarageArea'].fillna(0) > 0).astype(int)
    if 'TotalBsmtSF' in data.columns:
        data['HasBsmt'] = (data['TotalBsmtSF'].fillna(0) > 0).astype(int)

    return data

x = df.drop(columns=['SalePrice', 'LogSalePrice'])
y = df['LogSalePrice']

xTrain, xTest, yTrain, yTest = train_test_split(x, y, test_size=0.2, random_state=RANDOM_STATE)

print(f"Train shape: {xTrain.shape}, Test shape: {xTest.shape}")

VIF_THRESHOLD = 10.0
MAX_DROPS = 15

def computeVIF(numericDF: pd.DataFrame) -> pd.DataFrame:
    work = numericDF.select_dtypes(include=[np.number]).copy()
    work = work.replace([np.inf, -np.inf], np.nan)

    # Median imputation per column; if a column is fully missing, fallback to 0.
    for col in work.columns:
        colMedian = work[col].median()
        if pd.isna(colMedian):
            colMedian = 0.0
        work[col] = work[col].fillna(colMedian)
        
    vals = []
    cols = list(work.columns)

    for col in cols:
        yCol = work[col]
        xCols = work.drop(columns=[col])
        if xCols.shape[1] == 0:
            vals.append((col, 1.0))
            continue

        model = LinearRegression()
        model.fit(xCols, yCol)
        r2 = model.score(xCols, yCol)
        vif = 1.0 / max(1e-8, (1.0 - r2))
        vals.append((col, float(vif)))

    return  pd.DataFrame(vals, columns=["feature", "vif"]).sort_values('vif', ascending=False) if vals else pd.DataFrame(vals, columns=["feature", "vif"])

XTrainFe = engineerFeatures(xTrain)
numericVIFCols = XTrainFe.select_dtypes(include=[np.number]).columns.tolist()
vifDF = XTrainFe[numericVIFCols].copy()

droppedMulticollinear: List[str] = []

for _ in range(MAX_DROPS):
    currentVIF = computeVIF(vifDF)
    if currentVIF.empty:
        break
        
    worst = currentVIF.iloc[0]
    if worst['vif'] <= VIF_THRESHOLD:
        break
    dropFeature = worst["feature"]
    droppedMulticollinear.append(dropFeature)
    vifDF = vifDF.drop(columns=[dropFeature], errors="ignore")

finalVIF = computeVIF(vifDF)
finalVIFHead = finalVIF.head(15)
print("Dropped due to multicollinearity:")
print(droppedMulticollinear if droppedMulticollinear else "None")
display(finalVIFHead)

plt.figure(figsize=(10, 6))
sns.barplot(data=finalVIFHead, x="vif", y="feature", palette='mako', hue="feature", legend=False)
plt.xlabel("VIF (Variance Inflation Factor)")
plt.ylabel("Feature")
plt.title("Top Remaining VIF Values (Plot Mitigation)")
plt.tight_layout()
plt.show()

def dropSelectedColumns(data: pd.DataFrame, targets: List[str]) -> pd.DataFrame:
    return data.drop(targets, errors="ignore")

xTrainPrepared = dropSelectedColumns(engineerFeatures(xTrain), droppedMulticollinear)

numFeatures = xTrainPrepared.select_dtypes(include=[np.number]).columns.tolist()
catFeatures = xTrainPrepared.select_dtypes(exclude=[np.number]).columns.tolist()

numericTransformLinear = Pipeline(steps=[
    ("imputer", SimpleImputer(strategy="median"))
])

numericTransformScaled = Pipeline(steps=[
    ("imputer", SimpleImputer(strategy="median")),
    ("scaler", StandardScaler())
])

categoricalTransform = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='most_frequent')),
    ('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))
])

preprocessorLinear = ColumnTransformer(transformers=[
    ("num", numericTransformLinear, numFeatures),
    ("cat", categoricalTransform, catFeatures)
])

preprocessorScaled = ColumnTransformer(transformers=[
    ('num', numericTransformScaled, numFeatures),
    ('cat', categoricalTransform, catFeatures),
])

featureEngineeringStep = FunctionTransformer(engineerFeatures, validate=False)
dropMulticolStep = FunctionTransformer(lambda d: dropSelectedColumns(d, droppedMulticollinear), validate=False)

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

print(f'Numeric features: {len(numFeatures)}')
print(f'Categorical features: {len(catFeatures)}')

def rmse(yTrue, yPred):
    return np.sqrt(mean_squared_error(yTrue, yPred))

pipeLR = Pipeline(steps=[
    ("feat", featureEngineeringStep),
    ("drop_multi", dropMulticolStep),
    ("prep", preprocessorLinear),
    ("model", LinearRegression())
])

pipeSVR = Pipeline(steps=[
    ("feat", featureEngineeringStep),
    ("drop_multi", dropMulticolStep),
    ("prep", preprocessorScaled),
    ("model", SVR())
])

pipeANN = Pipeline(steps=[
    ("feat", featureEngineeringStep),
    ("drop_multi", dropMulticolStep),
    ("prep", preprocessorLinear),
    ("model", MLPRegressor(
        random_state = RANDOM_STATE,
        max_iter = 2000,
        early_stopping = True,
        validation_fraction=0.1
    ))
])

gridLR = {
    'model__fit_intercept': [True, False],
    'model__positive': [False, True],
}

randSVR = {
    'model__kernel': ['rbf'],
    'model__C': np.logspace(-1, 2.5, 30),
    'model__epsilon': np.linspace(0.01, 0.5, 30),
    'model__gamma': ['scale', 'auto'],
}

randANN = {
    'model__hidden_layer_sizes': [(64,), (128,), (64, 32), (128, 64)],
    'model__activation': ['relu', 'tanh'],
    'model__alpha': np.logspace(-6, -2, 20),
    'model__learning_rate_init': np.logspace(-4, -2, 20),
}

searches = {
    "Linear Regression": GridSearchCV(
        estimator= pipeLR,
        param_grid= gridLR,
        scoring= "neg_root_mean_squared_error",
        cv= cv,
        n_jobs= 1,
        return_train_score= True
    ),
    "SVR": RandomizedSearchCV(
        estimator= pipeSVR,
        param_distributions= randSVR,
        n_iter= 25,
        random_state= RANDOM_STATE,
        scoring= "neg_root_mean_squared_error",
        cv= cv,
        n_jobs= 1,
        return_train_score= True
    ),
    "ANN": RandomizedSearchCV(
        estimator= pipeANN,
        param_distributions= randANN,
        n_iter= 25,
        random_state= RANDOM_STATE,
        scoring= "neg_root_mean_squared_error",
        cv= cv,
        n_jobs= 1,
        return_train_score= True
    )
}

bestModels = {}
searchSummaries = []

for modelName, search in searches.items():
    print(f"\nTraining {modelName}...")
    start= time.perf_counter()
    search.fit(xTrain, yTrain)
    elapsed= time.perf_counter() - start

    bestModels[modelName] = search.best_estimator_
    searchSummaries.append({
        "model": modelName,
        "bestCVRmse": -search.best_score_,
        "bestTimeSec": elapsed,
        "MeanFitTimeBest": search.cv_results_['mean_fit_time'][search.best_index_],
        "bestParams": search.best_params_
    })

summaryDF = pd.DataFrame(searchSummaries).sort_values("bestCVRmse")
display(summaryDF[["model", "bestCVRmse", "bestTimeSec", "MeanFitTimeBest"]])
summaryDF

evaluationRows = []
foldScores = {}

for model_name, model in bestModels.items():
    fitStart = time.perf_counter()
    model.fit(xTrain, yTrain)
    fitTime = time.perf_counter() - fitStart

    predStart = time.perf_counter()
    predTest = model.predict(xTest)
    predTime = time.perf_counter() - predStart

    cvResult = cross_validate(
        model,
        xTrain,
        yTrain,
        cv=cv,
        scoring="neg_root_mean_squared_error",
        return_train_score=True,
        n_jobs=1
    )

    trainRmse = -cvResult["train_score"]
    valRmse = -cvResult["test_score"]
    foldScores[model_name] = valRmse

    evaluationRows.append({
        'model': model_name,
        'testRmse': rmse(yTest, predTest),
        'testMae': mean_absolute_error(yTest, predTest),
        'testR2': r2_score(yTest, predTest),
        'cvTrainRmseMean': trainRmse.mean(),
        'cvValRmseMean': valRmse.mean(),
        'biasVarianceGap': valRmse.mean() - trainRmse.mean(),
        'fitTimeSec': fitTime,
        'predTimeSec': predTime,
    })

evalDF = pd.DataFrame(evaluationRows).sort_values("testRmse")
display(evalDF)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sns.barplot(data=evalDF, x="model", y="testRmse", ax=axes[0], palette="crest", hue="testRmse", legend=True)
axes[0].set_title("Test RMSE (lower is better)")
axes[0].set_xlabel("Model")
axes[0].set_ylabel("Test RMSE")
axes[0].tick_params(axis="x", rotation=15)

sns.barplot(data=evalDF, x="model", y="biasVarianceGap", ax=axes[1], palette="flare", hue="biasVarianceGap", legend=True)
axes[1].axhline(0, linestyle="--", color="black", linewidth=1)
axes[1].set_title("Bias-Variance Gap: CV Val RMSE - CV Train RMSE")
axes[1].set_xlabel("Model")
axes[1].set_ylabel("Bias-Variance Gap")
axes[1].tick_params(axis="x", rotation=15)

plt.tight_layout()
plt.show()

modelNames = list(foldScores.keys())
stateRows = []

for i in range(len(modelNames)):
    for j in range(i + 1, len(modelNames)):
        a = modelNames[i]
        b = modelNames[j]
        tStat, pVal = stats.ttest_rel(foldScores[a], foldScores[b])
        stateRows.append({
            "modelA": a,
            "modelB": b,
            "meanRmseA": np.mean(foldScores[a]),
            "meanRmseB": np.mean(foldScores[b]),
            "tStat": tStat,
            "pValue": pVal,
        })

statsDF = pd.DataFrame(stateRows).sort_values("pValue")
display(statsDF)

print('Interpretation tip: pValue < 0.05 suggests statistically significant fold-level performance difference.')

fig, axes = plt.subplots(1, 3, figsize=(18, 5), sharey=True)
trainSizes = np.linspace(0.2, 1.0, 5)

for ax, (model_name, model) in zip(axes, bestModels.items()):
    sizes, trainScores, valScores = learning_curve(
        model,
        xTrain,
        yTrain,
        cv=cv,
        train_sizes=trainSizes,
        scoring="neg_root_mean_squared_error",
        n_jobs = 1
    )

    trainRmse = -trainScores.mean(axis=1)
    valRmse = -valScores.mean(axis=1)

    ax.plot(sizes, trainRmse, marker="o", label="Train RMSE")
    ax.plot(sizes, valRmse, marker="s", label="CV RMSE")
    ax.set_title(f"{model_name} learning curve")
    ax.set_xlabel("Training Samples")
    ax.grid(alpha=0.3)

axes[0].set_ylabel("RMSE")
axes[0].legend()
plt.tight_layout()
plt.show()

def getFeatureNamesFromPreprocessor(preprocessor: ColumnTransformer) -> np.ndarray:
    return preprocessor.get_feature_names_out()

impactFrames = []

 # Linear coefs:
LRModel = bestModels["Linear Regression"]
LRPrep = LRModel.named_steps["prep"]
LREst = LRModel.named_steps["model"]
LRFeatureNames = getFeatureNamesFromPreprocessor(LRPrep)

coefDF = pd.DataFrame({
    "feature": LRFeatureNames,
    "impact": LREst.coef_
})

coefDF["absImpact"] = coefDF["impact"].abs()
coefTop = coefDF.sort_values("absImpact", ascending=False).head(20)
coefTop["model"] = "Linear Regression (coefficient)"
impactFrames.append(coefTop[["model", "feature", "impact", "absImpact"]])

# Permutation importance for non linear models
for modelName in ["SVR", "ANN"]:
    mdl = bestModels[modelName]
    pi = permutation_importance(
        mdl,
        xTest,
        yTest,
        n_repeats=10,
        random_state=RANDOM_STATE,
        scoring="neg_root_mean_squared_error",
        n_jobs=1
    )

    piDF = pd.DataFrame({
        "feature": xTest.columns,
        "impact": pi.importances_mean
    })
    piDF["absImpact"] = piDF["impact"].abs()
    piTop = piDF.sort_values("absImpact", ascending=False).head(20)
    piTop["model"] = f"{modelName} (permutation)"
    impactFrames.append(piTop[["model", "feature", "impact", "absImpact"]])

impactDF = pd.concat(impactFrames, ignore_index=True)
display(impactDF.head(50))

fig, axes = plt.subplots(1, 3, figsize=(18, 8))
for ax, mdl in zip(axes, impactDF["model"].unique()):
    sub = impactDF[impactDF["model"] == mdl].sort_values("absImpact", ascending=True).tail(15)
    sns.barplot(data=sub, x="absImpact", y="feature", ax=ax, palette="cubehelix", hue="feature", legend=False)
    ax.set_xlabel("Absolute Impact")
    ax.set_ylabel("Feature")
    ax.set_title(mdl)
plt.tight_layout()
plt.show()

compareDF = evalDF.merge(summaryDF[['model', 'bestTimeSec']], on='model', how='left')
display(compareDF.sort_values('testRmse'))

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
sns.barplot(data=compareDF, x='model', y='testRmse', ax=axes[0], palette='Blues')
axes[0].set_title('Test RMSE')
axes[0].set_xlabel("Model")
axes[0].set_ylabel("Test RMSE")
axes[0].tick_params(axis='x', rotation=15)

sns.barplot(data=compareDF, x='model', y='testR2', ax=axes[1], palette='Greens')
axes[1].set_title('Test $R^2$')
axes[0].set_xlabel("Model")
axes[1].set_ylabel("Test R2")
axes[1].tick_params(axis='x', rotation=15)

sns.barplot(data=compareDF, x='model', y='bestTimeSec', ax=axes[2], palette='Reds')
axes[2].set_title('Hyperparameter Search Time (sec)')
axes[0].set_xlabel("Model")
axes[1].set_ylabel("Search Time")
axes[2].tick_params(axis='x', rotation=15)

plt.tight_layout()
plt.show()

TEST_DATASET_PATH = "/kaggle/input/competitions/house-prices-advanced-regression-techniques/train.csv"
testDF = pd.read_csv(TEST_DATASET_PATH)
submissionIds = testDF["Id"].copy()
xProd = testDF.drop(columns=["Id"])

svrModel = bestModels["SVR"]
predLog = svrModel.predict(xProd)
predPrice = np.expm1(predLog)

submissionDF = pd.DataFrame({
    "Id": submissionIds,
    "SalePrice": predPrice
})

submissionDF.to_csv("/kaggle/working/submission_svr.csv", index=False)