Repository Public Research preview
1#!/usr/bin/env python32# -*- coding: utf-8 -*-3"""Révision scientifique 2026-10-17 (editorial board) — réanalyses locales, sans appel API.45 (a) langue de la description (FR / EN / mixte, règle de stop-words comme WP 2) et indicatrice de langue ;6 (b) fuite d'information : prix demandé et évaluation municipale présents dans le texte lu par le LLM7 (100 % des descriptions DuProprio contiennent « Asking Price $… ») -> gains de la couche texte par source,8 par présence d'un montant en dollars ;9 (c) effets fixes plus fins (RTA = 3 premiers caractères du code postal) et validation spatiale10 (GroupKFold par RTA) vs 5-fold aléatoire ;11 (d) bootstrap en grappes (municipalités) des coefficients IA du modèle A ;12 (e) accord inter-sources des extractions (LLM vs mots-clés / champs structurés ; vision vs texte / structuré) ;13 (f) sensibilité du blend au poids w et vérification du 50/50.1415Résultats -> results/revision_2026-10-17.json ; les tables LaTeX (paper/tables/rob_*.tex, agreement.tex) sont ensuite16générées par scripts/19_make_revision_tables.py à partir du JSON.17Exécution : ~3 min sur un Mac Apple Silicon M-series (9 validations croisées GBM, 299 réplications bootstrap).18"""19import json, re, time, warnings20import numpy as np, pandas as pd21import statsmodels.formula.api as smf22from sklearn.linear_model import LinearRegression23from sklearn.ensemble import HistGradientBoostingRegressor24from sklearn.model_selection import KFold, GroupKFold, cross_val_predict25from sklearn.metrics import r2_score26warnings.filterwarnings("ignore")27T0 = time.time()2829# pipeline v5 jusqu'avant la CV : df, y, data, specs, facade, AMEN, rv, gv, vv, NUM_v4, NUM_v5, cats30exec(open("scripts/09_v5_ml_text_vision.py").read().split("kf=KFold")[0])31amen_ = list(AMEN.keys())32y32 = y.copy() # float32, exactly as in scripts 09/12 (keeps the GBM fits identical to the paper's tables)33y = y.astype(np.float64)34OUT = {}3536# ------------------------------------------------------------------ texte, langue, code postal37meta = {}38for d in data:39 uid = f"{d.get('source')}:{d.get('listingId')}"40 meta[uid] = ((d.get("description") or ""), ((d.get("address") or {}).get("postalCode") or ""))41desc = df["uid"].map(lambda u: meta[u][0])42pc = df["uid"].map(lambda u: meta[u][1]).str.upper().str.replace(" ", "", regex=False)43df["fsa"] = pc.str[:3].where(pc.str.len() >= 3, "")4445EN = set("the and with this that for from has have is are of to in on your you its it be by".split())46FR = set("le la les et avec cette ce pour des une un du au aux est sont vous votre dans sur qui que".split())47def detect_language(t):48 w = re.findall(r"[a-zàâçéèêëîïôûùüÿœ']+", t.lower())49 e = sum(1 for x in w if x in EN); f = sum(1 for x in w if x in FR)50 if e == 0 and f == 0: return "mixed"51 if e >= 1.5 * f: return "en"52 if f >= 1.5 * e: return "fr"53 return "mixed"54df["lang"] = [detect_language(t) for t in desc]55df["lang_en"] = (df["lang"] == "en").astype(int); df["lang_mixed"] = (df["lang"] == "mixed").astype(int)56df["desc_len"] = desc.str.len().clip(upper=6000) / 1000.05758DOLLAR = re.compile(r"\d[\d\s,.]{3,}\s?\$|\$\s?\d[\d,.]{3,}")59ASSESS = re.compile(r"(évaluation municipale|evaluation municipale|valeur municipale|municipal assessment|municipal (evaluation|valuation)|évalu[ée]e? à|assessed at)", re.I)60ASKING = re.compile(r"Asking Price\s*\$?\s*([\d,]+)", re.I)61df["price_in_text"] = [int(bool(DOLLAR.search(t))) for t in desc]62df["assess_in_text"] = [int(bool(ASSESS.search(t))) for t in desc]63ask = [ASKING.search(t) for t in desc]64df["asking_in_text"] = [int(m is not None) for m in ask]65eq = sum(1 for m, p in zip(ask, df["price"]) if m and int(m.group(1).replace(",", "")) == int(p))66OUT["text_leakage"] = {67 "n": int(len(df)),68 "asking_price_in_text": int(df["asking_in_text"].sum()),69 "asking_price_equals_list_price": int(eq),70 "asking_by_source": df.groupby("source")["asking_in_text"].sum().astype(int).to_dict(),71 "dollar_amount_in_text": int(df["price_in_text"].sum()),72 "dollar_by_source": df.groupby("source")["price_in_text"].mean().round(4).to_dict(),73 "assessment_mention": int(df["assess_in_text"].sum()),74 "assessment_by_source": df.groupby("source")["assess_in_text"].mean().round(4).to_dict(),75 "realtor_dollar_share": float(df.loc[df.source == "realtor", "price_in_text"].mean()),76}77print("leakage:", OUT["text_leakage"])7879lang_tab = (df.assign(lnp=np.log(df["price"]))80 .groupby(["source", "lang"]).agg(n=("price", "size"), median_price=("price", "median"),81 mean_lnp=("lnp", "mean"), desc_len=("desc_len", "mean"))82 .reset_index())83OUT["language"] = {"by_source": lang_tab.to_dict(orient="records"),84 "share_en_overall": float(df["lang_en"].mean()),85 "share_en_realtor": float(df.loc[df.source == "realtor", "lang_en"].mean()),86 "share_by_lang": df["lang"].value_counts(normalize=True).round(4).to_dict()}87print(lang_tab)8889# ------------------------------------------------------------------ coverage of AI layers90df["has_specs"] = df["uid"].map(lambda u: int(u in specs)); df["has_facade"] = df["uid"].map(lambda u: int(u in facade))91OUT["coverage"] = {"specs": int(df.has_specs.sum()), "facade": int(df.has_facade.sum()),92 "is_exterior": int(df.v_is_ext.sum()), "n": int(len(df))}93print("coverage:", OUT["coverage"])9495# ------------------------------------------------------------------ OLS helpers96df_cv = df.copy() # pristine copy (NaN kept) for the tree models; the OLS block below imputes df in place97df["ln_price"] = np.log(df["price"])98df["ln_sqft"] = np.log(df["sqft_final"].fillna(df.groupby("ptype")["sqft_final"].transform("median")).fillna(df["sqft_final"].median()))99df["ln_lot"] = np.log(df["lot_sqft"].fillna(df["lot_sqft"].median())); df["has_lot"] = df["lot_sqft"].notna().astype(int)100df["has_age"] = df["age"].notna().astype(int)101for c in ["beds", "baths", "storeys", "age", "age2"]:102 df[c] = df[c].fillna(df.groupby("ptype")[c].transform("median")).fillna(df[c].median())103df["parking"] = df["parking"].fillna(0)104vc = df["fsa"].value_counts(); keepf = set(vc[(vc >= 20)].index) - {""}105df["fsa_fe"] = df["fsa"].where(df["fsa"].isin(keepf), "__other__")106OUT["fsa"] = {"n_fsa_raw": int(df.loc[df.fsa != "", "fsa"].nunique()), "n_fsa_fe": int(df["fsa_fe"].nunique()),107 "n_city_fe": int(df["city_fe"].nunique())}108print("FSA:", OUT["fsa"])109110CORE = ("ln_lot+has_lot+beds+baths+parking+storeys+age+age2+has_age+C(source)+C(ptype_fe)+"111 + "+".join(amen_) + "+n_images_c+virtual_tour+" + "+".join(rv))112TEXT = "+".join(gv); VIS = "+".join(vv)113def form(fe="C(city_fe)", area=True, text=True, vis=True, extra=""):114 f = "ln_price ~ " + ("ln_sqft+" if area else "") + CORE115 if text: f += "+" + TEXT116 if vis: f += "+" + VIS117 if extra: f += "+" + extra118 return f + "+" + fe119KEYS = ["v_curb", "v_reno", "v_luxury", "v_garage_vis", "v_mat_brique", "v_mat_pierre", "v_mat_vinyle", "v_set_rural",120 "v_land", "g_cond", "g_renovated_kitchen", "g_hardwood", "g_income_property", "lang_en", "lang_mixed"]121def grab(m):122 return {k: {"b": float(m.params[k]), "se": float(m.bse[k]), "p": float(m.pvalues[k])} for k in KEYS if k in m.params.index}123def fit(f, d, groups):124 return smf.ols(f, data=d).fit(cov_type="cluster", cov_kwds={"groups": d[groups]})125def ablation(d, fe="C(city_fe)", area=True):126 r = {}127 for nm, (t, v) in {"base": (False, False), "plus_text": (True, False), "plus_text_vision": (True, True)}.items():128 r[nm] = float(smf.ols(form(fe, area, t, v), data=d).fit().rsquared)129 r["n"] = int(len(d)); r["d_text"] = r["plus_text"] - r["base"]; r["d_vision"] = r["plus_text_vision"] - r["plus_text"]130 return r131132dfA = df[df["has_sqft"] == 1].copy()133134# ------------------------------------------------------------------ (a) language indicator135mA = fit(form(), dfA, "city_fe")136mA_lang = fit(form(extra="lang_en+lang_mixed+desc_len"), dfA, "city_fe")137mB = fit(form(area=False), df, "city_fe")138mB_lang = fit(form(area=False, extra="lang_en+lang_mixed+desc_len"), df, "city_fe")139OUT["ols"] = {"A": {"r2": mA.rsquared, "n": int(mA.nobs), "coef": grab(mA)},140 "A_lang": {"r2": mA_lang.rsquared, "n": int(mA_lang.nobs), "coef": grab(mA_lang)},141 "B": {"r2": mB.rsquared, "n": int(mB.nobs), "coef": grab(mB)},142 "B_lang": {"r2": mB_lang.rsquared, "n": int(mB_lang.nobs), "coef": grab(mB_lang)}}143print("lang_en A:", OUT["ols"]["A_lang"]["coef"]["lang_en"], "| B:", OUT["ols"]["B_lang"]["coef"]["lang_en"])144# raw language gap (no controls, city FE only)145m_gap = smf.ols("ln_price ~ lang_en + lang_mixed + C(source) + C(city_fe)", data=df).fit()146OUT["language"]["raw_gap_en_cityFE"] = float(m_gap.params["lang_en"])147OUT["language"]["raw_gap_en_nocontrols"] = float(smf.ols("ln_price ~ lang_en + lang_mixed", data=df).fit().params["lang_en"])148149# in-sample ablation by language / source / price-in-text (model A, city FE)150OUT["ablation_A"] = {"all": ablation(dfA),151 "fr": ablation(dfA[dfA.lang == "fr"]), "en": ablation(dfA[dfA.lang == "en"]),152 "mixed": ablation(dfA[dfA.lang == "mixed"]),153 "realtor": ablation(dfA[dfA.source == "realtor"]), "duproprio": ablation(dfA[dfA.source == "duproprio"]),154 "realtor_no_dollar": ablation(dfA[(dfA.source == "realtor") & (dfA.price_in_text == 0)]),155 "realtor_dollar": ablation(dfA[(dfA.source == "realtor") & (dfA.price_in_text == 1)]),156 "with_lang": {k: float(v) for k, v in zip(["base", "plus_text", "plus_text_vision"],157 [smf.ols(form(text=t, vis=v, extra="lang_en+lang_mixed+desc_len"), data=dfA).fit().rsquared158 for t, v in [(False, False), (True, False), (True, True)]])}}159for k, v in OUT["ablation_A"].items(): print("ablation", k, {a: round(b, 4) for a, b in v.items()})160161# ------------------------------------------------------------------ (c) FSA fixed effects162mA_fsa = fit(form(fe="C(fsa_fe)"), dfA, "fsa_fe"); mB_fsa = fit(form(fe="C(fsa_fe)", area=False), df, "fsa_fe")163mA_fsa_lang = fit(form(fe="C(fsa_fe)", extra="lang_en+lang_mixed+desc_len"), dfA, "fsa_fe")164OUT["ols"]["A_fsa"] = {"r2": mA_fsa.rsquared, "n": int(mA_fsa.nobs), "coef": grab(mA_fsa)}165OUT["ols"]["A_fsa_lang"] = {"r2": mA_fsa_lang.rsquared, "n": int(mA_fsa_lang.nobs), "coef": grab(mA_fsa_lang)}166OUT["ols"]["B_fsa"] = {"r2": mB_fsa.rsquared, "n": int(mB_fsa.nobs), "coef": grab(mB_fsa)}167OUT["ablation_A"]["fsa_fe"] = ablation(dfA, fe="C(fsa_fe)")168# Realtor-only with FSA FE + language (the cleanest specification: no price in text)169dR = dfA[dfA.source == "realtor"]170mA_realtor = fit(form(fe="C(fsa_fe)", extra="lang_en+lang_mixed+desc_len"), dR, "fsa_fe")171OUT["ols"]["A_realtor_fsa_lang"] = {"r2": mA_realtor.rsquared, "n": int(mA_realtor.nobs), "coef": grab(mA_realtor)}172print("FSA A R2=%.4f B R2=%.4f realtor-FSA-lang R2=%.4f" % (mA_fsa.rsquared, mB_fsa.rsquared, mA_realtor.rsquared))173174# ------------------------------------------------------------------ (d) cluster (pairs) bootstrap, model A175import patsy176yA, XA = patsy.dmatrices(form(), data=dfA, return_type="dataframe")177XA = XA.values; yA = yA.values.ravel(); cols = list(patsy.dmatrices(form(), data=dfA, return_type="dataframe")[1].columns)178groups = dfA["city_fe"].values; ug = np.unique(groups); gidx = {g: np.where(groups == g)[0] for g in ug}179rng = np.random.RandomState(20261017); B = 299; keep_ix = [cols.index(k) for k in KEYS if k in cols]180boot = np.full((B, len(keep_ix)), np.nan)181for b in range(B):182 samp = rng.choice(ug, size=len(ug), replace=True)183 rows = np.concatenate([gidx[g] for g in samp])184 Xb = XA[rows]; yb = yA[rows]185 nz = Xb.std(0) > 0; nz[0] = True186 beta = np.full(XA.shape[1], np.nan)187 beta[nz] = np.linalg.lstsq(Xb[:, nz], yb, rcond=None)[0]188 boot[b] = beta[keep_ix]189OUT["cluster_bootstrap_A"] = {"B": B, "n_clusters": int(len(ug)),190 "ci95": {k: [float(np.nanpercentile(boot[:, j], 2.5)), float(np.nanpercentile(boot[:, j], 97.5))]191 for j, k in enumerate([k for k in KEYS if k in cols])},192 "se": {k: float(np.nanstd(boot[:, j], ddof=1)) for j, k in enumerate([k for k in KEYS if k in cols])}}193print("bootstrap CI:", {k: [round(a, 4) for a in v] for k, v in OUT["cluster_bootstrap_A"]["ci95"].items()})194195# ------------------------------------------------------------------ (e) agreement between extraction sources196def agree(a, b):197 a = np.asarray(a).astype(int); b = np.asarray(b).astype(int)198 po = float(np.mean(a == b)); pe = float(np.mean(a) * np.mean(b) + (1 - np.mean(a)) * (1 - np.mean(b)))199 kappa = (po - pe) / (1 - pe) if pe < 1 else np.nan200 tp = int(((a == 1) & (b == 1)).sum()); fp = int(((a == 1) & (b == 0)).sum()); fn = int(((a == 0) & (b == 1)).sum())201 prec = tp / (tp + fp) if tp + fp else np.nan; rec = tp / (tp + fn) if tp + fn else np.nan202 return {"agreement": po, "kappa": kappa, "precision_vs_ref": prec, "recall_vs_ref": rec,203 "share_a": float(np.mean(a)), "share_ref": float(np.mean(b)), "n": int(len(a))}204ds = df[df.has_specs == 1]205OUT["agreement"] = {"llm_vs_keyword": {206 "garage": agree(ds.g_garage, ds.garage), "pool": agree(ds.g_pool, ds.pool), "waterfront": agree(ds.g_waterfront, ds.waterfront),207 "fireplace": agree(ds.g_fireplace, ds.fireplace), "central_ac": agree(ds.g_central_ac, ds.central_air),208 "finished_basement": agree(ds.g_finished_basement, ds.finished_bsmt)}}209dv = df[(df.has_facade == 1) & (df.v_is_ext == 1) & (df.has_specs == 1)]210OUT["agreement"]["vision_vs_text"] = {"garage_visible_vs_llm_garage": agree(dv.v_garage_vis, dv.g_garage),211 "garage_visible_vs_keyword_garage": agree(dv.v_garage_vis, dv.garage),212 "waterfront_setting_vs_keyword": agree(dv.v_set_bord_eau, dv.waterfront)}213# vision storeys vs structured storeys ; vision building type vs listing type214sv = dv["uid"].map(lambda u: facade[u].get("storeys_visual")); st = dv["storeys"]215m = sv.isin(["1", "1.5", "2", "3+"]) & st.notna()216def st_cat(x): return "1" if x < 1.25 else "1.5" if x < 1.75 else "2" if x < 2.75 else "3+"217OUT["agreement"]["vision_storeys_vs_structured"] = {"n": int(m.sum()),218 "exact": float(np.mean(sv[m].values == st[m].map(st_cat).values)),219 "within_half_storey": float(np.mean(np.abs(sv[m].map({"1": 1, "1.5": 1.5, "2": 2, "3+": 3}).values - st[m].clip(upper=3).values) <= 0.5))}220bt = dv["uid"].map(lambda u: facade[u].get("building_type_visual"))221row_types = {"Row / Townhouse", "Townhouse"}222is_row = dv["ptype"].isin(row_types).astype(int); vis_row = (bt == "rangee").astype(int)223OUT["agreement"]["vision_rowhouse_vs_listing_type"] = agree(vis_row, is_row)224OUT["agreement"]["vision_building_type_dist"] = bt.value_counts(normalize=True).round(3).to_dict()225print("agreement:", json.dumps(OUT["agreement"], indent=0, default=float)[:1500])226227# ------------------------------------------------------------------ CV: random vs spatial (RTA) holdout, by subgroup228kf = KFold(n_splits=5, shuffle=True, random_state=42)229grp = df["fsa"].where(df["fsa"] != "", df["city_fe"]).values230gkf = GroupKFold(n_splits=5)231Xcat = pd.get_dummies(df[cats], dtype=float).reset_index(drop=True)232dfi = df.copy()233for c in ["sqft_final", "lot_sqft", "beds", "baths", "storeys", "age", "age2"]:234 dfi[c] = dfi.groupby("ptype")[c].transform(lambda s: s.fillna(s.median())); dfi[c] = dfi[c].fillna(dfi[c].median())235dfi["parking"] = dfi["parking"].fillna(0)236def stats(pred, mask=None):237 mask = np.ones(len(y), bool) if mask is None else mask238 e = np.abs(np.exp(pred[mask] - y[mask]) - 1)239 return {"r2": float(r2_score(y[mask], pred[mask])), "mdape": float(np.median(e) * 100), "n": int(mask.sum())}240def cv_ols(numcols, cv, groups=None):241 X = pd.concat([dfi[numcols].reset_index(drop=True), Xcat], axis=1); X = X.fillna(X.median())242 return cross_val_predict(LinearRegression(), X, y, cv=cv, groups=groups, n_jobs=-1)243def cv_gbm(numcols, cv, groups=None):244 X = pd.concat([df_cv[numcols].reset_index(drop=True), Xcat], axis=1)245 m = HistGradientBoostingRegressor(max_iter=600, learning_rate=0.05, max_leaf_nodes=63, l2_regularization=1.0, random_state=42)246 return cross_val_predict(m, X, y32, cv=cv, groups=groups, n_jobs=-1).astype(np.float64)247L1 = ["sqft_final", "has_sqft", "lot_sqft", "beds", "baths", "parking", "storeys", "age", "age2"] + amen_248L2 = L1 + ["n_images_c", "virtual_tour"] + rv; L3 = L2 + gv; L4 = L3 + vv249L4L = L4 + ["lang_en", "lang_mixed", "desc_len"]250masks = {"all": np.ones(len(y), bool), "fr": (df.lang == "fr").values, "en": (df.lang == "en").values,251 "mixed": (df.lang == "mixed").values, "realtor": (df.source == "realtor").values,252 "duproprio": (df.source == "duproprio").values,253 "realtor_no_dollar": ((df.source == "realtor") & (df.price_in_text == 0)).values,254 "realtor_dollar": ((df.source == "realtor") & (df.price_in_text == 1)).values}255OUT["cv"] = {"random": {}, "spatial_fsa": {}}256for lname, cols in [("L2", L2), ("L3", L3), ("L4", L4), ("L4+lang", L4L)]:257 po = cv_ols(cols, kf); pg = cv_gbm(cols, kf)258 OUT["cv"]["random"][lname] = {"ols": {k: stats(po, m) for k, m in masks.items()}, "gbm": {k: stats(pg, m) for k, m in masks.items()}}259 print(f"random {lname}: OLS {OUT['cv']['random'][lname]['ols']['all']} GBM {OUT['cv']['random'][lname]['gbm']['all']} [{time.time()-T0:.0f}s]")260for lname, cols in [("L2", L2), ("L3", L3), ("L4", L4)]:261 po = cv_ols(cols, gkf, grp); pg = cv_gbm(cols, gkf, grp)262 OUT["cv"]["spatial_fsa"][lname] = {"ols": {k: stats(po, m) for k, m in masks.items()}, "gbm": {k: stats(pg, m) for k, m in masks.items()}}263 print(f"spatial {lname}: OLS {OUT['cv']['spatial_fsa'][lname]['ols']['all']} GBM {OUT['cv']['spatial_fsa'][lname]['gbm']['all']} [{time.time()-T0:.0f}s]")264OUT["cv"]["spatial_groups"] = int(len(np.unique(grp)))265266# ------------------------------------------------------------------ (f) blend weights267g5 = cv_gbm(NUM_v5, kf); b5 = np.load("data/processed/v5_blend_oos.npy").astype(np.float64); n5 = 2 * b5 - g5268g4 = cv_gbm(NUM_v4, kf); n4 = np.load("data/processed/nn_robust_oos.npy").astype(np.float64)269OUT["blend"] = {"gbm_v5_r2": float(r2_score(y, g5)), "nn_v5_r2_recovered": float(r2_score(y, n5)),270 "blend_v5_r2_file": float(r2_score(y, b5)), "gbm_v4_r2": float(r2_score(y, g4)), "nn_v4_r2": float(r2_score(y, n4)),271 "grid_v5": {str(w): float(r2_score(y, w * g5 + (1 - w) * n5)) for w in [0.0, 0.3, 0.4, 0.5, 0.6, 0.7, 1.0]},272 "grid_v4": {str(w): float(r2_score(y, w * g4 + (1 - w) * n4)) for w in [0.0, 0.3, 0.4, 0.5, 0.6, 0.7, 1.0]},273 "blend_v5_by_group": {k: stats(b5, m) for k, m in masks.items()}}274# optimal stacking weight by nested CV (weight chosen on 4 folds, applied to the 5th)275wgrid = np.linspace(0, 1, 21); pred_nested = np.zeros_like(y)276for tr, te in kf.split(y):277 wbest = wgrid[np.argmax([r2_score(y[tr], w * g5[tr] + (1 - w) * n5[tr]) for w in wgrid])]278 pred_nested[te] = wbest * g5[te] + (1 - wbest) * n5[te]279OUT["blend"]["nested_weight_r2"] = float(r2_score(y, pred_nested))280print("blend:", {k: v for k, v in OUT["blend"].items() if k != "blend_v5_by_group"})281282json.dump(OUT, open("results/revision_2026-10-17.json", "w"), indent=1, default=float)283print(f"-> results/revision_2026-10-17.json ({time.time()-T0:.0f}s)")284print("Next: python3 scripts/19_make_revision_tables.py")285