Easy to type

非正規労働者の職業訓練記録です。ボーナスと福利厚生を勝ち取る夢を持っています。

Stanで生存時間解析(Weibull 回帰)

生存時間解析とは?

生存時間解析は、イベントの時間を解析するための手法です。例えば、

など、様々な応用例が考えられます。生存時間解析はノンパラメトリックな手法で行なうことが多いです。具体的には、

  • カプランマイヤー推定量で生存時間を推定
  • 生存時間の違いをログランク検定で判断
  • 共変量の違いをCoxの比例ハザードモデルやAelanの加法モデルで解析

といった辺りが挙げられるでしょうか。その原因としては、生存時間のデータへ打ち切りという特徴があることが挙げられます。結果として分布にあてはめるだけの解析では、打ち切りを説明することが難しく上手くフィッティング出来ません。並べて、最尤法を使った分布への当てはめでは、値の信頼性が出ません。これに対してノンパラメトリックなカプランマイヤー法を使えば、推定値の範囲(信頼性)も出せますので論文も書きやすいです。

しかし仮に分布を仮定しても、ベイズ的にモデリングすれば、これらの問題を解決することができます。今回はStanを使って生存時間解析を行なうことで、打ち切りを含めてモデリングをしてみました。何分不勉強なもので、ご指導ご鞭撻をいただければ幸いです。

テストデータ

とりあえず解析に使うデータを用意します。今回はpythonlifeline 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単位への丸めも発生しているのですが、本題ではないので飛ばします。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

打ち切りの原理を踏まえて、今回の現象のモデルを次のように考えました。

「再逮捕の区間を表すパラメータと分布に従って、再逮捕までの時間 yが得られる。ただし 52 \lt yである場合は、 y = 52とする 」

分布はワイブル分布を利用します。ワイブル分布は次のような特徴をもち、時間のモデル化へ向いています。

  •  0 \lt t \lt \inftyにおける累積密度関数の和が1になる(イベントの発生確率を無限時間まで考えられる)
  • 現象の傾向を示す形状パラメータ mと、現象の時間スケールを示す尺度パラメータ \etaで表現される
    •  m \lt 1なら、時間とともにハザード関数値(時間 tまでにイベントが発生しなかった場合、その瞬間にイベントが起こる確率)が低くなる
    •  m = 1なら、時間に依らずハザード関数値が一定
    •  m \gt 1なら、時間にとともにハザード関数値が高くなる

ワイブル分布を使って、値に対する尤度を考えます。観測値が52で無い場合は

weibull\_lpdf(y|m, \eta)

として得られます。肝心なのは52の場合です。この時尤度は


\begin{align}
Prob[52 \lt y]   &= \int_{52}^{\infty}weibull(52|m, \eta) \\\
    &= 1 - \int_0^{52} weibull(52|m, \eta) \\\
    &= 1 - weibull\_cdf(52|m, \eta) \\\
    &= weibull\_ccdf(52|m, \eta)\end{align}

となります。

以上を踏まえて、まず共変量に依らず、全ての人が共通の現象に従っているとして、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

まず集計表です。全パラメータで \hat{R} \lt 1.1が確認できましたので、収束は問題ないようです。続いて、得られたパラメータ値を用いて可視化を行います。

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")

f:id:ajhjhaf:20170805173622p:plain

f:id:ajhjhaf:20170805173634p:plain

一枚目のグラフが、確率密度関数とイベント発生のカウントです。MCMCEAPが青線、scale parameterの90%信用区間を用いて描かれた範囲がピンク色になっています。

二枚目のグラフが生存時間関数(survival function)です。見方は1枚目と同様です。元論文のP7には、Coxの比例ハザードモデルで描かれた場合の結果が載っています。今回のWeibull分布を用いて描いた結果と相違ありません。

打ち切り補正をしない場合は?

当然ここで気になってくるのは、補正の効果です。先程の1枚目のグラフは52辺りに最頻値があります。打ち切りをしている効果はどの程度でしょうか?これを検証するために、打ち切りを考慮しないモデルと比較しました。

f:id:ajhjhaf:20170805175603p:plain

打ち切りを考慮しない場合は、めちゃめちゃ過剰適合しているのがわかります。打ち切りのメカニズムをモデルに取り入れることで、52以上でもイベントが生起することを踏まえてモデル化できてるわけですね。ちなみに打ち切りが考慮されない場合、WAICも3411.15613467となり、モデルとしても性能が良くないことが分かります。

比例ハザードモデル

次に共変量の効果を調べます。一番単純な比例ハザードモデルで取り組んでみましょう。比例ハザードモデルは、全時間を通してハザード間数値を一定として、共変量に従ってハザード関数の値が大きくなることを仮定するモデルです。Weibull分布のハザード関数は


\begin{align}
    h(t|m, \eta) &= \frac{m}{\eta}(\frac{x}{\eta})^{m-1} \\\
    &= \frac{m}{\eta^m}x^{m-1} \\\
    &= m \lambda x^{m-1} \\\
\end{align}

と示されます。先行研究[1][2]では、尺度パラメータ \etaへ共変量の効果を加えています。この理由の僕個人による定性的な理解としては、形状パラメータの頑強性に帰するものと考えています。形状パラメータの値はそもそもの現象に依存します。共変量ぐらいでは、中々現象の傾向までは大きく変わりません。

更に、尺度に対する効果も尺度パラメータ \etaへダイレクトに加えるのではなく、 \lambdaへ加えるほうが多いようです。即ち基底状態に対するパラメータを \beta_0、共変量の係数を \beta、共変量のダミー変数をFとして、


    \lambda = exp(\beta_0 + \beta^TF )

ということになります。この式を \etaを中心に変形すれば


    \eta = exp(-\frac{\beta_0 + \beta^TF}{m})

です。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)) ;
        }
    }
}

スクリプトは殆ど同様なので飛ばします。結果として次のグラフが得られました。

f:id:ajhjhaf:20170805185919p:plain

f:id:ajhjhaf:20170805185927p:plain

2枚目のグラフが分かりやすいですね。経済的支援を受けた場合は、再犯率が低くなっています。WAICも1395.3474841となって低くなり、予測力が上がっています。更に言うと、経済的支援の効果を表すパラメータ \beta_1に対する、exp(\beta_1)の95%信用区間(片側)は 0.502 \lt exp(\beta_1) \lt 0.948でした。この区間に1(即ち効果が無い)が入っていませんので、経済的支援に効果があったことが期待されます。

まとめ

Python + Stan + Weibull分布を使って、打ち切りを含めたモデル化を行い、共変量の効果の信用区間まで計算しました。…あれ、これはCoxの比例ハザードモデルで十分なのでは(死)。

ベイズ的にモデル化することの利点は自由性です。パッと思いつくものとして、次のような応用が挙げられます。

  • 加法モデルやオッズ比モデルの利用
  • 値の丸めを加味したモデル
  • 時間tに対して、ハザード関数が変化するモデル(比例ハザードモデルでは一定とする)
  • リンク関数の中身へ正規化項を追加して、スパースにする(Lasso + weibull 回帰?)
  • 単純な2群比較のコホートデザインではない場合(利用できる共変量が他にある)

ふぅー、これでやった意味が見いだせました。ぜひ皆様も、Weibull分布に親しんで下さい。

参考

書籍

実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-

2節「ワイブル分布」がドンピシャでした。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

打ち切りデータのモデリングについて勉強させて貰いました。

Webリソース