はじめに

　はじめまして。某大学で計量経済学を専攻している大学院生のごりわくと申します。今回は、最低限の統計の知識（線形回帰、ロジスティック回帰くらいまで）の知識を有している方を対象に、生存時間解析をpythonで実装しながら勉強していきたいと思います。初学者の方でも理解できるように丁寧に書くつもりなので、知識のある方には冗長かとは思われますがご了承ください。また、使用したコードやファイルは私のgit(https://github.com/goriwaku/survival_analysis )にも載せるので、そちらをダウンロードしてくださってもよいと思いますし、python初心者の方は自分で手を動かして入力してみるのもよいかもしれません。拙い文章となるかもしれませんが、よろしくお願いいたします。

生存時間解析とは

　生存時間解析とは、あるイベント（死亡、病気の発現など）が発生するまでにかかる時間に対して影響を与える要因を特定する医療統計学の手法です。例えば、急性心筋梗塞に伴う入院からの生存時間に対して、(1)女性は男性に比べて生存時間が長いか、(2)年齢は生存時間にどのような影響を与えるか、などを考えることが目的です。例としては医療データをあげましたが、その用途は医療分野にとどまりません。例えば、ある機械が故障するまでの時間に対する影響であったり、あるベンチャー企業が倒産するまでにかかる時間に対する影響だったりを考えることも可能で、その一番の特徴は期間を追って要因の効果を検証することができるというところです。今回は、この生存時間解析を実際にpythonで実装しながら勉強していきたいと思います。

使用するデータセット

　今回の分析では、Worcester Heart Attack Survey(WHAS)というMassachusetts州Worcesterの居住者を対象に、急性心筋梗塞を発症した被験者について収集したデータを用います。これはUCLAのINTRODUCTION TO SURVIVAL ANALYSIS IN SASから得たもので、そのままだとSAS用のファイル(sas7bdat)となってしまうため、このサイトを参考にcsvファイルに変換しました。私のgitにコードとcsvファイルは載せてあるのですが、念のため以下に手順を記します。
まず、pipコマンドで必要なモジュールをインストールします。

pip install sas7bdat

インストールが終わったら、pythonを用いて以下のコマンドを実行します。pythonをsas7bdatファイルのあるディレクトリで実行していない方は、相対パスでなく絶対パスで指定してください。

from sas7bdat import SAS7BDAT
with SAS7BDAT('whas500.sas7bdat', skip_header=False) as reader:
    df = reader.to_data_frame()
df.to_csv('whas500.csv')

ここで、df.to_csv()はpandasのよく使われるメソッドの一つで、データフレーム(df)をcsvとして出力する関数です。引数として生成するcsvファイルにつけたい名前を入れます。ここではそのままwhas500.csvとしました。実行後、この名前のcsvファイルができてているかを確認してください。

whas500は500のサンプルに対し計22個の観測変数を有しており、その一覧と意味を以下の表にまとめます。（このデータは完全なwhasデータからランダムに500サンプル抽出したものであり、完全データを対象とした解析結果と比較不可能なことに注意してください。）

先ほどcsvファイルとして保存したwhas500をpandasのDataFrameに読み込み、それにwhasという名前を付けることにします。

import pandas as pd

whas = pd.read_csv('whas500.csv').drop('Unnamed: 0', axis=1)

ファイルが読み込めたらwhas.head()などを使って先ほど挙げた変数が含まれていることを確認してみてください。なお、.drop()は不必要な変数を落とすために使用する関数で、キーワード引数axis=1のとき該当する列を、axis=0(デフォルト)のときに該当する行を削除します。次節ではこのデータを用いて生存時間データについて解説していきます。

生存時間データ

　この節では、生存時間解析に使用されるデータについて解説していきます。生存時間データに特有のデータとして打ち切りと呼ばれるものがあります。これは、最後の追跡日まで関心のあるイベントが発生しなかった標本のことを表します。例えば、先ほどの急性心筋梗塞に伴う入院からの生存時間の例において、観察を終了とした時点でもなお生存している被験者がいた場合、その被験者の生存時間は観察期間を超過したことはわかりますが、実際にいつイベントが発現するかは不明です（観察終了翌日かもしれませんし、20年ほど生き続けるかもしれません）。したがって、この標本の観測値は不完全なものとなります。このような観測値のことを右側打ち切りと呼びます。ここでは、pythonのmatplotlibと呼ばれる描画ライブラリを用いて、この打ち切りについてみていきます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

whas = pd.read_csv('whas500.csv').drop('Unnamed: 0', axis=1)
for i in range(10):
    x = (0, whas.LENFOL[i])
    y = (whas.ID[i], whas.ID[i])
    plt.plot(x, y, color='black')
    if whas.FSTAT[i] == 1:
        plt.plot(whas.LENFOL[i], whas.ID[i], 'x', color='red')
    else:
        plt.plot(whas.LENFOL[i], whas.ID[i], 'o', color='blue')
plt.show()

実務上、生存時間は観察開始日と観察最終日の差を数えることによって得ることができます。よって、被験者の観察がいつ始まっていつ終わるのかを日数だけで判断し、観察開始日をt=0と定義します。上のコードによって描画されるグラフが以下になります。

横軸が生存時間、縦軸が被験者IDです。右端の青の点と赤のバツ印は被験者が生存か死亡かを表しています。ここに見られる青の点が右側打ち切りの例です。生存時間解析における不完全な観測値の原因は打ち切りと切断の2つがあります。打ち切りとは、例えば被験者が引っ越すため試験が継続できなくなるように、興味のあるイベントが観測される前に観測が打ち切られてしまう状況です。それに対し切断とは、実験のデザイン上、観察終了日を定めて実験を行う際に、観察終了日まで興味のあるイベントが発生しなかったサンプルのことを指します。この時、観察開始もランダムであるため、打ち切りも切断も同じ右側打ち切りとして扱うことは妥当であると考えられます。また、この記事では深くは触れませんが、打ち切りには左側打ち切りという開始時点以前に興味のあるイベントが発現しているケース、および連続的に観察できず、2か月などの一定の時間ごとに観察結果を得られる区間打ち切りと呼ばれる打ち切りが存在します。

今回は生存時間データについてpythonで描画することによってみることができました。次回は、生存時間解析における単変量解析の手法、Kaplan-Meier法についてみていきたいと思います。

次回の記事はこちら
Pythonで学ぶ生存時間解析2 -Kaplan-Meier推定量
https://qiita.com/Goriwaku/items/767cdc2640e29fddf7cc

参考文献・参照リンク

生存時間解析入門　Hosmer DW, Lemeshow S, May S
Introduction to Survival Analysis in SAS https://stats.idre.ucla.edu/sas/seminars/sas-survival/
sas7bdat 2.2.3 https://pypi.org/project/sas7bdat/
京都大学OCW　京都大学大学院医学研究科 聴講コース 臨床研究者のための生物統計学「生存時間解析の基礎」　https://www.youtube.com/watch?v=NmZaY2tDKSA&feature=emb_title