山本陵平 のすべての投稿

Clinical Journal Club 2022年8月17日

オオサンショウウオ先生の糖尿病論文で学ぶ医療統計セミナー 

第12講 臨床試験のアウトカム(担当:松村)

BMJ抄読会

Liang C. et al. Infertility, recurrent pregnancy loss, and risk of stroke: pooled analysis of individual patient data of 618 851 women. BMJ 2022; 377 doi: https://doi.org/10.1136/bmj-2022-070603
第3査読者と第4査読者のコメントおよび著者の対応

https://youtu.be/WQENR-hvfWs(参加者限定)

Linux: nkfを使ってshiftjisをUTF-8に変換する

macやlinuxを利用していると、microsoft excelで作成されたcsvファイルの文字コードがshiftjisのため文字化けして読めない状態によく遭遇します。そういう時はnkf (network kanji filter)で文字コードを変換してしまいましょう。

Rocky Linux 8にインストールする

terminalで下記コマンドを入力すればインストールできます。

[root@localhost]# dnf install nkf -y

nkfでUTF-8に文字コードを変換する

ユーザーのホームディレクトリのCSVファイルの文字コードを知りたい場合

[username@localhost]# nkf --guess ~/*csv

上記をまとめて UTF-8に変更したい場合

[username@localhost]# nkf --overwrite -wLu ~/*csv

再度–guessすれば無事変換されていることを確認できます。

[username@localhost]# nkf --guess ~/*csv

Linux: Rocky Linux 8のインストール

UbuntuとRの相性が悪いみたいで、CentOSでは機能していたRのプログラムが全く機能しません。CentOSの開発者が新たに開発したRockey Linuxに移行することにしました。

Rocky Linxuのホームページ https://rockylinux.org/ja/

Rockey Linuxのインストーラー.isoをダウンロードして、USBメディア等に保存して、インストール用メディアを作成します。作成方法はGoogle等で適宜検索して下さい。

インストール

BIOSを起動して、インストーラーが保存されているUSBメモリから起動するように設定します。TOSHIBA Dynabook R734の場合、起動ボタンを押してからF12を数回押せば、起動用デバイスの選択画面に移行します。

Test this media & install Rochy Linux 8 を選択します。

WELCOME TO ROCKY LINUX8において、言語としてEnglish > English (United States) を選択します。

INSTALLATION SUMMARYにおいて、下記の設定を行います。

  1. Root Password を選択して、パスワードを入力します。
  2. User Creationを選択して、普段のログイン時に利用するユーザー名とパスワードを設定します。”Make this user administrator”にもチェックしておくと、
  3. Installation Destinationを選択する。
    • Local Standard Disksからインストール先を選択する。
    • クリーンインストールを行う場合、Storage Configurationの”I would like to make additional space”をチェックし、左上の”Done”をクリックして、RECLAIM DISK SPACEが表示する。右下の”Delete all”をクリックしてから、”Reclaim space”をクリックする。
  4. Begin Installationをクリックして、インストールを開始する

インストールが完了したら、右下の”Reboot System”をクリックします。

初期設定

LICENSINGを選択して、左下の”I accept the liscence agreement”をチェックして、左上の”Done”をクリックします。右下の”FINISH CONFIGURATION”をクリックします。

ユーザーを選択してパスワードを入力します。

Welcome画面でメインで利用する言語を選択します。私は”English”を選択しました。右上の”Next”をクリックします。

利用しているキーボードを選択します。私は英語配列キーボードを利用しているので”English (US)”を選択しました。右上の”Next”をクリックします。

Wi-Fi画面で接続するルーターを利用します。私は有線LANしか利用しないので、右上の”Skip”をクリックしました。

Privacy画面で現在位置情報の利用の有無を設定します。私は”OFF”に設定しました。右上の”Next”をクリックします。

Ready to Go画面で”Start Using Rocky Linux”ボタンをクリックします。

以上でインストールは完了です。

Impact Factor 2021

Impact Factor 2020

一般・疫学系

202.731 Lacnet

176.079 New England Journal of Medicine

157.335 JAMA

93.467 BMJ

72.427 Lancet Publich Health

51.598 Annals of Internal Medicine

44.460 JAMA Internal Medicine

37.746 Lancet Child & Adolescent Health

17.033 eClinicalMedicine

16.882 Canadian Medical Association Journal

14.040 QJM

13.366 JAMA Network Open

13.068 Journal of Internal Medicine

12.776 Medical Journal of Australia

11.613 PLoS Medicine

11.150 BMC Medicine

11.104 Mayo Clinic Proceedings (≤3000単語, ≤5図表)

9.685 International Journal of Epidemiology

8.559 The Lancet Regional Health – Western Pacific

