「stata」タグアーカイブ

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

Stata: discrete time hazard model

追跡期間が、連続変数ではなく、カテゴリ変数(1-2ヶ月、3-6ヶ月、7-12ヶ月等)である場合、通常のCox比例ハザードモデルではなく、discrete time hazard modelを用いた方がよい場合があります。

Discrete time hazard modelの概念については、下記の論文がよくまとまっています。

Richardson DB. Discrete time hazards models for occupational and environmental cohort analyses. Occup Envirn Med 2010;67:67-71

link関数をlogitにする場合とcomplementary log-logにする場合がありますが、Cox比例ハザードモデルに近いcomplementary log-log関数を利用した方がよさそうです。

2010年JAMAにdiscrete time hazard model with a complementary log-log linkを利用した論文が実際に発表されていました。

The NS, et al. Association of adolescent obesity with risk of severe obesity in adulthood. JAMA 2010; 304:2042–7

Stataのcloglogコマンドを使って、実際にdiscrete time hazard model with complementary log-log linkを行ってみましょう。

参考文献:Rabe-Hesketh S and Skrondal A. Multilevel and longitudinal model using Stata, second edition, Stata Corp LP

サンプルデータをダウンロードする

サンプルデータmortality.dtaをダウンロードします。1974-76年にガテマラで行われた15-49歳の女性とその子供を対象にした後ろ向き研究であり、母親の年齢や子供の出生順番が、子供の生存率に及ぼす影響を評価しています。1行が子供一人分のデータです。

. use http://www.stata-press.com/data/mlmus3/mortality,clear
(Child mortality in Guatemala)
  • kidid = 子供ID
  • momid = 母親ID
  • time = 子供の生存期間(月、最大60ヶ月)
  • death = 子供の死亡(1 = 死亡、0 = 打ち切り)
  • mage = 母親の出産時年齢
  • border = 子供の出生順番
  • p0014, p1523, p2435, p36up = 前の子供が生まれてからの期間(0-14ヶ月、15-23ヶ月、24-35ヶ月、36ヶ月以上)のダミー変数(兄の出生後16ヶ月で妹が生まれれば、p0014 = 0、p1523 = 1、p2435 = 0、p36up = 0)
  • pdead = 前の子供の死産
  • f0011とf1223 = 次の子供が生まれるまでの期間(0-11ヶ月、12-23ヶ月)のダミー変数

Discrete time hazard model用のデータシートを作成する

listコマンドを使って、momid==1のデータを確認してみます。

. list kidid momid time death mage border f0011 f1223 if momid==1, clean

        kidid   momid   time   death   mage   border   f0011   f1223  
   1.     101       1     60       0     27        4       0       0  
   2.     102       1     60       0     30        5       0       0  
   3.     103       1     60       0     33        6       0       1  
   4.     104       1     60       0     35        7       0       0  

子供生存期間timeの分布を確認してみましょう。最大値60が最も多く、不規則な分布です。

. histogram time, xsize(3.7) ysize(3.7)
(bin=41, start=.25, width=1.4573171)

追跡期間time(月)を、0〜1未満、1〜6未満、6〜12未満、12〜24未満、24〜60のカテゴリに区分したdiscrete変数を作成します。

. egen discrete = cut(time), at(0 1 6 12 24 61) icodes

. table discrete, contents(min time max time)

----------------------------------
 discrete |  min(time)   max(time)
----------+-----------------------
        0 |        .25         .25
        1 |          1           5
        2 |          6          11
        3 |         12          23
        4 |         24          60
----------------------------------

discrete 0〜4を1〜5に変換します。

. replace discrete = discrete + 1
(3,120 real changes made)

ltableコマンドを使って、discrete別に生存率を算出します。noadjustオプションを利用すると、Kaplan-Meier法による生存率が出力されるので、便利です。

. ltable discrete death, noadjust

                 Beg.                                 Std.
   Interval     Total   Deaths   Lost    Survival    Error     [95% Conf. Int.]
-------------------------------------------------------------------------------
    1     2      3120      109      9     0.9651    0.0033     0.9580    0.9710
    2     3      3002       92     94     0.9355    0.0044     0.9263    0.9436
    3     4      2816       66    126     0.9136    0.0051     0.9031    0.9230
    4     5      2624       94    238     0.8808    0.0059     0.8687    0.8919
    5     6      2292       42   2250     0.8647    0.0063     0.8518    0.8765
