Stata分析機能例題集
生存分析

データの確認

最初に、基本となるグラフを作成します。
下記のコマンドで例題用のサンプルデータ「drug2」を入手します。

サンプルデータには、以下のデータが入力されています。

Sutdytime :死亡または打ち切りまでの月数
died ： 1=死亡、0=打ち切り
drug ： 1=服薬あり、0=服薬なし
age ： 観察開始時の患者の年齢

投薬の有無（drug）の違いによる生存率の時間変化をグラフにします。

このデータは、すべての患者のデータが観察開始時間（t=0）から始まっています。
t=0以降に観察を始める患者のデータを追加すること（遅延組入）はしていません。
そこで、次はt=0以降に新規患者のデータを追加していく場合のグラフを作成します。

遅延組入のあるサンプルデータ「drug2b」を入手します。

1つ目のサンプルデータと違って、t=5やt=8から観察を開始した患者のデータがあります。
先程と同様にグラフを作成します。

lostオプションは打ち切り人数を表示するオプションなので、遅延組入の人数は負の値として表示されます。
この「-1」は、例えば「遅延組入1人」や「遅延組入2人、打ち切り1人」のような状態を表します。

enterオプションを使用すると、遅延組入と打ち切りを分けて表示できます。

遅延組入と打ち切りを分けて表示する方法は、常に最善とは限りません。
次のサンプルデータ「drug2c」は、最初の「drug2」と同じ内容のデータです。
しかし、記述方法が異なり、時間経過によって変化する共変量（投薬量や体重など）を併記できるように、1人の患者に対して複数行のデータがあります。
遅延組入はありません。

ID ： 患者ID

サンプルデータ「drug2c」を入手し、lostオプションを付けてグラフを作成します。

このグラフは、「drug2」のグラフと同じです。データの内容も、形式が違うだけで理論的には同じ事象を表しています。
しかし、enterオプションを使用すると、たくさんの数値が記入されてしまいます。

同じ患者IDのデータが複数行あり、一行ごとに遅延組入・打ち切りとして処理されるためです。
遅延組入と打ち切りが区別できない場合はlostとenterのオプションを使うと良かったのですが、このような共変量のあるデータの場合は適しません。

データの確認

cumhazオプションを使用すると、Nelson-Alen累積ハザード関数を描くことができます。
最初と同じサンプル「drug2」を使用して、投薬の有無による累積ハザード関数をグラフにします。

打ち切りの人数を表示するグラフは、lostオプションで作成できます。

stsコマンドは、ハザード関数の推定にも使用することができます。
下記のグラフは、weightedカーネル平滑化を基にハザード関数を推定しています。
カーネル関数とバンド幅の選び方によって結果が変わるため、注意が必要です。

次のコマンドでは、カーネル平滑化とバンド幅を調整しています。
titileは、グラフタイトルを変更するコマンドです。

リスクテーブルを追加する

以下のコマンドで、リスクテーブルをグラフの下に追加します。

書式オプションは下記の通りです。
凡例の位置を変更します。

リスクテーブルのラベル列を「Placebo」と「Test drug」に変更します。

ラベル列の文字を左揃えにします。

リスクテーブルのタイトルをとラベル列の行頭を揃えます。

データの確認

stcoxコマンドを使用して、比例Coxハザードモデルを適用します。
打ち切りのないデータでCox回帰を行います。
新型のベアリングを搭載した非常用発電機の耐久実験を行います。
この実験では、保護回路を無効にして発電機が発火するまで負荷をかけて運転を行います。

サンプルデータを入手し、データを確認します。

failtime ： 故障するまでの時間（h）
lost ： 負荷（kVA）
bearing ： 1=新型のベアリング、0=旧型のベアリング

12台の発電機があり、内訳は新型ベアリング搭載機と旧型ベアリング搭載機が各6台ずつです。
例えば、1行目のデータは旧式ベアリング搭載機に15 lVAの負荷をかけて運転したところ、100時間後に故障したことを意味します。
Cox比例ハザードモデルに当てはめて、発電機の故障は負荷の程度に影響を受けるのか、それともベアリングの型によるのかを調べます。
ベアリングの型と負荷の大きさは、ハザード関数の形全体には影響しないものとします。