6.604 American Journal of Preventive Medicine

6.461 Frontiers in Public Health

5.919 Epidemiology and Health ($800, 韓国疫学会)

5.100 International Journal of Public Health (open)

5.058 Frontiers in Medicine

4.996 Scientific Report (掲載料$1990)ref

4.964 Journal of Clinical Medicine (MDPI)

4.860 Epidemiology

4.637 Preventive Medicine

4.614 International Journal of Environmental Research and Public Health (MDPI)

4.395 Environmental Health and Preventive Medicine (日本衛星学会)

4.135 BMC public health (open)

3.752 PloS One (open)

3.571 Journal of the American Nutrition Association (Journal of American College of Nutrition)

3.160 Healthcare (MDPI)

3.017 BMJ Open

2.956 American Journal of Health Promotion

2.395 Journal of American College Health

Digital Health系

36.615 Lancet Digital Health

7.076 Journal of Medical Internet Research

腎臓系

18.998 Kidney International

14.978 Journal of American Society of Nephrology

11.072 American Journal of Kidney Disease (<3500単語)ref

10.614 Clinical Journal of American Society of Nephrology (<3000単語) (ref)

7.186 Nephrology Dialysis Transplantation (ref)

6.234 Kidney International Report (APC $2100, <4000単語)

5.860 Clinical Kidney Journal (APC $2100, <3000単語)

4.605 American Journal of Nephrology

6.055 Frontiers in Endocrinology (Renal Endocrinlogy section, APC $2950)

4.393 Journal of Nephrology (イタリア腎臓学会雑誌, <3000単語)

4.360 CardioRenal Medicine

4.354 Journal of Renal Nutrition (NKF)

4.172 Kidney Research and Clinical Practice (APC $500, 韓国腎臓学会雑誌, <4000単語, ≤40文献)

3.651 Pediatric Nephrology

3.457 Nephron

3.096 Kidney & blood pressure research

3.084 Nephrologia

3.000 Kidney Diseases (Karger, APC $1710)

2.879 Peritoneal Dialysis International

2.617 Clinical and Experimental Nephrology

2.585 BMC Nephrology

2.358 Nephrology

糖尿病・代謝系

44.867 Lancet diabetes & endocrinology

17.152 Diabetes care

10.460 Diabetologia

10.460 Diabetes

8.949 Cardiovascular diabetology

8.254 Diabetes & metabolism

8.180 Diabetes research and clinical practice

6.408 Diabetes obesity & metaboblism

栄養系

9.298 Obesityref

6.408 American Journal of Clinical Nutrition

8.023 Current Obesity Reports

7.643 Clinical Nutrition

6.706 Nutrients

6.410 Diabetes Obesity & Metabolism

6.590 Frontiers in Nutrition

5.923 Annals of Nutrition and Metabolism

5.551 International Journal of Obesity (ref)

5.234 Journal of the Academy of Nutrition and Dietetics

5.214 Obesity Research & Clinical Practice

4.865 European Journal of Nutrition

4.807 Obesity Fact

4.725 Nutrition and Diabetes

4.687 Journal of Nutrition

4.344 Nutrition Journal

老年医学系

7.538 Journal of the American Geriatric Society

歯科系

24.897 International Journal of Oral Science (Springer)

12.239 Periodontology 2000 (Wiley)

8.924 Journal of Dental Research (SAGE)

7.478 Journal of Clnical Periodontology (Wiley)

6.468 Japanese Dental Science Review (Elsevier)

4.068 Oral diseases (Wiley)

3.608 Clinical Oral Investigation (Springer)

2.489 Community Dentistry and Oral Epidemiology (Wiley)

アルコール

3.346 Journal of Studies on Alcohol and Drugs

2.362 Substance Use and Misuse

Impact factor 2020

一般・疫学系

8.989 Journal of Internal Medicine

8.483 JAMA Network Open (採択率16%)

7.616 Mayo Clinic Proceedings (≤3000単語, ≤5図表)

5.043 American Journal of Preventive Medicine

4.379 Scientific Report (掲載料$1990)ref

4.018 Preventive Medicine

3.709 Frontiers in Public Health

3.674 Environmental Health and Preventive Medicine

3.390 International Journal of Environmental Research and Public Health (MDPI)

3.38 International Journal of Public Health (open)

3.295 BMC public health (open)

3.282 Epidemiology and Health ($800, 韓国疫学会)

3.240 PloS One (open)

3.087 Journal of American College Health

2.870 American Journal of Health Promotion

2.692 BMJ Open

2.645 Healthcare (MDPI)

腎臓系

10.612 Kidney International

10.121 Journal of American Society of Nephrology

8.860 American Journal of Kidney Disease (<3500単語, 投稿時施設入力不要)ref