-------------------------------------------------------------------------------

expand変数を使って、追跡期間(discrete)の分だけレコードを複製します。

. list momid kidid time discrete death if momid==2, clean

        momid   kidid   time   discrete   death  
   5.       2     201     60          5       0  
   6.       2     202     12          4       1  
   7.       2     203     60          5       0  
   8.       2     204      1          2       1  

. expand discrete
(10,734 observations created)

. sort momid kidid

. list momid kidid time discrete death if momid==2, clean

         momid   kidid   time   discrete   death  
   21.       2     201     60          5       0  
   22.       2     201     60          5       0  
   23.       2     201     60          5       0  
   24.       2     201     60          5       0  
   25.       2     201     60          5       0  
   26.       2     202     12          4       1  
   27.       2     202     12          4       1  
   28.       2     202     12          4       1  
   29.       2     202     12          4       1  
   30.       2     203     60          5       0  
   31.       2     203     60          5       0  
   32.       2     203     60          5       0  
   33.       2     203     60          5       0  
   34.       2     203     60          5       0  
   35.       2     204      1          2       1  
   36.       2     204      1          2       1  

同一kididのレコードに、連番の変数(interval)を追加します。

. bysort kidid: generate interval = _n

. list momid kidid time discrete interval death if momid==2, clean

         momid   kidid   time   discrete   interval   death  
   21.       2     201     60          5          1       0  
   22.       2     201     60          5          2       0  
   23.       2     201     60          5          3       0  
   24.       2     201     60          5          4       0  
   25.       2     201     60          5          5       0  
   26.       2     202     12          4          1       1  
   27.       2     202     12          4          2       1  
   28.       2     202     12          4          3       1  
   29.       2     202     12          4          4       1  
   30.       2     203     60          5          1       0  
   31.       2     203     60          5          2       0  
   32.       2     203     60          5          3       0  
   33.       2     203     60          5          4       0  
   34.       2     203     60          5          5       0  
   35.       2     204      1          2          1       1  
   36.       2     204      1          2          2       1  

従属変数y=0を作成します。最終追跡時(discrete == interval)ではy=deathに設定します。

. generate y = 0

. replace y = death if interval == discrete
(403 real changes made)

. list momid kidid discrete interval death y if momid==2, clean

         momid   kidid   discrete   interval   death   y  
   21.       2     201          5          1       0   0  
   22.       2     201          5          2       0   0  
   23.       2     201          5          3       0   0  
   24.       2     201          5          4       0   0  
   25.       2     201          5          5       0   0  
   26.       2     202          4          1       1   0  
   27.       2     202          4          2       1   0  
   28.       2     202          4          3       1   0  
   29.       2     202          4          4       1   1  
   30.       2     203          5          1       0   0  
   31.       2     203          5          2       0   0  
   32.       2     203          5          3       0   0  
   33.       2     203          5          4       0   0  
   34.       2     203          5          5       0   0  
   35.       2     204          2          1       1   0  
   36.       2     204          2          2       1   1  

これでdiscrete time hazard model用のデータセットがほぼ完成です。とりあえず、各intervalにおける死亡率を算出してみます。

. tabulate interval y, row

+----------------+
| Key            |
|----------------|
|   frequency    |
| row percentage |
+----------------+

           |           y
  interval |         0          1 |     Total
-----------+----------------------+----------
         1 |     3,011        109 |     3,120 
           |     96.51       3.49 |    100.00 
-----------+----------------------+----------
         2 |     2,910         92 |     3,002 
           |     96.94       3.06 |    100.00 
-----------+----------------------+----------
         3 |     2,750         66 |     2,816 
           |     97.66       2.34 |    100.00 
-----------+----------------------+----------
         4 |     2,530         94 |     2,624 
           |     96.42       3.58 |    100.00 
-----------+----------------------+----------
         5 |     2,250         42 |     2,292 
           |     98.17       1.83 |    100.00 
-----------+----------------------+----------
     Total |    13,451        403 |    13,854 
           |     97.09       2.91 |    100.00 

Complementary log-log modelを作成するためには、intervalのダミー変数を作成する必要があります。tabulate, generate()で作成します。

. tabulate interval, generate(int)

   interval |      Freq.     Percent        Cum.
------------+-----------------------------------
          1 |      3,120       22.52       22.52
          2 |      3,002       21.67       44.19
          3 |      2,816       20.33       64.52
          4 |      2,624       18.94       83.46
          5 |      2,292       16.54      100.00
