[Stata][生存時間解析]StataでのCox比例ハザードモデル

graphs display on an ipad stata
Photo by Burak The Weekender on Pexels.com

stsetとstxcox

今回は以下の2つのPDFをざっと読んで見る。

https://www.stata.com/manuals/ststset.pdf#ststset

https://www.stata.com/manuals/ststcox.pdf

stset

まずはstsetから。stとはsurvival-timeのことで、Stataでは変数が生存時間データであることを宣言する役割を持つ。stsetはメモリ上のデータを生存時間データ(stデータ)として宣言できる。例えば、特定のイベント(例:死亡)が起こった時刻を表す変数tvarを使い、

stset tvar

という形でコマンドを指定する。上記のコマンドで、一人あたりのレコードについての生存時間データを設定する。stsetしたデータは保存することができ、その後は再度 stset する必要がない。Stataは設定を記憶する。

Quick start

無修正データでの失敗時刻の指定: 失敗時刻を記録したtvar変数を使い、生データに対して

stset tvar

と指定する。

イベントの指定: event == 2 が失敗を表す場合、

stset tvar, failure(event == 2)

として、特定のイベントを失敗として指定する。また、以下のように複数のイベント(例:event == 2event == 3)を失敗として指定することも可能である。

stset tvar, failure(event == 2 3)

また、”death (0/1)”のように、変数で指定することも可。

stset tvar, failure(death)

いつから被験者がリスクに晒されているかを以下のように指定する。

stset tvar, failure(fail) origin(time torig)

以下のコードは、被験者が時刻0にリスクとなるが、時刻tenterに研究に参加することを指定する。

stset tvar, fail(failure) enter(time tenter)

以下のその組み合わせで、被験者がtorig時点でリスクになり、tenter時点で試験に参加することを指定する。

stset tvar, fail(failure) origin(time torig) enter(time tenter)

上記と同じだが、時間変数が日単位で測定される場合、分析時間を年単位で指定する。

stset tvar, fail(failure) origin(time torig) enter(time tenter) scale(365.25)

分析時間の単位を日数に戻すのは以下のコード。

streset, scale(1)

以前のst設定を表示し、データの変更が設定に対応していることを確認する。

stset

被験者ごとの複数記録の生存データがある場合、以下のようにIDを指定する。

stset tvar, failure(fail) id(idvar)

他は同様。

stcox

stsetは難しくない。stsetを使って、生存時間データに変換したものを用いて、Cox比例ハザードモデルを解析する。こちらもQuick startから。

Quick start

stsetデータを用いた共変量x1とx2によるCox比例ハザードモデル

stcox x1 x2

上記と同じだが、エフロン方式を使用する。

stcox x1 x2, efron

svarのレベルによって定義された層におけるベースラインハザードの違いを見るには以下。

stcox x1 x2, strata(svar)

svysetおよびstsetデータを使用して複雑な調査デザインを調整するときは以下。

svy: stcox x1 x2

構文

以下にoptionsと説明を記載する。下線の部分だけ入力すれば、省略してコマンドを走らせることができる。

Model説明
\(\underline{\text{est}}\text{imate}\)共変量なしの適合モデル
strata(変数名)層別化
\(\underline{\text{sh}}\text{ared(varname)}\)共有脆弱性ID変数
\(\underline{\text{off}}\text{set(varname)}\)varnameを係数1に制約されたモデルに含める
\(\underline{\text{bre}}\text{slow}\)ブレスローメソッドを使って失敗を処理する(これがdefault)。
\(\underline{\text{efr}}\text{on}\)エフロンメソッドで失敗を処理する
exactm厳密な限界尤度法で失敗を扱う
exactp厳密な偏尤度法で失敗を扱う
Time varying
tvc(varlist)時間の関数と相互作用させる共変量を指定する。
texp(exp)時間の関数を指定; デフォルトは texp(_t)
SE/Robust
vce(vcetype)vcetype は \(\text{oim}\), \(\underline{\text{r}}\text{obust}\), \(\underline{\text{cl}}\text{uster}\) \(\text{clustvar}\), \(\underline{\text{boot}}\text{strap}\), \(\underline{\text{jack}}\text{knife}\) のいずれか。
noadj(ust)標準的な自由度調整を使用しない
Reporting
\(\underline{\text{l}}\text{evel(#)}\)信頼水準を設定する; デフォルトは level(95)
nohrハザード比ではなく係数を報告する
\(\underline{\text{nosh}}\text{ow}\)ST設定情報を表示しない
stcoxのオプション

生存曲線の描写

このあとPDFではいくつか具体的なデータを使って例が表示されている。Cox比例ハザードモデル自体はそれほど難しくないので、それは飛ばして図の描写に移る。

stset tvar, failure(event)
* tvarはeventまでの時間、eventはevent == 2 3などでも指定可能

stcox exposure covar1 covar2
* exposureは主要アウトカム、covar1, covar2は共変量

stcurve, survival at1(exposure = 1) at2(exposure = 2) scheme(s2mono) xtitle(time) ytitle(Survival rate) title("Survival Curve") xlabel(#) legend(order(1 "Intervention" 2 "Control"))
* リスクテーブルなし

これでグラフが描写できる。リスクテーブルを下に記載したいときは、以下のようにする。

stset tvar, failure(event)
* tvarはeventまでの時間、eventはevent == 2 3などでも指定可能

sts graph, by(exposure) risktable(, order 1 "Intervention" 2 "Control") scheme(s2mono) xtitle(time) ytitle(Survival rate) title("Survival Curve") xlabel(#) legend(order(1 "Intervention" 2 "Control"))
* リスクテーブルあり

また、Log-rank検定を行い、その結果を図に含めたいときは、以下のようにすれば良い。

stset tvar, failure(event)
* tvarはeventまでの時間、eventはevent == 2 3などでも指定可能

sts test event
* Log-rank検定

sts graph, by(exposure) risktable(, order 1 "Intervention" 2 "Control") scheme(s2mono) xtitle(time) ytitle(Survival rate) title("Survival Curve") xlabel(#) legend(order(1 "Intervention" 2 "Control")) ttext(x y "Log-rank test" " " "p = ...")
* リスクテーブルあり、Log-rank結果あり

\(x, y\)でグラフ中の座標を指定する。

下2つの解析は、共変量の調整ができないことに注意する。

関連記事

Stataでのdummy変数の生成
時間の扱い
marginsplotの調整方法
Stataでのマルチレベル分析について
そもそもマルチレベル分析とはなにか

コメント

タイトルとURLをコピーしました