HARRY CHANG
存活分析：Kaplan-Meier 與 Cox 回歸 - 職災風險窗口與設備失效時間
—

存活分析：Kaplan-Meier 與 Cox 回歸 - 職災風險窗口與設備失效時間

一般模型問「會不會出事」，存活分析問「多久內會出事」。本篇介紹醫學統計的鎮山之寶——Kaplan-Meier 曲線與 Cox 比例風險模型，用 600 名新進作業員的首次工傷資料，找出職災風險窗口、量化速成訓練的代價，並揭露一個天真統計會踩死的陷阱：刪失資料。

WRITTEN BY

  • avatar
    Name
    Harry Chang
    Twitter

下一篇

電腦視覺：Vision Transformer (ViT) - 廠房裂縫影像分類，CNN 王座保衛戰
← 回到列表

核心貢獻者

Edward Kaplan 與 Paul Meier 於 1958 年各自向期刊投稿了同一個方法,編輯要求兩人合寫——這篇合著成為統計學史上引用數最高的論文之一 (超過 6 萬次),Kaplan-Meier 曲線從此成為醫學論文的標配。Sir David Cox (大衛·考克斯爵士) 則在 1972 年提出比例風險模型 (Proportional Hazards Model),天才之處在於「半母數」設計:不必知道基準風險長什麼樣子,照樣能估計每個因子的影響倍數。這兩件武器統治生物醫學統計半世紀,但在工廠裡,它們同樣鋒利。

為什麼工安需要存活分析？

前面幾篇的模型都在回答「會不會」:這班次會不會出事 (Day 41)、這操作是不是異常 (Day 40)。存活分析換一個更貼近管理決策的問法:「多久之內」——

  1. 職災風險窗口:新進員工進廠後,前幾個月風險最高?加強帶班該帶多久?
  2. 設備失效時間:這批設備平均撐多久?保養週期訂 3 個月還是 6 個月?
  3. 措施效果量化:速成訓練比完整訓練「危險幾倍」?——不是差幾個百分點,是幾倍。

而它存在的根本理由,是一種普通回歸與分類都處理不了的資料型態——刪失 (Censoring)。

1. 資料集來源

資料集：合成新進作業員「到首次工傷的時間」

備註：延續本系列做法。同一套方法把「員工」換成「設備」、「工傷」換成「故障」,一個字的程式碼都不用改。

資料集特色與欄位介紹:

600 名新進作業員,追蹤 24 個月,記錄「到職後多久發生第一次工傷」:

  • duration (月):觀察到的時間長度。
  • event:1 = 發生工傷;0 = 刪失——追蹤結束時還平安。
  • fast_track (速成訓練):0 = 完整安全訓練,1 = 速成上線 (佔 40%)。
  • high_risk_dept (高風險部門):沖壓/成型 vs 包裝。
  • overtime (加班程度)、age (年齡):皆標準化;年齡被我故意設計成毫無影響,看模型會不會誠實說出來。

刪失:這份資料的靈魂

600 人中只有 366 人 (61%) 真的出事,其餘 234 人 (39%) 是刪失——原因有二:追蹤滿 24 個月還平安,或中途離職失聯。

刪失資料示意
  • 紅色 X:事件發生 (工傷)。
  • 藍色圈:刪失——我們只知道「他至少平安了這麼久」,不知道之後會怎樣。

天真統計的陷阱:如果只拿 366 個出事者算平均,得到「平均 5.86 個月出事」——嚴重高估風險!因為撐得久的人大多還沒出事就被觀察期切掉了,根本沒進到你的平均數裡。刪失資料不能丟掉 (丟掉就是只看倒楣鬼),也不能當成事件 (他們明明還平安)。存活分析就是為了正確地「榨出刪失資料裡的資訊」而生。

2. 原理

2.0 兩個核心函數

  • 存活函數 S(t)S(t):撐過時間 t 還平安的機率。從 1 開始一路往下掉。
  • 風險函數 h(t)h(t):已經撐到 t 的前提下,下一瞬間出事的瞬時風險。「風險窗口」找的就是 h(t) 的高峰。

2.1 Kaplan-Meier:一步一腳印的條件機率

KM 估計的邏輯樸素到優美——把「活過 t」拆成一連串條件機率的連乘:

S^(t)=∏ti≤t(1−dini)\hat{S}(t) = \prod_{t_i \le t} \left(1 - \frac{d_i}{n_i}\right)

  • tit_i:每個「有事件發生」的時間點;did_i:當時出事的人數;nin_i:當時還在風險中的人數。
  • 刪失的人怎麼處理:他們待在分母 nin_i 裡直到刪失那一刻,然後安靜退場——貢獻了「他平安的那段時間」的資訊,又不污染事件計數。這就是 KM 的全部魔法。
  • Log-rank 檢定:比較兩條 KM 曲線是否真的不同 (虛無假設:兩組風險相同)。

2.2 Cox 比例風險模型:量化每個風險因子的倍數

KM 只能一次比一個因子,Cox 把所有因子放進同一個模型:

