接續上一篇内容，如果只有 8 年的城市指標數據的情況下想做趨勢分析，線性迴歸對小樣本和離群值都不友善。這時候我們可以使用 Mann–Kendall，可以套用在這樣的應用場景 —— 非參數、不假設常態、配上 Sen's slope 估斜率比 OLS 還抗噪。

在此我們先討論為什麼線性迴歸在這裡不夠用

首先，線性迴歸研究城市 8 年發展，最直覺的做法是評判「有沒有上升趨勢」

from scipy import stats
slope, intercept, r, p, se = stats.linregress(years, values)

p < 0.05 = 趨勢顯著。是否就能下結論呢，但是會遇到三個問題：

問題 1：常態性假設
線性迴歸的 p 值依賴「殘差常態分佈」。在 8 個點的基礎下沒辦法檢驗這個假設。

問題 2：對離群值敏感
建築數量數據裡，2019 年因為 OSM 社群活躍度低，數值反常掉到 12,771（其他年都在 16,000+）。OLS 估出來的斜率會被這個點拉低。

問題 3：假設線性
城市發展常常是「先慢後快」（前期積累、後期爆發）。而線性擬合會嚴重低估後期的成長速率。

Mann–Kendall 解決前兩個問題（非參數、用秩相關所以抗離群），Sen's slope 解決第三個（用所有點對的中位數斜率，不假設線性）。

二、Mann–Kendall 的直覺

Mann–Kendall 檢定（Mann–Kendall Test，簡稱 MK 檢定）是一種用於分析時間序列資料的非參數統計方法。它的主要目的是判斷一組資料隨著時間演進，是否存在單調趨勢（也就是持續上升、持續下降，或是沒有明顯趨勢）

Mann–Kendall 在城市分析中提出問題：

  「任取一對年份，後面那年的值比前面那年大的次數多嗎？」

計算公式：

S = Σᵢ<ⱼ sign(xⱼ - xᵢ)

對所有 i < j 的點對，後面比前面大記 +1、小記 -1、相等記 0。最後加總。

  • S > 0 → 大部分點對是「後面比前面大」 → 上升趨勢
  • S < 0 → 下降趨勢
  • S ≈ 0 → 沒有趨勢

關鍵：只看「大小關係」，不看具體數值。所以離群值不會影響 S（離群值再大、跟其他點的「大小關係」也就是 ±1）。

程式碼如下：

S = 0
for i in range(n - 1):
    for j in range(i + 1, n):
        diff = x[j] - x[i]
        if diff > 0:
            S += 1
        elif diff < 0:
            S -= 1

O(n²) 但 n=8 不會慢。對於 n > 100 可以用 numpy 矩陣化加速，但在本案例中沒必要。

三、S 的分佈與 Z 檢定

在虛無假設 H₀（值是隨機排列、無趨勢）下，S 的期望是 0，方差有封閉解：

Var[S] = [n(n-1)(2n+5) - Σ tᵢ(tᵢ-1)(2tᵢ+5)] / 18

其中 tᵢ 是「第 i 個重複值出現的次數」（連結值校正）。

對應的 Z 分數：

       (S - 1) / √Var[S]     若 S > 0
Z =    0                     若 S = 0
       (S + 1) / √Var[S]     若 S < 0

那個 ±1 叫連續性校正，因為 S 是離散值但 Z 用連續分佈近似。

程式碼

unique_vals, tie_counts = np.unique(x, return_counts=True)
tie_term = sum(t * (t - 1) * (2 * t + 5) for t in tie_counts if t > 1)
Var_S = (n * (n - 1) * (2 * n + 5) - tie_term) / 18.0

if Var_S <= 0:
    Z_MK = 0.0
elif S > 0:
    Z_MK = (S - 1) / np.sqrt(Var_S)
elif S < 0:
    Z_MK = (S + 1) / np.sqrt(Var_S)
else:
    Z_MK = 0.0

p_value = float(2 * (1 - stats.norm.cdf(abs(Z_MK))))

雙尾 p < 0.05 即為顯著趨勢。

四、Sen's slope：抗噪的斜率估計

使用 Mann–Kendall 可以獲取「是否有趨勢」、「趨勢是正還是負」，但無法得知「每年漲多少」。所以要用 Sen's slope 補上這塊

Sen's slope（中文常稱為 Sen 斜率估計 或 泰爾－森估算，Theil–Sen estimator）是一種用來計算時間序列資料「變化趨勢斜率」的非參數統計方法。

β_Sen = median { (xⱼ - xᵢ) / (j - i) : i < j }

對所有 i < j 點對，算出 (n choose 2) 個斜率，取中位數。

OLS（Ordinary Least Squares）是線性回歸最常用的方法，但在面對現實世界不完美的數據時，OLS 有很大的局限性。

跟 OLS 作比較：

  • OLS 是用「最小化殘差平方和」找斜率，離群值貢獻平方誤差，影響很大
  • Sen 是「所有點對斜率的中位數」，離群點對只貢獻幾個極端斜率，被中位數過濾掉
  • 結果：對污染穩健、不要求殘差常態