------------+-----------------------------------
      Total |     13,854      100.00

. list momid kidid interval int1-int5 death y if momid==2, clean

         momid   kidid   interval   int1   int2   int3   int4   int5   death   y  
   21.       2     201          1      1      0      0      0      0       0   0  
   22.       2     201          2      0      1      0      0      0       0   0  
   23.       2     201          3      0      0      1      0      0       0   0  
   24.       2     201          4      0      0      0      1      0       0   0  
   25.       2     201          5      0      0      0      0      1       0   0  
   26.       2     202          1      1      0      0      0      0       1   0  
   27.       2     202          2      0      1      0      0      0       1   0  
   28.       2     202          3      0      0      1      0      0       1   0  
   29.       2     202          4      0      0      0      1      0       1   1  
   30.       2     203          1      1      0      0      0      0       0   0  
   31.       2     203          2      0      1      0      0      0       0   0  
   32.       2     203          3      0      0      1      0      0       0   0  
   33.       2     203          4      0      0      0      1      0       0   0  
   34.       2     203          5      0      0      0      0      1       0   0  
   35.       2     204          1      1      0      0      0      0       1   0  
   36.       2     204          2      0      1      0      0      0       1   1  

Discrete time hazard model with complementary log-log model

cloglogコマンドを利用して、complementary log-logモデルを作成します。まずは暴露因子を含まないモデルを作成します。オプションに切片0およびハザード比表示を指定します。

. cloglog y int1-int5, noconstant eform

Iteration 0:   log likelihood = -1811.7312  
Iteration 1:   log likelihood = -1811.6791  
Iteration 2:   log likelihood = -1811.6791  

Complementary log-log regression                Number of obs     =     13,854
                                                Zero outcomes     =     13,451
                                                Nonzero outcomes  =        403

                                                Wald chi2(5)      =    4943.80
Log likelihood = -1811.6791                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
           y |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        int1 |   .0355608   .0034063   -34.83   0.000     .0294738    .0429048
        int2 |   .0311257   .0032452   -33.28   0.000      .025373    .0381826
        int3 |   .0237165   .0029194   -30.40   0.000     .0186326    .0301876
        int4 |   .0364806   .0037629   -32.10   0.000     .0298031    .0446541
        int5 |   .0184946   .0028538   -25.86   0.000     .0136678    .0250259
------------------------------------------------------------------------------

次に、年齢・性別補正モデルを作成します。

. cloglog y int1-int5 sex mage, noconstant eform

Iteration 0:   log likelihood =  -1811.344  
Iteration 1:   log likelihood = -1811.2917  
Iteration 2:   log likelihood = -1811.2917  

Complementary log-log regression                Number of obs     =     13,854
                                                Zero outcomes     =     13,451
                                                Nonzero outcomes  =        403

                                                Wald chi2(7)      =    4941.87
Log likelihood = -1811.2917                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
           y |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        int1 |   .0297238   .0067244   -15.54   0.000     .0190781    .0463097
        int2 |   .0260221   .0059782   -15.88   0.000     .0165879    .0408218
        int3 |   .0198284    .004736   -16.41   0.000     .0124158    .0316663
        int4 |   .0304998   .0069904   -15.23   0.000     .0194629    .0477956
        int5 |   .0154664   .0039604   -16.28   0.000     .0093633    .0255477
         sex |   1.026905   .1024613     0.27   0.790     .8445013    1.248707
        mage |    1.00634   .0076432     0.83   0.405     .9914704    1.021432
------------------------------------------------------------------------------

母親によるクラスターであると考えれば、Stataではvec(cluster cluster_var)オプションを利用して、robust standard error for clustered dataを求めることができます。

. cloglog y int1-int5 sex mage, noconstant eform vce(cluster momid)

Iteration 0:   log pseudolikelihood =  -1811.344  
Iteration 1:   log pseudolikelihood = -1811.2917  
Iteration 2:   log pseudolikelihood = -1811.2917  

Complementary log-log regression                Number of obs     =     13,854
                                                Zero outcomes     =     13,451
                                                Nonzero outcomes  =        403

                                                Wald chi2(7)      =    4328.62
Log pseudolikelihood = -1811.2917               Prob > chi2       =     0.0000

                                (Std. Err. adjusted for 851 clusters in momid)
