Add statistical models: chi-square tests and fatal-severity logistic regression

- Chi-square independence tests (rain / time period / first-party vehicle
group vs fatal outcome) with Cramer's V effect sizes
- Accident-level logistic regression on fatality: heavy vehicles OR 6.7,
early morning OR 2.7, male OR 2.4, +10km/h speed limit OR 1.3 (all p<0.05);
rain not significant for severity
- Outputs to 分析結果/模型結果/ with rare-event caveat in model summary

src/taipei_traffic/models.py

@@ -0,0 +1,154 @@
+"""統計檢定與事故嚴重度模型。
+
+1. 卡方獨立性檢定：雨天／時段／第一當事人車種大類 × 是否為死亡事故
+2. 邏輯迴歸：以事故層級預測「是否為死亡事故」（24小時內或2-30日內有人死亡）
+
+死亡事故僅 72 件（佔 0.32%），屬罕見事件：迴歸結果著重方向與相對風險
+（勝算比），不宜直接作為機率預測器使用。
+"""
+
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import statsmodels.api as sm
+from scipy import stats as scipy_stats
+
+from .config import DB_PATH, MODEL_DIR
+from .db import connect
+
+_MODEL_QUERY = """
+SELECT a.is_fatal, a.is_rain, a.time_period, a.weekday, a.speed_limit,
+       p.age, p.gender_code, p.vehicle_category
+FROM accidents a
+JOIN parties p ON p.accident_id = a.accident_id AND p.party_seq = 1
+"""
+
+# 死亡事故稀少，車種合併為大類以避免完全分離
+_VEHICLE_GROUPS = {
+    "機車": "機車",
+    "小客車": "小客車",
+    "小貨車": "小貨車",
+    "大客車": "大型車",
+    "大貨車": "大型車",
+    "聯結車": "大型車",
+    "曳引車": "大型車",
+    "慢車": "慢車",
+    "人（行人/乘客）": "其他",
+    "其他車": "其他",
+    "特種車": "其他",
+    "軍車": "其他",
+}
+
+_VEHICLE_REF = "小客車"
+_PERIOD_REF = "下午 (12:00-17:59)"
+
+
+def load_model_data(db_path: Path = DB_PATH) -> pd.DataFrame:
+    with connect(db_path) as conn:
+        df = pd.read_sql_query(_MODEL_QUERY, conn)
+    df["vehicle_group"] = df["vehicle_category"].map(_VEHICLE_GROUPS)
+    df["is_weekend"] = (df["weekday"] >= 5).astype(int)
+    return df
+
+
+def _cramers_v(table: np.ndarray, chi2: float) -> float:
+    n = table.sum()
+    k = min(table.shape) - 1
+    return float(np.sqrt(chi2 / (n * k))) if k > 0 else float("nan")
+
+
+def chi_square_tests(df: pd.DataFrame) -> pd.DataFrame:
+    """對是否死亡事故做三組卡方獨立性檢定。"""
+    tests = {
+        "雨天 × 是否死亡事故": df["is_rain"],
+        "時段 × 是否死亡事故": df["time_period"],
+        "第一當事人車種大類 × 是否死亡事故": df["vehicle_group"],
+    }
+    rows = []
+    for name, factor in tests.items():
+        table = pd.crosstab(factor, df["is_fatal"]).to_numpy()
+        chi2, p, dof, _ = scipy_stats.chi2_contingency(table)
+        rows.append(
+            {
+                "檢定": name,
+                "卡方值": round(chi2, 3),
+                "自由度": dof,
+                "p值": round(p, 5),
+                "Cramér's V": round(_cramers_v(table, chi2), 4),
+                "顯著(α=0.05)": "是" if p < 0.05 else "否",
+            }
+        )
+    return pd.DataFrame(rows)
+
+
+def _design_matrix(df: pd.DataFrame) -> tuple[pd.DataFrame, pd.Series]:
+    d = df[
+        df["gender_code"].isin([1, 2])
+        & df["age"].notna()
+        & df["speed_limit"].notna()
+        & df["vehicle_group"].notna()
+    ].copy()
+
+    X = pd.DataFrame(index=d.index)
+    X["雨天"] = d["is_rain"]
+    X["週末"] = d["is_weekend"]
+    X["年齡(每+10歲)"] = d["age"] / 10
+    X["男性"] = (d["gender_code"] == 1).astype(int)
+    X["速限(每+10km/h)"] = d["speed_limit"] / 10
+
+    period_dummies = pd.get_dummies(d["time_period"], prefix="時段").astype(int)
+    period_dummies = period_dummies.drop(columns=f"時段_{_PERIOD_REF}")
+    X = pd.concat([X, period_dummies], axis=1)
+
+    vehicle_dummies = pd.get_dummies(d["vehicle_group"], prefix="車種").astype(int)
+    vehicle_dummies = vehicle_dummies.drop(columns=f"車種_{_VEHICLE_REF}")
+    X = pd.concat([X, vehicle_dummies], axis=1)
+
+    X = sm.add_constant(X)
+    return X, d["is_fatal"]
+
+
+def severity_logit(df: pd.DataFrame):
+    """配適死亡事故邏輯迴歸，回傳 (fit結果, 勝算比表)。"""
+    X, y = _design_matrix(df)
+    fit = sm.Logit(y, X.astype(float)).fit(disp=False)
+
+    conf = fit.conf_int()
+    odds = pd.DataFrame(
+        {
+            "變項": fit.params.index,
+            "係數": fit.params.round(4).values,
+            "勝算比(OR)": np.exp(fit.params).round(3).values,
+            "OR 95%CI下界": np.exp(conf[0]).round(3).values,
+            "OR 95%CI上界": np.exp(conf[1]).round(3).values,
+            "p值": fit.pvalues.round(5).values,
+        }
+    )
+    odds = odds[odds["變項"] != "const"].reset_index(drop=True)
+    return fit, odds
+
+
+def run_all(db_path: Path = DB_PATH, out_dir: Path = MODEL_DIR) -> None:
+    out_dir.mkdir(parents=True, exist_ok=True)
+    df = load_model_data(db_path)
+
+    chi = chi_square_tests(df)
+    chi.to_csv(out_dir / "卡方檢定結果.csv", index=False, encoding="utf-8-sig")
+    print("✓ 卡方檢定結果.csv")
+    print(chi.to_string(index=False))
+
+    fit, odds = severity_logit(df)
+    odds.to_csv(out_dir / "邏輯迴歸_勝算比.csv", index=False, encoding="utf-8-sig")
+    print("\n✓ 邏輯迴歸_勝算比.csv")
+    print(odds.to_string(index=False))
+
+    summary_path = out_dir / "邏��迴歸_模型摘要.txt"
+    n_events = int(df["is_fatal"].sum())
+    caveat = (
+        "註：死亡事故為罕見事件（{} 件 / {} 件，約 {:.2%}），本模型用於解讀風險\n"
+        "因子方向與相對強度（勝算比），非機率預測器。參考組別：時段=下午、\n"
+        "車種=小客車。\n\n"
+    ).format(n_events, len(df), n_events / len(df))
+    summary_path.write_text(caveat + str(fit.summary()), encoding="utf-8")
+    print(f"✓ {summary_path.name}（McFadden pseudo R² = {fit.prsquared:.3f}）")