下記のコマンドでモデルを適用します。

負荷の程度の影響を制御すると、新型ベアリングはハザード率が6.36％と低く正常作動時間も長くなることが分かります。
一度stcoxを適用すると、以降は引数を省略できます。

続けて、上記のモデルにハザード比ではなく係数を表示オプションnohrを使用して再表示させます。

打ち切りのあるデータでCox回帰を行う

下記のような48人の癌患者の服薬データがあります。
21人は薬物療法を行っていて、20人は行っていません。
患者の年齢は47才から67才までです。
このデータを使って、患者が死亡するまでの月数について解析を行います。

sutdytime ： 死亡または打ち切りまでの月数
died ： 1=死亡、0=打ち切り
drug ： 1=服薬あり、0=服薬なし
age ： 観察開始時の患者の年齢

データを入手してstsetでデータを整え、サマリーを確認します。

各変数のデータ数、平均、標準偏差、最小値、最大値が表示されます。
Coxモデルを適用します。

年齢を制御すると、投薬ありの場合のハザード率が10.5％と低く、生存時間が長いことが分かります。
ハザード率は、変数が1単位変わる時の変化を表しています。
上記の例では、年齢が1才上がるごとにハザードが12％増加しています。
5才年齢が変わる場合のリスクを計算します。
年齢を5才で1単位とするために、新しい変数を作成します。

再度Coxモデルを適用します。

年齢が5才上がると、ハザードは76％増加します。

タイ（同じ値）のあるデータを扱う

ハザードモデルを構築するには、タイのない連続的な値を扱うのが理想的です。
タイが発生した場合には、部分尤度を計算する必要があります。
タイのあるデータに対してCox部分尤度を計算するために、Stataには4つのオプション（brslow, efron, exactm, exactp）が用意されています。

時間で変化する共変量のあるデータでCox回帰を行う

スタンフォード心臓移植プログラムの心臓移植データを例にします。
これには、実際に心臓移植を受けた患者と受けられなかった患者（合計103人）のデータが含まれています。
心臓移植を待っている間に死亡した患者や待機をやめた患者もいますが、67％の患者は移植を受けました。

サンプルデータを入手し、データを確認します。

id ： 患者ID
year ： 心臓移植プログラムに参加した年
age ： 患者の年齢
died ： 1=死亡、0=打ち切り
stime ： 生存時間（日）
surgey ： 1=外科手術あり、0=外科手術なし（CABGなど）
transplant ： 1=心臓移植実施、0=心臓移植なし
wait ： 移植までの待機時間
posttran ： post-transplant 1=新しい心臓を受け取った、0=受け取らなかった

stsetでデータを整え、Coxモデルを適用します。

高齢の患者ほどハザード率が高く、手術を受けた患者のハザード率が低いことが分かります。
また、最終的に心臓移植を受けたかどうかによって、ハザード率は大きく変わりません。

継続的に時間変化する共変量のあるデータでCox回帰を行う

基本的なCox回帰は、次の式で表されます。

h(t)=h_0 (t)exp(β_1 x_1+⋯+β_k+x_k )

変数がz_i (t)=z_i g(t)に変化したとすると、式は下記のようになります。

h(t)=h_0 (t)exp{├ β_1 x_1+⋯+β_k+x_k ┤+g(t)┤ ├ (γ_1 z_1+⋯+γ_m z_m )}

z_1,⋯⋯,z_mが時間変化する共変量し推定に影響を及ぼす場合、回帰係数γ_1と共変量g(t) z_iは現在の関数です。

変数z_1,⋯⋯,z_mを決定するにはtvc(varlist)オプションを使用します。
g(t)におけるtは分析時間なので、g(t)はtexp(exp)オプションで計算できます。
例えば、g(t)=log(t)を計算する場合は下記のように入力します。_tはtsetによって計算されています。

Cox回帰はイベント発生時の部分尤度を基に計算されているので、イベント発生時にstsplitを使って分割し、手動で時間変化のある共変量を生成することでも計算できます。
stsplitを使ってイベント発生の多い大きなデータセットを処理すると大量のメモリを消費するため、その場合はtvc()やtexp()を使用します。