------------------------------------------------------------------------------
             |               Robust
           y |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        int1 |   .0297238    .007154   -14.61   0.000     .0185453    .0476403
        int2 |   .0260221   .0069143   -13.73   0.000     .0154586     .043804
        int3 |   .0198284   .0053308   -14.58   0.000      .011707    .0335837
        int4 |   .0304998   .0077332   -13.76   0.000     .0185557    .0501323
        int5 |   .0154664   .0043305   -14.89   0.000     .0089343    .0267745
         sex |   1.026905   .1055175     0.26   0.796     .8395895    1.256012
        mage |    1.00634   .0087612     0.73   0.468     .9893139    1.023659
------------------------------------------------------------------------------

母親をランダム変数として組み込んだ、random-interceptモデルを作成します。xtcloglogコマンドを利用します。

. xtcloglog y int1-int5 sex mage, noconstant eform i(momid)

Fitting comparison model:

Iteration 0:   log likelihood =  -1811.344  
Iteration 1:   log likelihood = -1811.2917  
Iteration 2:   log likelihood = -1811.2917  

Fitting full model:

tau =  0.0     log likelihood = -1911.7661
tau =  0.1     log likelihood = -1893.1909
tau =  0.2     log likelihood = -1874.1809
tau =  0.3     log likelihood = -1856.4094
tau =  0.4     log likelihood = -1842.2824
tau =  0.5     log likelihood = -1834.7269
tau =  0.6     log likelihood = -1837.5067

Iteration 0:   log likelihood = -1841.9925  
Iteration 1:   log likelihood = -1808.6147  
Iteration 2:   log likelihood = -1806.9026  
Iteration 3:   log likelihood = -1806.7965  
Iteration 4:   log likelihood = -1806.7963  
Iteration 5:   log likelihood = -1806.7963  

Random-effects complementary log-log model      Number of obs     =     13,854
Group variable: momid                           Number of groups  =        851

Random effects u_i ~ Gaussian                   Obs per group:
                                                              min =          1
                                                              avg =       16.3
                                                              max =         40

Integration method: mvaghermite                 Integration pts.  =         12

                                                Wald chi2(7)      =    2450.04
Log likelihood  = -1806.7963                    Prob > chi2       =     0.0000

------------------------------------------------------------------------------
           y |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        int1 |   .0277072   .0065873   -15.08   0.000     .0173869    .0441532
        int2 |   .0245443   .0059068   -15.40   0.000     .0153144    .0393369
        int3 |   .0188305    .004693   -15.94   0.000     .0115538    .0306903
        int4 |   .0291856   .0069987   -14.74   0.000     .0182411    .0466969
        int5 |   .0149082   .0039561   -15.85   0.000     .0088623    .0250786
         sex |   1.033765   .1049692     0.33   0.744     .8472088    1.261402
        mage |   1.003026   .0081413     0.37   0.710     .9871959    1.019111
-------------+----------------------------------------------------------------
    /lnsig2u |  -1.254711   .3784523                     -1.996464   -.5129578
-------------+----------------------------------------------------------------
     sigma_u |   .5340022   .1010472                      .3685305    .7737713
         rho |   .1477433   .0476529                      .0762683    .2668511
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.
LR test of rho=0: chibar2(01) = 8.99                   Prob >= chibar2 = 0.001

Stata: Kaplan-Meir曲線

Stataの特徴は比較的簡単に美しいグラフの作成できることです。Stataのサンプルデータdrug2を使って、Kaplan-Meier曲線(白黒)を作成してみましょう。

use https://www.stata-press.com/data/r16/drug2
sts graph, by(drug)

このKaplan-Meier曲線に設定を追加して、下図を作成します。

この図を作成するためには、以下の魔法を唱える必要があります。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)    xtitle(" " "Observational period (year)", size(small))  xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small))  xtick(0(6)48)   title("A", size("medlarge") span position(11))    plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black))    legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 "Non-user (n=20)" 2 "User (n=28)") size(small) title("Drug use", size("small") color("black")))    risktable(0(12)48, size(small) order(1 "Non-user" 2 "User") rowtitle(, justification(left)) title(, size(small)))    text(1 48 "P<0.001", placement(sw) size(small))    graphregion(color(white) lcolor(white))

それでは、この魔法が作られた過程を確認していきましょう。

グラフのサイズを変更する

A4サイズ(8.27 × 11.69インチ)に2 x 3個のグラフが収まるサイズにしたいので、図のサイズを3.7 × 3.7 インチに設定します。