分析結果/模型結果/卡方檢定結果.csv

@@ -0,0 +1,4 @@
+檢定,卡方值,自由度,p值,Cramér's V,顯著(α=0.05)
+雨天 × 是否死亡事故,0.615,1,0.43284,0.0052,否
+時段 × 是否死亡事故,9.194,3,0.02682,0.0203,是
+第一當事人車種大類 × 是否死亡事故,71.658,5,0.0,0.0567,是

分析結果/模型結果/邏輯迴歸_勝算比.csv

@@ -0,0 +1,14 @@
+變項,係數,勝算比(OR),OR 95%CI下界,OR 95%CI上界,p值
+雨天,-0.327,0.721,0.344,1.513,0.38723
+週末,0.0099,1.01,0.575,1.775,0.97268
+年齡(每+10歲),0.1913,1.211,1.047,1.401,0.01006
+男性,0.8547,2.351,1.096,5.041,0.02812
+速限(每+10km/h),0.2888,1.335,1.01,1.763,0.0421
+時段_上午 (06:00-11:59),0.3189,1.376,0.794,2.384,0.25571
+時段_夜晚 (18:00-23:59),-0.1838,0.832,0.409,1.694,0.61225
+時段_清晨 (00:00-05:59),0.9892,2.689,1.13,6.396,0.02527
+車種_其他,1.8253,6.205,2.067,18.627,0.00114
+車種_大型車,1.8983,6.674,3.192,13.957,0.0
+車種_小貨車,1.1825,3.263,1.434,7.42,0.00479
+車種_慢車,0.3315,1.393,0.321,6.049,0.65812
+車種_機車,0.1141,1.121,0.622,2.02,0.70397

分析結果/模型結果/邏輯迴歸_模型摘要.txt

@@ -0,0 +1,31 @@
+註：死亡事故為罕見事件（72 件 / 22368 件，約 0.32%），本模型用於解讀風險
+因子方向與相對強度（勝算比），非機率預測器。參考組別：時段=下午、
+車種=小客車。
+
+                           Logit Regression Results                           
+==============================================================================
+Dep. Variable:               is_fatal   No. Observations:                22071
+Model:                          Logit   Df Residuals:                    22057
+Method:                           MLE   Df Model:                           13
+Date:                Sat, 11 Jul 2026   Pseudo R-squ.:                 0.06586
+Time:                        20:44:03   Log-Likelihood:                -452.22
+converged:                       True   LL-Null:                       -484.11
+Covariance Type:            nonrobust   LLR p-value:                 1.103e-08
+=======================================================================================
+                          coef    std err          z      P>|z|      [0.025      0.975]
+---------------------------------------------------------------------------------------
+const                  -9.0942      0.902    -10.085      0.000     -10.862      -7.327
+雨天                     -0.3270      0.378     -0.865      0.387      -1.068       0.414
+週末                      0.0099      0.288      0.034      0.973      -0.554       0.574
+年齡(每+10歲)               0.1913      0.074      2.574      0.010       0.046       0.337
+男性                      0.8547      0.389      2.196      0.028       0.092       1.618
+速限(每+10km/h)            0.2888      0.142      2.033      0.042       0.010       0.567
+時段_上午 (06:00-11:59)     0.3189      0.281      1.137      0.256      -0.231       0.869
+時段_夜晚 (18:00-23:59)    -0.1838      0.363     -0.507      0.612      -0.895       0.527
+時段_清晨 (00:00-05:59)     0.9892      0.442      2.237      0.025       0.123       1.856
+車種_其他                   1.8253      0.561      3.254      0.001       0.726       2.925
+車種_大型車                  1.8983      0.376      5.044      0.000       1.161       2.636
+車種_小貨車                  1.1825      0.419      2.821      0.005       0.361       2.004
+車種_慢車                   0.3315      0.749      0.443      0.658      -1.137       1.800
+車種_機車                   0.1141      0.300      0.380      0.704      -0.475       0.703
+=======================================================================================