8.237 Clinical Journal of American Society of Nephrology (<3000単語) (ref)

5.992 Nephrology Dialysis Transplantation (91%) (ref)

5.555 Frontiers in Endocrinology (Renal Endocrinlogy section, APC $2950)

4.452 Clinical Kidney Journal (APC $2100, <3000単語, 著者全員の利益相反提出)

4.164 Kidney International Report (APC $2100, <4000単語)

3.902 Journal of Nephrology (イタリア腎臓学会雑誌, <3000単語)

3.754 American Journal of Nephrology

3.714 Pediatric Nephrology

3.667 Kidney Research and Clinical Practice (APC $500, 韓国腎臓学会雑誌, <4000単語, ≤40文献)

3.655 Journal of Renal Nutrition (NKF)

3.222 Kidney Diseases (Karger, APC $1710, 採択率37%)

2.847 Nephron

2.801 Clinical and Experimental Nephrology (62%)

2.687 Kidney & blood pressure research

2.506 Nephrology

2.388 BMC Nephrology

2.238 Nephron Exp Nephrol

1.511 Exp Nephrol

0.975 Clin Nephrol

糖尿病・代謝系

32.069 Lancet diabetes & endocrinology

19.112 Diabetes care

10.122 Diabetologia

9.951 Cardiovascular diabetology

9.461 Diabetes

8.694 Metabolsim – clinical and experimental

6.664 European Journal of Endocrionology

6.577 Diabetes obesity & metaboblism

6.041 Diabetes & metabolism

5.958 Journal of clinical endocrinology & metabbolism

5.829 Aging Male

5.602 Diabetes research and clinical practice

5.376 Diabetes & Metabolism Journal

4.359 Diabetic Medicine

4.280 Acta Diabetologica

栄養系

7.324 Clinical Nutrition

7.045 American Journal of Clinical Nutrition

6.919 Current Obesity Reports

6.577 Diabetes Obesity & Metabolism

6.576 Frontiers in Nutrition

5.719 Nutrients (87%)

5.614 European Journal of Nutrition

5.095 Nutrition and Diabetes

5.095 International Journal of Obesity (ref)

5.002 Obesityref

4.910 Journal of the Academy of Nutrition and Dietetics

4.599 Journal of Nutrition

3.942 Obesity Fact

3.374 Annals of Nutrition and Metabolism

3.359 Nutrition Journal

2.288 Obesity Research & Clinical Practice

睡眠系

5.849 Sleep

5.346 Nature and Science of Sleep

4.450 Sleep Health

4.062 Journal of Clinical Sleep Medicine

3.981 Journal of sleep research

3.492 Sleep Medicine

2.964 Behavioral Sleep Medicine

2.816 Sleep and Breathing

1.186 Sleep and Biological Rhythms

アルコール系

6.525 Addiction

4.492 Drug and Alcohol Dependence

3.455 Alcoholism: clinical and experimental research

2.582 Journal of Studies on Alcohol and Drugs

歯科系

6.116 Journal of Dental Research

4.642 Journal of Prosthodontic Research

3.511 Oral Diseases

3.383 Community Dentistry and Oral Epidemiology

2.980 Gerodontology

2.757 BMC Oral Health

老年医学系

2.730 Geriatrics & Gerontology International

日本の学会雑誌

7.527 Journal of Gastroenterology

4.928 Journal of Atherosclerosis and Thrombosis

4.823 International Immunology

4.228 Hepatology Research

4.005 Journal of Dermatology

3.872 Hypertension Research (日本高血圧学会)

3.674 Environmental Health and Preventive Medicine (日本衛生学会)

3.374 Annals of Nutrition and Metabolism

3.369 International Journal of Urology

3.211 Journal of Epidemiology (日本疫学会)

3.159 Journal of Cardiology

3.023 Modern Rheumatology (日本リウマチ学会)

2.993 Circulation Journal

2.801 Clinical and Experimental Nephrology (日本腎臓学会)

2.730 Geriatrics & Gerontology International

2.612 Journal of Bone and Mineral Metabolism

2.549 Surgery Today

2.490 International Journal of Hematology

2.349 Endocrine Journal

1.848 Tohoku Journal of Experimental Medicine

1.730 Journal of Obstetrics and Gynaecology Research

1.726 Therapeutic Apheresis and Dialysis

1.271 Internal Medicine (日本内科学会)

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))

ホームページ、始めました。

大阪大学キャンパスライフ支援センターの山本陵平です。最近、いろいろなところで疫学研究の基本を教える機会が増えてきたので、ホームページにまとめて公開することにしました。統計、研究デザイン、バイアス、プログラミング等の疫学研究を行う上で役に立ちそうな内容を、書き連ねていこうと思います。