sts graph, xsize(3.7) ysize(3.7)
sts graph, xsize(3.7) ysize(3.7)

カテゴリ別の累積アウトカム発症率を算出する

Y軸を、累積生存率ではなく、累積アウトカム発症率に変更し、durgカテゴリ別のKaplan-Meier曲線を描きます。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure
sts graph, xsize(3.7) ysize(3.7) by(drug) failure

Y軸を調整する

Y軸名(ytitle)を”Cumulative probability of incidence of outcome”にします。デフォルトの文字サイズは大きいので、小さくします(small)。デフォルト設定のままだと軸と軸名がくっつきすぎるので、Y軸名に空白行を1行挿入して、軸と軸名の間隔を少し広げます。

Y軸の数値(ylabel)の向き(angle)を水平にして、その間隔を0.2にしますが、目盛り線(tick)の間隔は0.1にします。ラベルの文字を少し小さくして(small)、小数点1桁表示にします。ガイド線(グリッド)を消します。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)
sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle(“Cumulative probability of incidence of outcome” ” “, size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6)

X軸を調整する

X軸名(xtitle)を”Observational period (year)”にします。デフォルトの文字サイズは大きいので、小さくします(small)。デフォルト設定のままだと軸と軸名がくっつきすぎるので、X軸名に空白行を1行挿入して、軸と軸名の間隔を少し広げます。

追跡期間データの単位は月(0〜48)ですが、X軸の単位を年に変換するために、第0, 12, 24, 36, 48月に0, 1, 2, 3, 4を表示させます。目盛り線(tick)の間隔を6ヶ月にします。ラベルの文字を少し小さくします(small)。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)    xtitle(" " "Observational period (year)", size(small))  xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small))  xtick(0(6)48)
sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle(“Cumulative probability of incidence of outcome” ” “, size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6) xtitle(” ” “Observational period (year)”, size(small)) xlabel(0 “0” 12 “1” 24 “2” 36 “3” 48 “4”, labsize(small)) xtick(0(6)48)

タイトルを変更する

グラフのタイトル”A”を、グラフ全体の左上にやや大きめ(medlarge)で追加します。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)    xtitle(" " "Observational period (year)", size(small))  xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small))  xtick(0(6)48)   title("A", size("medlarge") span position(11))    
sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle(“Cumulative probability of incidence of outcome” ” “, size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6) xtitle(” ” “Observational period (year)”, size(small)) xlabel(0 “0” 12 “1” 24 “2” 36 “3” 48 “4”, labsize(small)) xtick(0(6)48) title(“A”, size(“medlarge”) span position(11))

生存曲線のレイアウトを変更する

drugカテゴリ別の生存曲線のレイアウトを変更します。カラーイラストは論文掲載料が追加されることが多いので、生存曲線の色を黒に統一して、drug=0を点線、drug=1を実線に設定します。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle("Cumulative probability of incidence of outcome" " ", size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6) xtitle(" " "Observational period (year)", size(small)) xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small)) xtick(0(6)48) title("A", size("medlarge") span position(11)) plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black))
sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle(“Cumulative probability of incidence of outcome” ” “, size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6) xtitle(” ” “Observational period (year)”, size(small)) xlabel(0 “0” 12 “1” 24 “2” 36 “3” 48 “4”, labsize(small)) xtick(0(6)48) title(“A”, size(“medlarge”) span position(11)) plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black))

凡例を修正する

凡例を、生存曲線のプロット領域内(ring(0))の右下(5時方向)に、縦1列(cols(1))で表示する。枠線を白に設定して見えなくする(region(lcolor(white)))。カテゴリ名を”Drug use”にして、drugカテゴリ0には”Non-user (n=20)”を、drugカテゴリ1には”User (n=28)”を、黒文字で小さめに表示する。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)    xtitle(" " "Observational period (year)", size(small))  xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small))  xtick(0(6)48)   title("A", size("medlarge") span position(11))    plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black))    legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 "Non-user (n=20)" 2 "User (n=28)") size(small) title("Drug use", size("small") color("black")))
sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle(“Cumulative probability of incidence of outcome” ” “, size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6) xtitle(” ” “Observational period (year)”, size(small)) xlabel(0 “0” 12 “1” 24 “2” 36 “3” 48 “4”, labsize(small)) xtick(0(6)48) title(“A”, size(“medlarge”) span position(11)) plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black)) legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 “Non-user (n=20)” 2 “User (n=28)”) size(small) title(“Drug use”, size(“small”) color(“black”)))

