Kaggle Competition - Flood Prediction EDA

# Standard library
import datetime
from collections import defaultdict

# Data manipulation
import numpy as np
import pandas as pd

# Visualisation
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
import seaborn as sns

# Machine learning and modelling
from sklearn.base import clone
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.ensemble import ExtraTreesRegressor, RandomForestRegressor  # Optional: if used later
from sklearn.pipeline import make_pipeline, Pipeline
from sklearn.model_selection import KFold
from sklearn.preprocessing import StandardScaler, PolynomialFeatures, SplineTransformer
from sklearn.metrics import r2_score
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE

# UMAP for dimensionality reduction
import umap

# Statistical modelling
import statsmodels.api as sm

train = pd.read_csv('train.csv', index_col='id')
test = pd.read_csv('test.csv', index_col='id')

initial_features = list(test.columns)
train

bin_edges = np.arange(0.2825, 0.7300, 0.005)

plt.figure(figsize=(8, 3))
plt.hist(train['FloodProbability'], bins=bin_edges, density=True, edgecolor='black', linewidth=0.05)
plt.title('Distribution of FloodProbability')
plt.xlabel('FloodProbability')
plt.ylabel('Probability Mass')
plt.grid(axis='y', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator

# Plot feature distributions for train vs test
fig, axs = plt.subplots(5, 4, figsize=(14, 12))  # Slightly wider for better spacing

for col, ax in zip(initial_features, axs.ravel()):
    # Relative frequency in training set
    train_dist = train[col].value_counts(normalize=True).sort_index()
    ax.bar(train_dist.index, train_dist.values, label='Train', alpha=1.0)

    # Relative frequency in test set
    test_dist = test[col].value_counts(normalize=True).sort_index()
    ax.bar(test_dist.index, test_dist.values, label='Test', alpha=0.6)

    ax.set_title(col)
    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
    ax.set_ylabel('Proportion')

# Add legend only once, outside the grid
handles, labels = ax.get_legend_handles_labels()
fig.legend(handles, labels, loc='upper right', fontsize='medium')

plt.tight_layout()
plt.show()

# Prepare the list of features including the target
corr_features = initial_features + ['FloodProbability']

# Compute correlation matrix
correlation_matrix = np.corrcoef(train[corr_features], rowvar=False)

# Plot heatmap
plt.figure(figsize=(11, 11))
sns.heatmap(
    correlation_matrix,
    center=0,
    cmap='coolwarm',
    annot=True,
    fmt='.1f',
    xticklabels=corr_features,
    yticklabels=corr_features,
    square=True,
    linewidths=0.5,
    cbar_kws={'shrink': 0.8}
)
plt.title('Correlation Matrix of Features and Target', fontsize=14)
plt.tight_layout()
plt.show()

# Fit PCA on the feature set
pca = PCA()
pca.fit(train[initial_features])

# Plot cumulative explained variance
plt.figure(figsize=(4, 3))
plt.plot(np.arange(1, len(initial_features) + 1), pca.explained_variance_ratio_.cumsum(), marker='o')
plt.title('Principal Component Analysis', fontsize=12)
plt.xlabel('Principal Component', fontsize=10)
plt.ylabel('Cumulative Explained Variance', fontsize=10)
plt.xticks(fontsize=9)
plt.yticks([0, 0.5, 1.0], fontsize=9)
plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# Function to plot the embedding
def plot_embedding(embedding, target, title):
    plt.figure(figsize=(6, 5))
    plt.scatter(
        embedding[:, 0],
        embedding[:, 1],
        s=2,
        c=target,
        cmap='coolwarm',
        alpha=0.6
    )
    plt.gca().set_aspect('equal', 'datalim')
    plt.title(title, fontsize=14)
    plt.xlabel('UMAP Dimension 1')
    plt.ylabel('UMAP Dimension 2')
    plt.colorbar(label='FloodProbability')
    plt.tight_layout()
    plt.show()

# Sample a subset for faster projection
train_sample = train.sample(10_000)

# Fit UMAP
reducer = umap.UMAP()
embedding = reducer.fit_transform(train_sample[initial_features])

# Plot the result
plot_embedding(
    embedding,
    train_sample['FloodProbability'],
    'Unsupervised UMAP Projection of the Training Dataset'
)

# Sample a smaller subset for faster computation
train_sample = train.sample(10000, random_state=42)

# Fit Supervised UMAP using regression target
reducer = umap.UMAP(
    n_neighbors=50,         # smaller neighborhood for faster computation
    min_dist=0.5,           # moderate compression
    target_metric='manhattan',
    target_weight=0.6,
)

# Perform the supervised projection
embedding = reducer.fit_transform(
    train_sample[initial_features],
    y=train_sample['FloodProbability']
)

# Plotting function
def plot_embedding(embedding, target, title):
    plt.figure(figsize=(6, 5))
    scatter = plt.scatter(
        embedding[:, 0],
        embedding[:, 1],
        s=2,
        c=target,
        cmap='coolwarm',
        alpha=0.7
    )
    plt.gca().set_aspect('equal', 'datalim')
    plt.title(title, fontsize=14)
    plt.xlabel('UMAP Dimension 1')
    plt.ylabel('UMAP Dimension 2')
    cbar = plt.colorbar(scatter)
    cbar.set_label('FloodProbability')
    plt.tight_layout()
    plt.show()

# Visualise the result
plot_embedding(
    embedding,
    train_sample['FloodProbability'],
    'Supervised UMAP Projection of the Training Dataset'
)

# Sample a subset for computational efficiency
train_sample = train.sample(20_000, random_state=42)

# Fit and transform with t-SNE
reducer = TSNE(random_state=42, perplexity=30, max_iter=1000, learning_rate='auto')
embedding = reducer.fit_transform(train_sample[initial_features])

# Plot the result
plot_embedding(
    embedding,
    train_sample['FloodProbability'],
    '(Unsupervised) t-SNE Projection of the Training Dataset'
)

from collections import defaultdict

# Configuration
kf = KFold(n_splits=5, shuffle=True, random_state=1)
SINGLE_FOLD = True       # Toggle for full vs. single-fold CV
COMPUTE_TEST_PRED = True # Toggle to compute test set predictions
oof = defaultdict(lambda: np.full_like(train.FloodProbability, np.nan, dtype=float))
test_pred = {}

def cross_validate(model, label, features=initial_features, n_repeats=1):
    """
    Evaluate a model using cross-validation.
    
    Parameters:
    - model: scikit-learn compatible regressor
    - label: string identifier for the model
    - features: list of feature names to use
    - n_repeats: number of training repeats with different random seeds

    Outputs:
    - Stores out-of-fold predictions in `oof[label]`
    - Stores averaged test predictions in `test_pred[label]` if enabled
    """
    start_time = datetime.datetime.now()
    scores = []
    oof_preds = np.full_like(train.FloodProbability, np.nan, dtype=float)

    for fold, (idx_tr, idx_va) in enumerate(kf.split(train)):
        X_tr = train.iloc[idx_tr][features]
        X_va = train.iloc[idx_va][features]
        y_tr = train.iloc[idx_tr].FloodProbability
        y_va = train.iloc[idx_va].FloodProbability

        y_pred = np.zeros_like(y_va, dtype=float)

        for i in range(n_repeats):
            m = clone(model)
            if n_repeats > 1:
                mm = m
                if isinstance(mm, Pipeline):
                    mm = mm[-1]
                if isinstance(mm, TransformedTargetRegressor):
                    mm = mm.regressor
                mm.set_params(random_state=i)

            m.fit(X_tr, y_tr)
            y_pred += m.predict(X_va)

        y_pred /= n_repeats
        score = r2_score(y_va, y_pred)
        print(f"# Fold {fold}: R² = {score:.5f}")
        scores.append(score)
        oof_preds[idx_va] = y_pred

        if SINGLE_FOLD:
            break

    elapsed_time = datetime.datetime.now() - start_time
    avg_score = np.mean(scores)
    print(f"{Fore.GREEN}# Overall R²: {avg_score:.5f} | Model: {label}"
          f"{' (single fold)' if SINGLE_FOLD else ''} "
          f"| Time: {int(np.round(elapsed_time.total_seconds() / 60))} min{Style.RESET_ALL}")
    oof[label] = oof_preds

    # Optional: generate test set predictions
    if COMPUTE_TEST_PRED:
        X_tr = train[features]
        y_tr = train.FloodProbability
        y_pred = np.zeros(len(test), dtype=float)

        for i in range(n_repeats):
            m = clone(model)
            if n_repeats > 1:
                mm = m
                if isinstance(mm, Pipeline):
                    mm = mm[-1]
                if isinstance(mm, TransformedTargetRegressor):
                    mm = mm.regressor
                mm.set_params(random_state=i)

            m.fit(X_tr, y_tr)
            y_pred += m.predict(test[features])

        y_pred /= n_repeats
        test_pred[label] = y_pred

# Standard Linear Regression
model = make_pipeline(
    StandardScaler(),
    LinearRegression()
)
cross_validate(model, 'LinearRegression')

# Polynomial Features + Ridge Regression
model = make_pipeline(
    StandardScaler(),
    PolynomialFeatures(degree=2, include_bias=False),
    Ridge()
)
cross_validate(model, 'Poly-Ridge')

# Spline Transformation + Ridge Regression
model = make_pipeline(
    StandardScaler(),
    SplineTransformer(n_knots=5, degree=3),
    Ridge()
)
cross_validate(model, 'Spline-Ridge')

import statsmodels.api as sm

# Add intercept term
X = sm.add_constant(train[initial_features])

# Fit OLS model and display summary
model = sm.OLS(train['FloodProbability'], X, missing='error')
results = model.fit()
results.summary()

# Compute approximation using summed features
approx_prediction = train[initial_features].sum(axis=1) * 0.00565 - 0.053

# Evaluate goodness of fit
score = r2_score(train['FloodProbability'], approx_prediction)
print(f"R² Score (approximate model): {score:.5f}")