h(t∣x)=h0(t)⋅exp⁡(β1x1+β2x2+… )h(t \mid x) = h_0(t) \cdot \exp(\beta_1 x_1 + \beta_2 x_2 + \dots)

  • h0(t)h_0(t) 基準風險:所有人共享的時間形狀 (新人期高、老手期低)。Cox 的天才在於根本不去估計它——半母數設計,只估係數 β\beta。
  • 風險比 (Hazard Ratio, HR) =eβ= e^{\beta},整個模型最好懂的輸出:
    • HR = 2.0:這個因子讓風險翻倍。
    • HR = 1.0:沒有影響。
    • HR = 0.5:保護因子,風險砍半。
  • 比例風險假設:HR 不隨時間改變 (速成訓練的危害在第 1 個月和第 20 個月都是同一個倍數)。用之前要檢查,violated 時可用分層或時變係數處理。

3. 實戰

Python 程式碼實作

from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

# ---- Kaplan-Meier:按訓練組別各畫一條 ----
kmf = KaplanMeierFitter()
for val, label in [(0, 'full training'), (1, 'fast-track training')]:
    m = df.fast_track == val
    kmf.fit(df.duration[m], df.event[m], label=label)   # durations + 事件旗標
    kmf.plot_survival_function(show_censors=True)
    print(label, kmf.median_survival_time_)             # 中位存活時間

# ---- Log-rank 檢定:兩條曲線真的不同嗎 ----
res = logrank_test(df.duration[m0], df.duration[m1], df.event[m0], df.event[m1])

# ---- Cox 回歸:全部因子一起上 ----
cph = CoxPHFitter()
cph.fit(df, duration_col='duration', event_col='event')
cph.print_summary()          # exp(coef) 就是 Hazard Ratio

程式碼重點:

  • lifelines 是 Python 存活分析的標準套件,API 與 sklearn 神似。
  • 餵給模型的永遠是兩欄一組:duration + event——刪失資訊就藏在 event=0 裡。

4. 模型評估

先對答案:天真平均 vs KM 中位數

統計方式結果評語
只算出事者的平均5.86 個月錯! 刪失者被丟掉,風險被誇大
KM 中位存活 (完整訓練組)12.36 個月正確處理刪失後的答案
KM 中位存活 (速成訓練組)5.19 個月速成組的中位數只有完整組的 42%

Kaplan-Meier 曲線:風險窗口現形

KM 曲線比較
  • 兩組天差地遠:log-rank p = 1.2e-12,速成訓練組 (紅) 的曲線一路被完整訓練組 (藍) 壓在下面。中位數 5.19 vs 12.36 個月——速成訓練省下的時數,用員工的安全月數付了帳。
  • 職災風險窗口:注意兩條曲線都在前 3~6 個月掉得最陡,之後趨緩——這就是風險函數遞減的形狀,「新人前 90 天」正是加強帶班、增加巡檢的黃金窗口。曲線上的小豎線是刪失標記。

Cox 回歸:每個因子的風險倍數

風險因子HR (exp(coef))95% 信賴區間p 值解讀
速成訓練2.211.80 ~ 2.73顯著風險翻 2.2 倍,本資料最大元凶
高風險部門2.011.62 ~ 2.48顯著沖壓/成型的宿命,靠工程控制解
加班程度1.301.16 ~ 1.45顯著每多一個標準差,風險 +30%
年齡1.070.96 ~ 1.190.199不顯著——模型誠實還原了設計
Cox 風險比森林圖
  • 模型的 Concordance Index 為 0.650 (隨機猜是 0.5):對「誰先出事」的排序有中等鑑別力——個體命運仍有大量隨機性,但群體層級的風險倍數估得很準。
  • 年齡那一行是本篇彩蛋:資料生成時年齡係數設為 0,Cox 給出 HR 1.07、信賴區間跨過 1、p=0.199——模型沒有把雜訊說成因果。一個誠實說「我不知道」的模型,比什麼都顯著的模型可信。

同一套數學,直接搬去修設備

把「員工」換成「射出機」、「工傷」換成「故障停機」、「速成訓練」換成「延遲保養」,整套分析原封不動:KM 中位數就是中位失效時間,Cox 的 HR 告訴你「保養每延遲一個月,故障風險乘幾倍」,而 KM 曲線掉到某個存活率的時間點,就是預防保養週期的科學訂法。可靠度工程裡的 Weibull 分析,正是存活分析的工業孿生兄弟。

5. 總結

我們學習了存活分析——一般 ML 教材幾乎不教、但工安與設備管理天生需要的武器:

  • 刪失資料:不能丟、不能亂補,「他至少平安了這麼久」本身就是資訊。只算出事者的平均 (5.86 個月) 是會誤導決策的錯誤統計。
  • Kaplan-Meier:條件機率連乘的存活曲線;曲線最陡的區段就是風險窗口——新人前 3~6 個月。
  • Cox 回歸:不猜基準風險、只估風險倍數;速成訓練 HR 2.21,而且模型誠實地放過了無辜的年齡。
  • 工安啟示:存活分析把「要不要加強訓練」從理念之爭變成算術題——速成訓練讓風險翻 2.2 倍、中位平安時間砍半,這筆帳拿去跟管理層報告,比任何口號都有力。

下一篇回到電腦視覺,迎接新世代的挑戰者:Vision Transformer (ViT)——把影像切成 16×16 的拼圖丟進 Transformer,CNN 的王座還坐得穩嗎?

View on GitHub