spb/wp3_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3#4"""Spatial and robustness analysis (refactor of the original run_v3_spatial.py).56Produces results/v3_spatial_results.pkl with:7  1. ZIP3 fixed-effects OLS (in/out-of-sample R2)8  2. Moran's I on OLS residuals (3 random subsamples of 5,000, KNN k=8)9  3. XGBoost with latitude/longitude10  4. XGBoost without any geographic features11  5. Ablation analysis over 6 progressive feature sets1213Usage:  python scripts/01_spatial_robustness.py [--output results/v3_spatial_results.pkl]14"""1516import argparse17import pickle18import sys19import time20import warnings21from pathlib import Path2223import numpy as np24import pandas as pd25from sklearn.linear_model import LinearRegression26from sklearn.metrics import mean_squared_error, r2_score27from sklearn.model_selection import train_test_split2829sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))30from wp3 import config, data31from wp3.models import train_xgb3233warnings.filterwarnings("ignore")343536def zip3_fixed_effects_ols(df, X_full, y, idx_train, idx_test):37    """OLS with 3-digit-ZIP fixed effects, evaluated on the random holdout."""38    zip3 = df["address_zipcode"].fillna("000").str[:3]39    zip3_dummies = pd.get_dummies(zip3, prefix="zip3", drop_first=True).astype(np.int8)40    X_zip3 = pd.concat([X_full, zip3_dummies], axis=1)4142    ols = LinearRegression(n_jobs=-1)43    ols.fit(X_zip3.iloc[idx_train].values, y[idx_train])44    y_pred_train = ols.predict(X_zip3.iloc[idx_train].values)45    y_pred_test = ols.predict(X_zip3.iloc[idx_test].values)4647    return {48        "r2_insample": r2_score(y[idx_train], y_pred_train),49        "r2_outsample": r2_score(y[idx_test], y_pred_test),50        "rmse_outsample": float(np.sqrt(mean_squared_error(y[idx_test], y_pred_test))),51        "n_zip3": int(zip3.nunique()),52        "n_features": X_zip3.shape[1],53    }545556def morans_i_robustness(df, X_full, y, idx_train, n_subsamples=3, size=5000, k=8):57    """Moran's I on baseline-OLS residuals over independent random subsamples."""58    from esda.moran import Moran59    from libpysal.weights import KNN6061    ols = LinearRegression(n_jobs=-1)62    ols.fit(X_full.iloc[idx_train].values, y[idx_train])63    residuals = y - ols.predict(X_full.values)6465    coords = df[["latitude", "longitude"]].values66    values, pvalues = [], []67    for i in range(n_subsamples):68        rng = np.random.RandomState(config.RANDOM_STATE + i)69        idx_sub = rng.choice(len(df), size=size, replace=False)70        w = KNN.from_array(coords[idx_sub], k=k)71        w.transform = "r"72        mi = Moran(residuals[idx_sub], w)73        values.append(mi.I)74        pvalues.append(mi.p_sim)75        print(f"  Subsample {i + 1}: Moran's I = {mi.I:.4f}, p = {mi.p_sim:.4f}")7677    return {78        "values": values,79        "pvalues": pvalues,80        "mean": float(np.mean(values)),81        "std": float(np.std(values)),82        "pval_mean": float(np.mean(pvalues)),83        "n_subsamples": n_subsamples,84        "subsample_size": size,85        "k_neighbors": k,86    }878889def main():90    parser = argparse.ArgumentParser(description=__doc__)91    parser.add_argument("--output", type=Path,92                        default=config.RESULTS_DIR / "v3_spatial_results.pkl")93    args = parser.parse_args()9495    t0 = time.time()96    print("Loading analytical sample...")97    df = data.load_analytical_sample()98    X_full, parts = data.build_feature_matrix(df)99    y = df["ln_price"].values100    print(f"  {df.shape[0]:,} rows, {X_full.shape[1]} regressors")101102    idx_train, idx_test = train_test_split(103        np.arange(len(df)), test_size=0.2, random_state=config.RANDOM_STATE)104    geo_mask = df["address_state"].isin(config.HOLDOUT_STATES).values105    idx_geo_train, idx_geo_test = np.where(~geo_mask)[0], np.where(geo_mask)[0]106107    results = {}108109    print("\n1. ZIP3 fixed-effects OLS...")110    results["zip3_fe_ols"] = zip3_fixed_effects_ols(df, X_full, y, idx_train, idx_test)111    print(f"   R2 in={results['zip3_fe_ols']['r2_insample']:.4f} "112          f"out={results['zip3_fe_ols']['r2_outsample']:.4f}")113114    print("\n2. Moran's I robustness...")115    results["morans_i_robustness"] = morans_i_robustness(df, X_full, y, idx_train)116117    print("\n3. XGBoost with latitude/longitude...")118    X_latlon = pd.concat([X_full, df[["latitude", "longitude"]]], axis=1)119    results["xgb_with_latlon"] = train_xgb(120        X_latlon.values.astype(np.float32), y,121        idx_train, idx_test, idx_geo_train, idx_geo_test)122    print(f"   R2 random={results['xgb_with_latlon']['r2_random']:.4f} "123          f"geo={results['xgb_with_latlon']['r2_geo']:.4f}")124125    print("\n4. XGBoost without geographic features...")126    X_nogeo = pd.concat(127        [df[config.ALL_BASE_FEATS],128         parts["cat_dummies"][parts["pure_categorical_dummies"]]], axis=1)129    results["xgb_no_geography"] = train_xgb(130        X_nogeo.values.astype(np.float32), y,131        idx_train, idx_test, idx_geo_train, idx_geo_test)132    print(f"   R2 random={results['xgb_no_geography']['r2_random']:.4f} "133          f"geo={results['xgb_no_geography']['r2_geo']:.4f}")134135    print("\n5. Ablation analysis...")136    ablation_specs = [137        ("1_structural", config.STRUCTURAL_FEATS, []),138        ("2_plus_lot", config.STRUCTURAL_FEATS + config.LOT_FEATS, []),139        ("3_plus_amenities",140         config.STRUCTURAL_FEATS + config.LOT_FEATS + config.AMENITY_FEATS, []),141        ("4_plus_neighborhood",142         config.STRUCTURAL_FEATS + config.LOT_FEATS + config.AMENITY_FEATS143         + config.NEIGHBORHOOD_FEATS, []),144        ("5_plus_market",145         config.STRUCTURAL_FEATS + config.LOT_FEATS + config.AMENITY_FEATS146         + config.NEIGHBORHOOD_FEATS + config.MARKET_FEATS, []),147        ("6_full_model", config.ALL_BASE_FEATS, list(parts["cat_dummies"].columns)),148    ]149    ablation = {}150    for name, feats, dummy_cols in ablation_specs:151        X_abl = (pd.concat([df[feats], parts["cat_dummies"][dummy_cols]], axis=1)152                 if dummy_cols else df[feats])153        res = train_xgb(X_abl.values.astype(np.float32), y,154                        idx_train, idx_test, idx_geo_train, idx_geo_test)155        res["features"] = list(X_abl.columns)156        ablation[name] = res157        print(f"   {name}: R2 random={res['r2_random']:.4f} geo={res['r2_geo']:.4f}")158    results["ablation"] = ablation159160    args.output.parent.mkdir(parents=True, exist_ok=True)161    with open(args.output, "wb") as f:162        pickle.dump(results, f)163    print(f"\nSaved to {args.output}  ({time.time() - t0:.0f}s total)")164165166if __name__ == "__main__":167    main()168