1"""Phase C — effect-size harmonisation (protocol §3–4, CLAUDE.md §7). Every transformation is explicit and recorded in the2`formula` column; nothing is converted without the information the formula requires.34Families: pct_premium (binary attribute, log price), semi_elasticity (% change in price per unit of X, log price, X in level),5elasticity (log-log), monetary (level price, level X: currency per unit), partial_r (unit-free, from t and df).6"""7import math8import numpy as np, pandas as pd910AREA_TO_M2 = {"m2": 1.0, "sqm": 1.0, "sqft": 0.09290304, "ft2": 0.09290304, "100sqft": 9.290304, "1000sqft": 92.90304, "ha": 10000.0, "acre": 4046.8564224, "are": 100.0, "100m2": 100.0, "pyeong": 3.3058, "tsubo": 3.3058}11DIST_TO_KM = {"km": 1.0, "m": 0.001, "100m": 0.1, "mile": 1.609344, "miles": 1.609344, "ft": 0.0003048, "1000ft": 0.3048, "minutes": None}121314def pct_premium(beta, se, kennedy=True):15 """Binary X in a log-price model: 100[exp(b)-1]; Kennedy (1981) correction 100[exp(b - V/2)-1]; delta-method SE 100 exp(b) se."""16 if kennedy and se is not None and not np.isnan(se):17 est = 100 * (math.exp(beta - se ** 2 / 2) - 1); formula = "100*(exp(b - se^2/2) - 1)"18 else:19 est = 100 * (math.exp(beta) - 1); formula = "100*(exp(b) - 1)"20 se_out = 100 * math.exp(beta) * se if se is not None and not np.isnan(se) else np.nan21 return est, se_out, formula222324def semi_elasticity(beta, se, unit_scale=1.0):25 """log price, level X: 100*b per unit of X (after rescaling X to the harmonised unit: b_h = b / scale)."""26 return 100 * beta / unit_scale, 100 * se / unit_scale if se is not None else np.nan, f"100*b/{unit_scale}"272829def elasticity(beta, se):30 return beta, se, "b (log-log)"313233def partial_r(t, df):34 if t is None or df is None or np.isnan(t) or np.isnan(df) or df <= 2:35 return np.nan, np.nan, ""36 r = t / math.sqrt(t * t + df); se = math.sqrt((1 - r * r) / df)37 return r, se, "t/sqrt(t^2+df); se=sqrt((1-r^2)/df)"383940def recover_se(beta, se, t, p, ci_lo, ci_hi, df=None):41 """Recover a standard error transparently; returns (se, formula) or (nan, '')."""42 if se is not None and not np.isnan(se) and se > 0:43 return se, "reported"44 if t is not None and not np.isnan(t) and t != 0:45 return abs(beta / t), "|b/t|"46 if ci_lo is not None and ci_hi is not None and not np.isnan(ci_lo) and not np.isnan(ci_hi):47 return (ci_hi - ci_lo) / (2 * 1.959964), "(ci_hi-ci_lo)/(2*1.96)"48 if p is not None and not np.isnan(p) and 0 < p < 1 and beta != 0:49 from scipy.stats import norm50 z = norm.isf(p / 2); return abs(beta / z), "|b/z(p)| (two-sided normal)"51 return np.nan, ""525354UNIT_ALIASES = {"m2": "m2", "sqm": "m2", "sq m": "m2", "sq. m": "m2", "m²": "m2", "square meters": "m2", "square metres": "m2", "square meter": "m2", "sqft": "sqft", "sq ft": "sqft", "sq. ft": "sqft",55 "square feet": "sqft", "ft2": "sqft", "ft²": "sqft", "feet2": "sqft", "100sqft": "100sqft", "1000sqft": "1000sqft", "1,000 sqft": "1000sqft", "acres": "acre", "acre": "acre", "hectare": "ha", "ha": "ha",56 "km": "km", "kilometers": "km", "kilometres": "km", "kilometer": "km", "m": "m", "meters": "m", "metres": "m", "meter": "m", "metre": "m", "100m": "100m", "100 m": "100m", "miles": "miles", "mile": "mile",57 "ft": "ft", "feet": "ft", "foot": "ft", "1000ft": "1000ft", "1,000 ft": "1000ft", "minutes": "minutes", "min": "minutes", "years": "years", "year": "years", "yrs": "years", "count": "count", "number": "count",58 "units": "count", "index": "index", "share": "share", "percent": "percent", "%": "percent", "db": "db", "dba": "db", "db(a)": "db", "ug/m3": "ug/m3", "µg/m3": "ug/m3", "μg/m3": "ug/m3", "currency": "currency"}596061AREA_LEAVES = {"living_area", "lot_size", "garden"}62DIST_LEAVES = {"cbd_distance", "transit_distance", "school_distance", "university_distance", "hospital_distance", "highway_distance", "commercial_distance", "airport_distance", "water_proximity", "green_space", "hazard_site"}63COUNT_LEAVES = {"bedrooms", "bathrooms", "storeys", "floor_level", "garage", "parking", "fireplace", "terrace_balcony"}64AREA_UNITS = {"m2": 1.0, "sqft": 0.09290304, "1000sqft": 92.90304, "100sqft": 9.290304, "100m2": 100.0, "acre": 4046.8564224, "ha": 10000.0, "pyeong": 3.3058, "are": 100.0, "1000m2": 1000.0}65DIST_UNITS = {"km": 1.0, "m": 0.001, "100m": 0.1, "miles": 1.609344, "ft": 0.0003048, "1000ft": 0.3048, "100km": 100.0, "10km": 10.0, "100ft": 0.03048}66PLAUSIBLE = {"area_per_m2": 5.0, "dist_per_km": 60.0, "years": 10.0, "count": 60.0, "db": 5.0, "ug/m3": 10.0}676869def norm_unit(u):70 """Parse a free-text unit into a canonical token: strips parentheticals, spaces, plurals; recognises scaled units."""71 import re as _re72 u = (u or "").strip().lower(); u = _re.sub(r"\(.*?\)", "", u).strip(); u = u.replace("square", "sq").replace(",", "")73 u = _re.sub(r"\s+", " ", u)74 # scaled units like "1000 sqft", "100 m", "100 m2", "hundred years"75 m = _re.match(r"^(1000|100|10)\s*(sq ?ft|sqft|ft2|m2|m²|sqm|m|ft|km|feet)$", u)76 if m:77 n, base = m.group(1), m.group(2)78 base = {"sq ft": "sqft", "sqft": "sqft", "ft2": "sqft", "m2": "m2", "m²": "m2", "sqm": "m2", "feet": "ft", "ft": "ft", "m": "m", "km": "km"}[base]79 return f"{n}{base}"80 if u in ("hundred_years", "hundred years", "100 years", "centuries", "century"): return "hundred_years"81 if u.startswith("ping") or u.startswith("pyeong") or u.startswith("tsubo"): return "pyeong"82 if u in ("sq ft", "sqft", "sq. ft", "ft2", "ft²", "feet2", "sq feet", "sqf"): return "sqft"83 if u in ("m2", "m²", "sqm", "sq m", "sq. m", "sq meters", "sq metres", "sq meter", "square meters", "square metres", "sq.m"): return "m2"84 if u in ("acre", "acres"): return "acre"85 if u in ("hectare", "hectares", "ha"): return "ha"86 if u in ("km", "kilometer", "kilometers", "kilometre", "kilometres"): return "km"87 if u in ("m", "meter", "meters", "metre", "metres"): return "m"88 if u in ("mile", "miles"): return "miles"89 if u in ("ft", "feet", "foot"): return "ft"90 if u in ("year", "years", "yrs", "yr"): return "years"91 if u in ("count", "number", "units", "n", "integer", "rooms"): return "count"92 if u in ("db", "dba", "db(a)", "decibel", "decibels", "dB"): return "db"93 if u in ("ug/m3", "µg/m3", "μg/m3", "ug/m³", "micrograms per cubic meter"): return "ug/m3"94 if u in ("minutes", "min", "minute"): return "minutes"95 return UNIT_ALIASES.get(u, u)969798def unit_scale(tax, unit):99 """Return (scale to the harmonised unit, unit family, note) or (None, None, reason) if not harmonisable for this leaf."""100 if tax in AREA_LEAVES:101 return (AREA_UNITS[unit], "area_per_m2", f"X per {unit}->m2") if unit in AREA_UNITS else (None, None, "area unit not harmonisable")102 if tax in DIST_LEAVES:103 return (DIST_UNITS[unit], "dist_per_km", f"X per {unit}->km") if unit in DIST_UNITS else (None, None, "distance unit not harmonisable (counts, shares, indices excluded)")104 if tax == "age":105 return (1.0, "years", "per year") if unit == "years" else ((100.0, "years", "hundred years->year") if unit == "hundred_years" else (None, None, "age unit not harmonisable"))106 if tax in COUNT_LEAVES:107 return (1.0, "count", "per unit count") if unit == "count" else (None, None, "count expected")108 if tax == "noise":109 return (1.0, "db", "per dB") if unit == "db" else (None, None, "noise unit not dB")110 if tax == "air_pollution":111 return (1.0, "ug/m3", "per ug/m3") if unit == "ug/m3" else (None, None, "pollution unit not ug/m3")112 return (None, None, "no harmonised unit defined for this characteristic")113114115def harmonise_row(r):116 """r: dict with keys beta, se, t, p, ci_lo, ci_hi, dep_transform (log|level|boxcox|other), x_transform (level|log|binary|117 dummy|other), x_unit, taxonomy_id, n_obs, k_params, var_family (from taxonomy). Returns dict with family, es, es_se, formula."""118 beta = r["beta"]; se, se_formula = recover_se(beta, r.get("se"), r.get("t"), r.get("p"), r.get("ci_lo"), r.get("ci_hi"))119 dep = (r.get("dep_transform") or "").lower(); xt = (r.get("x_transform") or "").lower(); unit = norm_unit(r.get("x_unit"))120 if dep == "per_m2":121 dep = "level" # price per m2 in levels: monetary family (currency per m2 per unit of X)122 name = (r.get("var_name_printed") or "").lower()123 year_built = r.get("taxonomy_id") == "age" and any(k in name for k in ("built", "construct", "vintage", "year of", "yr built", "baujahr", "année"))124 fam = es = es_se = None; formula = ""125 if dep == "log" and xt in ("binary", "dummy"):126 es, es_se, formula = pct_premium(beta, se); fam = "pct_premium"127 elif dep == "log" and xt == "log":128 es, es_se, formula = elasticity(beta, se); fam = "elasticity"129 elif dep == "log" and xt == "level":130 sc, ufam, note = unit_scale(r.get("taxonomy_id"), unit)131 es, es_se, formula = semi_elasticity(beta, se, 1.0); fam = "semi_elasticity"132 if sc is None:133 formula = f"100*b (unit as reported: {unit or 'not stated'}; {note})"; unit_ok = False; ufam_ = None134 else:135 # b is per reported unit; per harmonised unit = b * (reported units per harmonised unit) = b / sc ... sc = harmonised units per reported unit136 es, es_se = es / sc, (es_se / sc if es_se is not None and not np.isnan(es_se) else np.nan); formula = f"100*b/{sc:.6g} ({note})"; unit_ok = True; ufam_ = ufam137 elif dep == "level" and xt in ("level", "binary", "dummy"):138 es, es_se, formula, fam = beta, se, "b (currency per unit, as reported)", "monetary"139 else:140 fam = "other"; es, es_se, formula = beta, se, "not harmonised"141 df = (r.get("n_obs") or np.nan) - (r.get("k_params") or 10)142 t = r.get("t") if r.get("t") is not None and not np.isnan(r.get("t")) else (beta / se if se and se > 0 else np.nan)143 pr, pr_se, pr_formula = partial_r(t, df)144 flags = []145 if year_built: flags.append("year_built_not_age")146 if fam == "semi_elasticity" and not locals().get("unit_ok", False): flags.append("unit_not_harmonisable")147 if fam == "semi_elasticity" and es is not None and locals().get("ufam_") and abs(es) > PLAUSIBLE.get(locals().get("ufam_"), 100): flags.append("implausible_magnitude")148 if fam == "pct_premium" and es is not None and abs(es) > 200: flags.append("implausible_magnitude")149 if fam == "elasticity" and es is not None and abs(es) > 5: flags.append("implausible_magnitude")150 if xt in ("squared", "other"): flags.append("nonstandard_x_transform")151 return {"family": fam, "es": es, "es_se": es_se, "es_formula": formula, "se_source": se_formula, "partial_r": pr, "partial_r_se": pr_se, "partial_r_formula": pr_formula, "t_used": t, "df_used": df,152 "x_unit_norm": unit, "harm_flags": ";".join(flags)}153