spb/wp5_uqo Public
UQO Working Paper No. 5 — Airbnb, residential rents and housing market pressure.
TeX 53.4%
Python 46.5%
1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""309_ml_robustness.py4-------------------5Model 6: ML Robustness Models.6Predict log_rent using Airbnb exposure and property characteristics.7Compare OLS, LASSO, Elastic Net, Random Forest, and Gradient Boosting.8Compute SHAP values (if available) or sklearn feature importances.910Outputs11-------12- results/tables/ml_comparison.tex13- results/tables/ml_comparison.csv14- figures/ml_predicted_vs_actual.pdf15- figures/feature_importance.pdf16- figures/shap_summary.pdf (if shap available)17"""1819import sys20import warnings21from pathlib import Path2223sys.path.insert(0, str(Path(__file__).resolve().parent.parent))2425from src.plotting import use_publication_style2627use_publication_style()2829import matplotlib.pyplot as plt # noqa: E40230import numpy as np # noqa: E40231import pandas as pd # noqa: E40232from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor # noqa: E40233from sklearn.linear_model import ElasticNetCV, LassoCV, LinearRegression # noqa: E40234from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score # noqa: E40235from sklearn.model_selection import train_test_split # noqa: E40236from sklearn.preprocessing import StandardScaler # noqa: E4023738from src.config import ( # noqa: E40239 MERGED_ANALYSIS,40 TABLE_DIR,41 FIG_DIR,42 RANDOM_STATE,43 require,44)4546warnings.filterwarnings("ignore")4748PROCESSED_HINT = "Run scripts/04_merge_data.py first."4950FEATURES = [51 "airbnb_count_500m",52 "airbnb_density_500m",53 "mean_airbnb_price_500m",54 "share_entire_home_500m",55 "bedrooms",56 "bathrooms",57 "lat",58 "lon",59]60TARGET = "log_rent"6162# Linear models use standardized features; tree-based models use raw features63USE_SCALED = {"OLS", "LASSO", "Elastic Net"}646566def train_models(X_train, X_test, y_train, y_test, X_train_sc, X_test_sc):67 """Train the five models and return (results DataFrame, fitted dict)."""68 models = {69 "OLS": LinearRegression(),70 "LASSO": LassoCV(cv=5, random_state=RANDOM_STATE, max_iter=10_000),71 "Elastic Net": ElasticNetCV(cv=5, random_state=RANDOM_STATE, max_iter=10_000),72 "Random Forest": RandomForestRegressor(73 n_estimators=500, max_depth=15, random_state=RANDOM_STATE, n_jobs=-174 ),75 "Gradient Boosting": GradientBoostingRegressor(76 n_estimators=500, max_depth=5, learning_rate=0.05,77 random_state=RANDOM_STATE78 ),79 }8081 results = []82 fitted = {}8384 for name, model in models.items():85 print(f"Training {name} ...")86 Xtr = X_train_sc if name in USE_SCALED else X_train87 Xte = X_test_sc if name in USE_SCALED else X_test8889 model.fit(Xtr, y_train)90 fitted[name] = model9192 y_pred_train = model.predict(Xtr)93 y_pred_test = model.predict(Xte)9495 r2_train = r2_score(y_train, y_pred_train)96 r2_test = r2_score(y_test, y_pred_test)97 rmse_test = np.sqrt(mean_squared_error(y_test, y_pred_test))98 mae_test = mean_absolute_error(y_test, y_pred_test)99100 results.append({101 "Model": name,102 "R2_train": r2_train,103 "R2_test": r2_test,104 "RMSE_test": rmse_test,105 "MAE_test": mae_test,106 })107 print(f" R2 train={r2_train:.4f} test={r2_test:.4f} "108 f"RMSE={rmse_test:.4f} MAE={mae_test:.4f}")109110 return pd.DataFrame(results), fitted111112113def save_comparison_table(res_df: pd.DataFrame) -> None:114 """Write ml_comparison.csv and ml_comparison.tex."""115 res_df.to_csv(TABLE_DIR / "ml_comparison.csv", index=False)116117 latex_rows = []118 latex_rows.append(r"\begin{tabular}{lcccc}")119 latex_rows.append(r"\toprule")120 latex_rows.append(r"Model & $R^2$ (Train) & $R^2$ (Test) & RMSE (Test) & MAE (Test) \\")121 latex_rows.append(r"\midrule")122 for _, row in res_df.iterrows():123 latex_rows.append(124 f"{row['Model']} & {row['R2_train']:.4f} & {row['R2_test']:.4f} "125 f"& {row['RMSE_test']:.4f} & {row['MAE_test']:.4f} \\\\"126 )127 latex_rows.append(r"\bottomrule")128 latex_rows.append(r"\end{tabular}")129 latex_rows.append(r"\begin{tablenotes}\small")130 latex_rows.append(r"\item \textit{Notes:} OLS, LASSO, and Elastic Net use "131 r"standardized features. Tree-based models use raw features. "132 r"Train/test split is 80/20 with random\_state=42.")133 latex_rows.append(r"\end{tablenotes}")134135 (TABLE_DIR / "ml_comparison.tex").write_text("\n".join(latex_rows) + "\n",136 encoding="utf-8")137 print(f"\nSaved {TABLE_DIR / 'ml_comparison.tex'}")138 print(f"Saved {TABLE_DIR / 'ml_comparison.csv'}")139140141def plot_predicted_vs_actual(res_df, fitted, X_test, X_test_sc, y_test) -> None:142 """Scatter of predicted vs actual log rent for the best model by test R2."""143 best_name = res_df.loc[res_df["R2_test"].idxmax(), "Model"]144 best_model = fitted[best_name]145 Xte_best = X_test_sc if best_name in USE_SCALED else X_test146 y_pred_best = best_model.predict(Xte_best)147148 fig, ax = plt.subplots()149 ax.scatter(y_test, y_pred_best, alpha=0.3, s=8, edgecolors="none", c="#2166ac")150 mn, mx = min(y_test.min(), y_pred_best.min()), max(y_test.max(), y_pred_best.max())151 ax.plot([mn, mx], [mn, mx], "k--", lw=1, label="45-degree line")152 ax.set_xlabel("Actual log(rent)")153 ax.set_ylabel("Predicted log(rent)")154 ax.set_title(f"Predicted vs Actual — {best_name} "155 f"(test R²={res_df.loc[res_df['R2_test'].idxmax(), 'R2_test']:.4f})")156 ax.legend(loc="upper left", frameon=False)157 fig.savefig(FIG_DIR / "ml_predicted_vs_actual.pdf")158 plt.close(fig)159 print(f"Saved {FIG_DIR / 'ml_predicted_vs_actual.pdf'}")160161162def plot_importances(fitted, X_test) -> None:163 """SHAP summary + importance plots, with sklearn fallback."""164 shap_available = False165 try:166 import shap167 shap_available = True168 print("\nSHAP library found — computing SHAP values ...")169 except ImportError:170 print("\nSHAP not available — using sklearn feature importances instead.")171172 if shap_available:173 # Use the best tree model for SHAP; prefer Gradient Boosting, fall back RF174 for shap_model_name in ["Gradient Boosting", "Random Forest"]:175 if shap_model_name in fitted:176 break177 shap_model = fitted[shap_model_name]178 Xte_shap = X_test # tree models use raw features179180 try:181 explainer = shap.TreeExplainer(shap_model)182 shap_values = explainer.shap_values(Xte_shap)183184 # Summary bee-swarm plot185 plt.figure()186 shap.summary_plot(187 shap_values, Xte_shap, feature_names=FEATURES, show=False188 )189 plt.tight_layout()190 plt.savefig(FIG_DIR / "shap_summary.pdf")191 plt.close()192 print(f"Saved {FIG_DIR / 'shap_summary.pdf'}")193194 # Bar plot of mean |SHAP|195 plt.figure()196 shap.summary_plot(197 shap_values, Xte_shap, feature_names=FEATURES,198 plot_type="bar", show=False199 )200 plt.tight_layout()201 plt.savefig(FIG_DIR / "feature_importance.pdf")202 plt.close()203 print(f"Saved {FIG_DIR / 'feature_importance.pdf'}")204205 except Exception as e:206 print(f"SHAP computation failed ({e}); falling back to sklearn importances.")207 shap_available = False208209 if not shap_available:210 # Fallback: sklearn feature importances from tree models211 for imp_model_name in ["Gradient Boosting", "Random Forest"]:212 if imp_model_name in fitted and hasattr(fitted[imp_model_name],213 "feature_importances_"):214 break215 imp_model = fitted[imp_model_name]216 importances = imp_model.feature_importances_217 idx = np.argsort(importances)[::-1]218219 fig, ax = plt.subplots(figsize=(7, 5))220 ax.barh(221 range(len(FEATURES)),222 importances[idx[::-1]],223 color="#2166ac",224 edgecolor="white",225 )226 ax.set_yticks(range(len(FEATURES)))227 ax.set_yticklabels([FEATURES[i] for i in idx[::-1]])228 ax.set_xlabel("Feature Importance")229 ax.set_title(f"Feature Importance — {imp_model_name}")230 fig.savefig(FIG_DIR / "feature_importance.pdf")231 plt.close(fig)232 print(f"Saved {FIG_DIR / 'feature_importance.pdf'}")233234235def main() -> None:236 # 1. Load and prepare data237 print("Loading data ...")238 df = pd.read_parquet(require(MERGED_ANALYSIS, PROCESSED_HINT))239240 cols = FEATURES + [TARGET]241 df_ml = df[cols].dropna()242 print(f"Sample size after dropping NaN: {len(df_ml):,}")243244 X = df_ml[FEATURES].values245 y = df_ml[TARGET].values246247 # 2. Train / test split248 X_train, X_test, y_train, y_test = train_test_split(249 X, y, test_size=0.20, random_state=RANDOM_STATE250 )251 print(f"Train: {len(X_train):,} | Test: {len(X_test):,}")252253 # 3. Standardize features254 scaler = StandardScaler()255 X_train_sc = scaler.fit_transform(X_train)256 X_test_sc = scaler.transform(X_test)257258 # 4. Define and train models259 res_df, fitted = train_models(X_train, X_test, y_train, y_test,260 X_train_sc, X_test_sc)261262 # 5. Save comparison table263 save_comparison_table(res_df)264265 # 6. Predicted vs actual scatter (best model by test R²)266 plot_predicted_vs_actual(res_df, fitted, X_test, X_test_sc, y_test)267268 # 7. SHAP values or sklearn feature importances269 plot_importances(fitted, X_test)270271 print("\n=== 09_ml_robustness.py complete ===")272273274if __name__ == "__main__":275 main()276