Stanで生存時間解析(Weibull 回帰)
生存時間解析とは?
生存時間解析は、イベントの時間を解析するための手法です。例えば、
- ソーシャルゲームやwebサービスなどに登録した人の利用継続時間(マーケティング)
- 投薬群と対照群(プラセボ)で、どれだけ長生きするか(医用統計)
- 新規材料が既存の材料とくらべて、どれぐらい耐久性があるか(信用工学)
など、様々な応用例が考えられます。生存時間解析はノンパラメトリックな手法で行なうことが多いです。具体的には、
- カプランマイヤー推定量で生存時間を推定
- 生存時間の違いをログランク検定で判断
- 共変量の違いをCoxの比例ハザードモデルやAelanの加法モデルで解析
といった辺りが挙げられるでしょうか。その原因としては、生存時間のデータへ打ち切りという特徴があることが挙げられます。結果として分布にあてはめるだけの解析では、打ち切りを説明することが難しく上手くフィッティング出来ません。並べて、最尤法を使った分布への当てはめでは、値の信頼性が出ません。これに対してノンパラメトリックなカプランマイヤー法を使えば、推定値の範囲(信頼性)も出せますので論文も書きやすいです。
しかし仮に分布を仮定しても、ベイズ的にモデリングすれば、これらの問題を解決することができます。今回はStanを使って生存時間解析を行なうことで、打ち切りを含めてモデリングをしてみました。何分不勉強なもので、ご指導ご鞭撻をいただければ幸いです。
テストデータ
とりあえず解析に使うデータを用意します。今回はpythonのlifeline packageに含まれる、rossiのデータを使います。データのソース論文はこちらです。このデータは逮捕された人が再逮捕されるまでの期間(週)と、その他の共変量が記されています。何はともあれ、データを見てみましょう。上から10件だけ表示します。
from lifelines.datasets import load_rossi df = load_rossi() df
index | week | arrest | fin | age | race | wexp | mar | paro | prio |
---|---|---|---|---|---|---|---|---|---|
0 | 20 | 1 | 0 | 27 | 1 | 0 | 0 | 1 | 3 |
1 | 17 | 1 | 0 | 18 | 1 | 0 | 0 | 1 | 8 |
2 | 25 | 1 | 0 | 19 | 0 | 1 | 0 | 1 | 13 |
3 | 52 | 0 | 1 | 23 | 1 | 1 | 1 | 1 | 1 |
4 | 52 | 0 | 0 | 19 | 0 | 1 | 0 | 1 | 3 |
5 | 52 | 0 | 0 | 24 | 1 | 1 | 0 | 0 | 2 |
6 | 23 | 1 | 0 | 25 | 1 | 1 | 1 | 1 | 0 |
7 | 52 | 0 | 1 | 21 | 1 | 1 | 0 | 1 | 4 |
8 | 52 | 0 | 0 | 22 | 1 | 0 | 0 | 0 | 6 |
9 | 52 | 0 | 0 | 20 | 1 | 1 | 0 | 0 | 0 |
10 | 52 | 0 | 1 | 26 | 1 | 0 | 0 | 1 | 3 |
詳細は論文の3.2節に書いてありますが、今回注目するのはweekとarrestです。arrestが1である人は再逮捕されています。0である人は観測期間中では再逮捕はされず、weekの値が52になっています。所謂打ち切りの発生です。
モデル
アヒル本P115 7.8 打ち切り、に従い、現象をモデル化してみます。なお今回のデータはweek単位への丸めも発生しているのですが、本題ではないので飛ばします。
打ち切りの原理を踏まえて、今回の現象のモデルを次のように考えました。
「再逮捕の区間を表すパラメータと分布に従って、再逮捕までの時間が得られる。ただしである場合は、とする 」
分布はワイブル分布を利用します。ワイブル分布は次のような特徴をもち、時間のモデル化へ向いています。
- における累積密度関数の和が1になる(イベントの発生確率を無限時間まで考えられる)
- 現象の傾向を示す形状パラメータと、現象の時間スケールを示す尺度パラメータで表現される
- なら、時間とともにハザード関数値(時間までにイベントが発生しなかった場合、その瞬間にイベントが起こる確率)が低くなる
- なら、時間に依らずハザード関数値が一定
- なら、時間にとともにハザード関数値が高くなる
ワイブル分布を使って、値に対する尤度を考えます。観測値が52で無い場合は
として得られます。肝心なのは52の場合です。この時尤度は
となります。
以上を踏まえて、まず共変量に依らず、全ての人が共通の現象に従っているとして、Stanを書いてみました。
data { int N ; int week[N] ; int arrest[N] ; } parameters { real shape ; real scale ; } model { for(n in 1:N){ if(arrest[n] == 0){ target += weibull_lccdf(week[n]| shape, scale) ; }else{ target += weibull_lpdf(week[n]| shape, scale) ; } } } generated quantities { real log_lik[N] ; for(n in 1:N){ if(arrest[n] == 0){ log_lik[n] = weibull_lccdf(week[n]| shape, scale) ; }else{ log_lik[n] = weibull_lpdf(week[n]| shape, scale) ; } } }
解析と可視化
データに対して、構築したStanモデルを当てはめてみます。
import numpy as np import os import pystan import pandas as pd import pickle from lifelines.datasets import load_rossi def stan_compile(stan_file_path, compiled_file_path, recompile=False): if os.path.exists(compiled_file_path) is False or recompile is True: model = pystan.StanModel(file=stan_file_path) with open(compiled_file_path, "wb") as f: pickle.dump(model, f) else: with open(compiled_file_path, "rb") as f: model = pickle.load(f) return model def get_summary_df(fit): summary = fit.summary(probs=prob) summary_df = pd.DataFrame(summary["summary"], index=summary["summary_rownames"], columns=summary["summary_colnames"]) return summary_df def get_waic(fit): log_lik = fit.extract("log_lik")["log_lik"] waic = (-2 * np.sum(np.log(np.mean(np.exp(log_lik), axis=0))) + 2 * np.sum(np.var(log_lik, axis=0))) return waic # Stanに入れるデータの準備 df = load_rossi() N = len(df) week = df["week"] arrest = df["arrest"] stan_data = {"N": N, "week": week, "arrest": arrest} prob = [0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975] # Stanをコンパイル stan_file = os.path.join("weibull.stan") stan_file_c = os.path.join("weibull.pkl") model = stan_compile(stan_file, stan_file_c) # StanへKick fit = model.sampling(data=stan_data, iter=3000, chains=5, warmup=1000, seed=1234) # 結果の集計 summary_df = get_summary_df(fit) waic = get_waic(fit) print(waic) summary_df
output
WAIC: 1397.22843377
parameter | mean | se_mean | sd | 2.5% | 5% | 25% | 50% | 75% | 95% | 97.5% | n_eff | Rhat |
---|---|---|---|---|---|---|---|---|---|---|---|---|
shape | 1.35177 | 0.00225134 | 0.123002 | 1.12027 | 1.15716 | 1.26817 | 1.347 | 1.43093 | 1.56233 | 1.61164 | 2985 | 1.00366 |
scale | 127.559 | 0.267877 | 14.5396 | 103.759 | 106.833 | 117.305 | 125.893 | 136 | 153.639 | 161.079 | 2946 | 1.00274 |
log_lik[0] | -5.26822 | 0.000938137 | 0.0870498 | -5.44863 | -5.41606 | -5.3248 | -5.26577 | -5.20816 | -5.13211 | -5.10419 | 8610 | 0.999675 |
log_lik[1] | -5.30895 | 0.000968553 | 0.0933687 | -5.5044 | -5.46746 | -5.36939 | -5.30494 | -5.2457 | -5.16148 | -5.13464 | 9293 | 0.999847 |
log_lik[2] | -5.21905 | 0.00110178 | 0.0852228 | -5.39228 | -5.36399 | -5.27516 | -5.21646 | -5.16062 | -5.08346 | -5.05936 | 5983 | 1.00011 |
… | … | … | … | … | … | … | … | … | … | … | … | … |
まず集計表です。全パラメータでが確認できましたので、収束は問題ないようです。続いて、得られたパラメータ値を用いて可視化を行います。
from scipy.stats import weibull_min import matplotlib.pyplot as plt # 可視化その1 fig, ax1 = plt.subplots() x = np.arange(0, 500, 1) y = weibull_min.pdf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "mean"]) y_min = weibull_min.pdf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "5%"]) y_max = weibull_min.pdf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "95%"]) # 確率密度関数 ax1.plot(x, y, "-") ax1.fill_between(x, y_min, y_max, facecolor="pink") ax1.set_ylabel("probability") ax1.tick_params("y") ax1.set_xlabel("Weeks") # 件数 ax2 = ax1.twinx() ax2.hist(df["week"]) ax2.set_ylabel("count") ax2.tick_params("y") plt.show() # 可視化その2 fig, ax1 = plt.subplots() x = np.arange(0, 50, 1) y = weibull_min.sf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "mean"]) y_min = weibull_min.sf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "5%"]) y_max = weibull_min.sf(x, c=summary_df.loc["shape", "mean"], scale=summary_df.loc["scale", "95%"]) # 生存時間関数 ax1.plot(x, y, "-") ax1.set_ylabel("Propotion Not Arrested") ax1.fill_between(x, y_min, y_max, facecolor="pink") ax1.set_xlabel("Weeks") ax1.tick_params("y")
一枚目のグラフが、確率密度関数とイベント発生のカウントです。MCMCのEAPが青線、scale parameterの90%信用区間を用いて描かれた範囲がピンク色になっています。
二枚目のグラフが生存時間関数(survival function)です。見方は1枚目と同様です。元論文のP7には、Coxの比例ハザードモデルで描かれた場合の結果が載っています。今回のWeibull分布を用いて描いた結果と相違ありません。
打ち切り補正をしない場合は?
当然ここで気になってくるのは、補正の効果です。先程の1枚目のグラフは52辺りに最頻値があります。打ち切りをしている効果はどの程度でしょうか?これを検証するために、打ち切りを考慮しないモデルと比較しました。
打ち切りを考慮しない場合は、めちゃめちゃ過剰適合しているのがわかります。打ち切りのメカニズムをモデルに取り入れることで、52以上でもイベントが生起することを踏まえてモデル化できてるわけですね。ちなみに打ち切りが考慮されない場合、WAICも3411.15613467となり、モデルとしても性能が良くないことが分かります。
比例ハザードモデル
次に共変量の効果を調べます。一番単純な比例ハザードモデルで取り組んでみましょう。比例ハザードモデルは、全時間を通してハザード間数値を一定として、共変量に従ってハザード関数の値が大きくなることを仮定するモデルです。Weibull分布のハザード関数は
と示されます。先行研究[1][2]では、尺度パラメータへ共変量の効果を加えています。この理由の僕個人による定性的な理解としては、形状パラメータの頑強性に帰するものと考えています。形状パラメータの値はそもそもの現象に依存します。共変量ぐらいでは、中々現象の傾向までは大きく変わりません。
更に、尺度に対する効果も尺度パラメータへダイレクトに加えるのではなく、へ加えるほうが多いようです。即ち基底状態に対するパラメータを、共変量の係数を、共変量のダミー変数をとして、
ということになります。この式をを中心に変形すれば
です。Stanで実装するときは、この式で書けばよいかと思います。
今回の解析では、逮捕後の経済的支援(financial aid)の効果を調べてみたいと思います。モデルは次の通りです。
data { int N ; int week[N] ; int arrest[N] ; int fin[N] ; } parameters { real shape ; real beta[2] ; } model { for(n in 1:N){ if(arrest[n] == 0){ target += weibull_lccdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ; }else{ target += weibull_lpdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ; } } } generated quantities { real log_lik[N] ; for(n in 1:N){ if(arrest[n] == 0){ log_lik[n] = weibull_lccdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ; }else{ log_lik[n] = weibull_lpdf(week[n]| shape, exp(- (beta[1] + fin[n] * beta[2]) / shape)) ; } } }
スクリプトは殆ど同様なので飛ばします。結果として次のグラフが得られました。
2枚目のグラフが分かりやすいですね。経済的支援を受けた場合は、再犯率が低くなっています。WAICも1395.3474841となって低くなり、予測力が上がっています。更に言うと、経済的支援の効果を表すパラメータに対する、の95%信用区間(片側)はでした。この区間に1(即ち効果が無い)が入っていませんので、経済的支援に効果があったことが期待されます。
まとめ
Python + Stan + Weibull分布を使って、打ち切りを含めたモデル化を行い、共変量の効果の信用区間まで計算しました。…あれ、これはCoxの比例ハザードモデルで十分なのでは(死)。
ベイズ的にモデル化することの利点は自由性です。パッと思いつくものとして、次のような応用が挙げられます。
- 加法モデルやオッズ比モデルの利用
- 値の丸めを加味したモデル
- 時間tに対して、ハザード関数が変化するモデル(比例ハザードモデルでは一定とする)
- リンク関数の中身へ正規化項を追加して、スパースにする(Lasso + weibull 回帰?)
- 単純な2群比較のコホートデザインではない場合(利用できる共変量が他にある)
ふぅー、これでやった意味が見いだせました。ぜひ皆様も、Weibull分布に親しんで下さい。
参考
書籍
- 作者:豊田 秀樹
- 出版社/メーカー: 朝倉書店
- 発売日: 2017/01/25
- メディア: 単行本(ソフトカバー)
2節「ワイブル分布」がドンピシャでした。
打ち切りデータのモデリングについて勉強させて貰いました。
Webリソース
- Modeling Language User's Guide and Reference Manual, Version 2.16.0
- Cox Proportional-Hazards Regression for Survival Data in R
- COMPARISON BETWEEN WEIBULL AND COX PROPORTIONAL HAZARDS MODELS
- Estimating trends in data from the Weibull and a generalized extreme value distribution
- Lifelines
- 比例ハザードモデルはとってもtricky!
- Bayesian Lassoで特徴選択