Number at riskを追加する

Number at riskを12ヶ月毎に追加する。drugカテゴリ名は小さめ(small)に左揃えで”Non-user”と”User”と表示する。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)    xtitle(" " "Observational period (year)", size(small))  xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small))  xtick(0(6)48)   title("A", size("medlarge") span position(11))    plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black))    legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 "Non-user (n=20)" 2 "User (n=28)") size(small) title("Drug use", size("small") color("black")))    risktable(0(12)48, size(small) order(1 "Non-user" 2 "User") rowtitle(, justification(left)) title(, size(small)))
sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle(“Cumulative probability of incidence of outcome” ” “, size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6) xtitle(” ” “Observational period (year)”, size(small)) xlabel(0 “0” 12 “1” 24 “2” 36 “3” 48 “4”, labsize(small)) xtick(0(6)48) title(“A”, size(“medlarge”) span position(11)) plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black)) legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 “Non-user (n=20)” 2 “User (n=28)”) size(small) title(“Drug use”, size(“small”) color(“black”))) risktable(0(12)48, size(small) order(1 “Non-user” 2 “User”) rowtitle(, justification(left)) title(, size(small)))

テキスト(P値)を追加する

XY座標(1, 48)の左下(sw = southwest)にテキスト”P<0.001″を追加する。

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)    xtitle(" " "Observational period (year)", size(small))  xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small))  xtick(0(6)48)   title("A", size("medlarge") span position(11))    plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black))    legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 "Non-user (n=20)" 2 "User (n=28)") size(small) title("Drug use", size("small") color("black")))    risktable(0(12)48, size(small) order(1 "Non-user" 2 "User") rowtitle(, justification(left)) title(, size(small)))    text(1 48 "P<0.001", placement(sw) size(small))
sts graph, xsize(3.7) ysize(3.7) by(drug) failure ytitle(“Cumulative probability of incidence of outcome” ” “, size(small)) ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid) ytick(0.0(0.1)0.6) xtitle(” ” “Observational period (year)”, size(small)) xlabel(0 “0” 12 “1” 24 “2” 36 “3” 48 “4”, labsize(small)) xtick(0(6)48) title(“A”, size(“medlarge”) span position(11)) plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black)) legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 “Non-user (n=20)” 2 “User (n=28)”) size(small) title(“Drug use”, size(“small”) color(“black”))) risktable(0(12)48, size(small) order(1 “Non-user” 2 “User”) rowtitle(, justification(left)) title(, size(small))) text(1 48 “P<0.001”, placement(sw) size(small))

テキストの配置(placement)には、sw以外にも下記を設定できます。テキストの表示位置に応じて適宜変更して下さい。

  • c = centered on the point, vertically and horizontally
  • n = above the point, centered
  • ne = above and to the right of the point
  • e = right of the point, vertically centered
  • se = below and to the right of the point
  • s = below point, centered
  • sw = below and to the left of the point
  • w = left of the point, vertically centered
  • nw = above and to the left of the point

Stata公式マニュアル https://www.stata.com/manuals/g-3added_text_options.pdf

背景色を変更する

グラフ領域の背景を白色に設定して完成!

sts graph, xsize(3.7) ysize(3.7) by(drug) failure    ytitle("Cumulative probability of incidence of outcome" " ", size(small))  ylabel(0.0(0.2)1, angle(0) labsize(small) format(%9.1f) nogrid)  ytick(0.0(0.1)0.6)    xtitle(" " "Observational period (year)", size(small))  xlabel(0 "0" 12 "1" 24 "2" 36 "3" 48 "4", labsize(small))  xtick(0(6)48)   title("A", size("medlarge") span position(11))    plot1opt(lpattern(dash)) plot2opt(lpattern(solid)) plotopt(lcolor(black))    legend(ring(0) position(5) cols(1) region(lcolor(white)) order(1 "Non-user (n=20)" 2 "User (n=28)") size(small) title("Drug use", size("small") color("black")))    risktable(0(12)48, size(small) order(1 "Non-user" 2 "User") rowtitle(, justification(left)) title(, size(small)))    text(1 48 "P<0.001", placement(sw) size(small))    graphregion(color(white) lcolor(white))