Stata: Cox比例ハザードモデル
Stataを用いて、Kaplan-Meier曲線を作成し、多変量Cox比例ハザードモデルを作成してみましょう。
今回は、サンプルデータCancer Drug Trialを利用します。
. use https://www.stata-press.com/data/r16/drugtr, clear
(Patient Survival in Drug Trial)Cancer Drug Trialデータセット (drugtr) には、下記の変数が含まれています。
- studytime: 追跡期間（月）
- died: アウトカム（0=生存、1=死亡）
- drug: 介入因子（1=薬剤介入群、0=プラセボ）
- age: 年齢
Kaplan-Meier曲線を作成する
生存解析を始める前に、stsetコマンドで追跡期間とアウトカムを設定します。
. stset studytime, failure(died==1)
failure event: died == 1
obs. time interval: (0, studytime]
exit on or before: failure
------------------------------------------------------------------------------
48 total observations
0 exclusions
------------------------------------------------------------------------------
48 observations remaining, representing
31 failures in single-record/single-failure data
744 total analysis time at risk and under observation
at risk from t = 0
earliest observed entry t = 0
last observed exit t = 39対象48人、アウトカム発症31人であることがわかります。
まずは、薬剤介入群とプラセボ群のKaplan-Meier曲線を作成してみましょう。sts graphコマンドを使います。図のサイズは、A4（8.27 × 11.69インチ）に2 x 3個のグラフが収まるサイズ（3.7 × 3.7 インチ）に設定します。
. sts graph, by(drug) xsize(3.7) ysize(3.7)
failure _d: died
analysis time _t: studytime
薬剤介入群の累積生存率は、プラセボ群よりもかなり良いみたいです。二群の累積生存率をLog-rank検定で比較してみましょう。sts testコマンドを利用します。
. sts test drug
failure _d: died
analysis time _t: studytime
Log-rank test for equality of survivor functions
| Events Events
drug | observed expected
------+-------------------------
0 | 19 7.25
1 | 12 23.75
------+-------------------------
Total | 31 31.00
chi2(1) = 28.27
Pr>chi2 = 0.0000Log-rank検定のP値 (Pr>chi2) < 0.001ですので、薬剤介入群とプラセボ群の累積生存率は、統計学的に有意に異なると結論づけることができます。
Cox比例ハザードモデルを作成する
それでは、Cox比例ハザードモデルを用いて、プラセボ群に対する薬剤介入群のハザード比を算出してみましょう。stcoxコマンドを利用します。
. stcox drug
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -88.254734
Iteration 2: log likelihood = -88.001551
Iteration 3: log likelihood = -88.00019
Refining estimates:
Iteration 0: log likelihood = -88.00019
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(1) = 23.82
Log likelihood = -88.00019 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
drug | .1327581 .0584002 -4.59 0.000 .0560555 .3144157
------------------------------------------------------------------------------プラセボ群（drug = 0）に対する薬剤介入群（drug = 1）のハザード比は、0.1327581 [0.0560555 – 0.3144157]であり、P = 0.000です。stcoxでは、変数の最も小さい値が自動的にreferenceに設定されます。任意の値をreferenceに設定したければ、変数の前にb.任意の値を設定します。薬剤介入群（drug = 1）をreferenceに設定したければ、b1.drugです。
. stcox b1.drug
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -88.254734
Iteration 2: log likelihood = -88.001551
Iteration 3: log likelihood = -88.00019
Refining estimates:
Iteration 0: log likelihood = -88.00019
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(1) = 23.82
Log likelihood = -88.00019 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
0.drug | 7.532496 3.313542 4.59 0.000 3.180502 17.83948
------------------------------------------------------------------------------続いて、本モデルが比例ハザード性を満たしているかを確認します。比例ハザード性を確認する方法はいろいろありますが、今回はSchoenfeld残差を利用した方法で確認してみましょう。estat phtestコマンドを利用して、直前に作成したCox比例ハザードモデルの比例ハザード性を検定します。
. estat phtest, detail
Test of proportional-hazards assumption
Time: Time
----------------------------------------------------------------
| rho chi2 df Prob>chi2
------------+---------------------------------------------------
0.drug | 0.00948 0.00 1 0.9583
1b.drug | . . 1 .
------------+---------------------------------------------------
global test | 0.00 1 0.9583
----------------------------------------------------------------P = 0.9582ですので、比例ハザード性が成立していないわけではないと解釈できます。
次は、年齢で補正した多変量Cox比例ハザードモデルを作成してみましょう。stcoxに年齢を追加します。
. stcox drug age
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -83.551879
Iteration 2: log likelihood = -83.324009
Iteration 3: log likelihood = -83.323546
Refining estimates:
Iteration 0: log likelihood = -83.323546
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(2) = 33.18
Log likelihood = -83.323546 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
drug | .1048772 .0477017 -4.96 0.000 .0430057 .2557622
age | 1.120325 .0417711 3.05 0.002 1.041375 1.20526
------------------------------------------------------------------------------年齢＋1歳あたりのハザード比は1.120325 [1.041375 – 1.20526]であり、P = 0.002です。比例ハザード性も確認してみましょう。
. estat phtest, detail
Test of proportional-hazards assumption
Time: Time
----------------------------------------------------------------
| rho chi2 df Prob>chi2
------------+---------------------------------------------------
drug | 0.00949 0.00 1 0.9603
age | -0.11758 0.42 1 0.5168
------------+---------------------------------------------------
global test | 0.43 2 0.8064
----------------------------------------------------------------いずれの変数もP > 0.05ですので、比例ハザードが成立していないわけではないということなります。
もし、P < 0.05となってしまって、比例ハザード性が成立していない可能性がでてきたときは、graphicalに評価してみましょう。大規模コホート研究では、検出力が高すぎて、しばしばP < 0.05になってしまいます。estat phtestコマンドに、plot(変数名)オプションを追加します。
estat phtest, plot(age) xsize(3.7) ysize(3.7)
lowessカーブが、時間に依存せずに横一直線になっていれば、比例ハザード性は成立していると考えます。対象数が少ないので、なんとも言えない感じですが、これくらいならばフラットと考えていいでしょう、多分。
このモデルでは、年齢1歳あたりハザード比が1.12 [1.04-1.20]だけ上昇することがわかりました。しかしながら、年齢1歳あたりでは幅が小さすぎて、臨床的感覚にマッチしません。そこで、年齢10歳あたりハザード比がどれくらい上がるかを算出してみます。そのためには、generateコマンドを使って、ageを10で割った変数age10を作成し、age10をモデルに投入します。
. generate age10 = age/10
. stcox drug age10
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -83.551877
Iteration 2: log likelihood = -83.324008
Iteration 3: log likelihood = -83.323544
Refining estimates:
Iteration 0: log likelihood = -83.323544
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(2) = 33.18
Log likelihood = -83.323544 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
drug | .1048772 .0477017 -4.96 0.000 .0430057 .2557622
age10 | 3.114865 1.161372 3.05 0.002 1.499927 6.46857
------------------------------------------------------------------------------年齢10歳あたりハザード比が3.11 [1.16-1.50]だけ上昇することがわかりました。かなりのリスク上昇ですね。