程式碼

slopes = []
for i in range(n - 1):
    for j in range(i + 1, n):
        if j != i:
            slopes.append((x[j] - x[i]) / (j - i))
Sen_slope = float(np.median(slopes)) if slopes else np.nan

五、執行套用城市數據

把面板數據裡每個指標的「年度均值時序」餵給 Mann–Kendall：

def run_mann_kendall(panel, indicators=None):
    if indicators is None:
        exclude = {"year", "lon", "lat"}
        indicators = [c for c in panel.select_dtypes(include="number").columns
                      if c not in exclude and not c.startswith("d_")]

    years = sorted(panel["year"].unique())
    results = []

    for ind in indicators:
        yearly_mean = panel.groupby("year")[ind].mean().reindex(years).values
        valid = ~np.isnan(yearly_mean)
        if valid.sum() < 4:
            continue

        mk = mann_kendall(yearly_mean[valid])
        results.append({
            "indicator": ind,
            "n_years": int(valid.sum()),
            "series_mean": float(np.nanmean(yearly_mean)),
            "series_std": float(np.nanstd(yearly_mean)),
            **mk,
        })

    return pd.DataFrame(results).sort_values("p_value")

執行結果如下：

預期輸出（按 p 值排序）：

指標 S Z_MK p 值 趨勢 Sen 斜率
building_count +28 4.04 0.0001 *** increasing +14.5
poi_diversity +27 3.89 0.0001 *** increasing +0.023
shop_count +25 3.60 0.0003 *** increasing +0.155
amenity_count +24 3.45 0.0006 *** increasing +0.27
texture_entropy +20 2.86 0.0042 ** increasing +0.010
edge_density +18 2.57 0.0103 * increasing +0.0042
building_coverage +14 1.97 0.0488 * increasing +0.0021
building_area_mean -6 -0.81 0.4186 no trend -3.5

讀法：

  • 建築數量 S = 28，幾乎所有點對都是後面比前面大 → 強上升趨勢 (p < 0.001)
  • POI 多樣性、商店、設施、紋理熵、邊緣密度、建築覆蓋率全都顯著上升
  • 建築平均面積 S = -6，無顯著趨勢 → 「每棟建築變小」不顯著，但「建築數量大幅增加」顯著
  • 這呼應前面結論：台中是「數量驅動型開發」

Sen 斜率的解讀：

  • 建築數量 +14.5/年 → 每點每年平均多 14.5 棟建築
  • POI 多樣性 +0.023/年 → 8 年累計 +0.18，跟 0.44→0.60 (+0.16) 吻合

六、Mann–Kendall vs 線性迴歸：兩個跑一遍對比

拿 building_count 的 8 年均值來比：

2018: 27.6
2019: 31.0
2020: 34.8
2021: 47.6
2022: 67.8
2023: 109.8
2024: 131.3
2025: 136.0

OLS 線性迴歸：斜率 = +17.4/年，R² = 0.86，p < 0.001
Sen's slope：+14.5/年

兩個結論一致都「強顯著上升」，但 Sen 斜率比 OLS 小 17%。原因：OLS 被後期 2023-2025 三年的爆發拉高了；Sen 取中位數，早期慢、後期快的「彎曲」對它沒有影響。

而這正是 Sen 抗噪的價值：當你不確定增長是線性還是加速時，Sen 給的是「中位增速」，比 OLS 的「平均增速」更保守、更可解釋。

七、在這裡我們也説明下什麼時候不該用 Mann–Kendall

雖然 MK 很強，但也有限制：

1. n < 4 直接無法計算

if n < 4:
    return {"S": np.nan, "Var_S": np.nan, "Z_MK": np.nan,
            "p_value": np.nan, "trend": "insufficient data", "Sen_slope": np.nan}

如果是 4 點以下，S 的取值範圍太窄（-6 到 +6），統計檢定沒意義。

2. 有強季節性的數據

如果你的時序有季節週期（譬如月度資料），直接 MK 會被季節性主導。要先 deseasonalize，或用 Seasonal Mann–Kendall（每個季節獨立算 S 再加總）。

3. 有自相關的時序

如果今年的值跟去年高度相關（譬如「今年建築數 = 去年 + 新增」），MK 會高估顯著性。這時要用 Modified Mann–Kendall（Hamed-Rao 校正）。

對城市發展這種「年度快照」資料，自相關沒那麼嚴重 —— OSM 每年快照彼此獨立。可以直接用標準 MK。

4. 想抓「轉折點」

MK 只告訴你「整體有沒有趨勢」，不告訴你「趨勢從哪一年改變」。這要用 Chow 結構突變點檢定。

九、總結

規則總結：

  • 大樣本（n > 30）、無離群、確定線性 → 使用 OLS 沒問題
  • 小樣本、有離群、不確定線性 → 則使用 MK + Sen
  • 同時跑兩個進行分析則結論一致時更有說服力

無標籤

關注作者：

新增評論