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
sts graph, by(drug) xsize(3.7) ysize(3.7)

薬剤介入群の累積生存率は、プラセボ群よりもかなり良いみたいです。二群の累積生存率を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.0000

Log-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)
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]だけ上昇することがわかりました。かなりのリスク上昇ですね。

次の記事

Impact factor 2020