第13講 ランダム化（大槻）
Fu E. Timing of dialysis initiation to reduce mortality and cardiovascular events in advanced chronic kidney disease: nationwide cohort study. BMJ 2021;375:e066306. doi: https://doi.org/10.1136/bmj-2021-066306, Published on 2021-11-29
Introduction, methos, and results of oritinal article
第12講 臨床試験のアウトカム（担当：松村）
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（参加者限定）
macやlinuxを利用していると、microsoft excelで作成されたcsvファイルの文字コードがshiftjisのため文字化けして読めない状態によく遭遇します。そういう時はnkf (network kanji filter)で文字コードを変換してしまいましょう。
terminalで下記コマンドを入力すればインストールできます。
[root@localhost]# dnf install nkf -yユーザーのホームディレクトリのCSVファイルの文字コードを知りたい場合
[username@localhost]# nkf --guess ~/*csv上記をまとめて UTF-8に変更したい場合
[username@localhost]# nkf --overwrite -wLu ~/*csv再度–guessすれば無事変換されていることを確認できます。
[username@localhost]# nkf --guess ~/*csvUbuntuと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において、下記の設定を行います。
インストールが完了したら、右下の”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”ボタンをクリックします。
以上でインストールは完了です。
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
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
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を用いて、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) には、下記の変数が含まれています。
生存解析を始める前に、stsetコマンドで追跡期間とアウトカムを設定します。
. stset studytime, failure(died==1)
failure event: died == 1
obs. time interval: (0, studytime]
exit on or before: failure
------------------------------------------------------------------------------
48 total observations
0 exclusions
------------------------------------------------------------------------------
48 observations remaining, representing
31 failures in single-record/single-failure data
744 total analysis time at risk and under observation
at risk from t = 0
earliest observed entry t = 0
last observed exit t = 39対象48人、アウトカム発症31人であることがわかります。
まずは、薬剤介入群とプラセボ群のKaplan-Meier曲線を作成してみましょう。sts graphコマンドを使います。図のサイズは、A4（8.27 × 11.69インチ）に2 x 3個のグラフが収まるサイズ（3.7 × 3.7 インチ）に設定します。
. sts graph, by(drug) xsize(3.7) ysize(3.7)
failure _d: died
analysis time _t: studytime
薬剤介入群の累積生存率は、プラセボ群よりもかなり良いみたいです。二群の累積生存率をLog-rank検定で比較してみましょう。sts testコマンドを利用します。
. sts test drug
failure _d: died
analysis time _t: studytime
Log-rank test for equality of survivor functions
| Events Events
drug | observed expected
------+-------------------------
0 | 19 7.25
1 | 12 23.75
------+-------------------------
Total | 31 31.00
chi2(1) = 28.27
Pr>chi2 = 0.0000Log-rank検定のP値 (Pr>chi2) < 0.001ですので、薬剤介入群とプラセボ群の累積生存率は、統計学的に有意に異なると結論づけることができます。
それでは、Cox比例ハザードモデルを用いて、プラセボ群に対する薬剤介入群のハザード比を算出してみましょう。stcoxコマンドを利用します。
. stcox drug
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -88.254734
Iteration 2: log likelihood = -88.001551
Iteration 3: log likelihood = -88.00019
Refining estimates:
Iteration 0: log likelihood = -88.00019
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(1) = 23.82
Log likelihood = -88.00019 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
drug | .1327581 .0584002 -4.59 0.000 .0560555 .3144157
------------------------------------------------------------------------------プラセボ群（drug = 0）に対する薬剤介入群（drug = 1）のハザード比は、0.1327581 [0.0560555 – 0.3144157]であり、P = 0.000です。stcoxでは、変数の最も小さい値が自動的にreferenceに設定されます。任意の値をreferenceに設定したければ、変数の前にb.任意の値を設定します。薬剤介入群（drug = 1）をreferenceに設定したければ、b1.drugです。
. stcox b1.drug
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -88.254734
Iteration 2: log likelihood = -88.001551
Iteration 3: log likelihood = -88.00019
Refining estimates:
Iteration 0: log likelihood = -88.00019
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(1) = 23.82
Log likelihood = -88.00019 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
0.drug | 7.532496 3.313542 4.59 0.000 3.180502 17.83948
------------------------------------------------------------------------------続いて、本モデルが比例ハザード性を満たしているかを確認します。比例ハザード性を確認する方法はいろいろありますが、今回はSchoenfeld残差を利用した方法で確認してみましょう。estat phtestコマンドを利用して、直前に作成したCox比例ハザードモデルの比例ハザード性を検定します。
. estat phtest, detail
Test of proportional-hazards assumption
Time: Time
----------------------------------------------------------------
| rho chi2 df Prob>chi2
------------+---------------------------------------------------
0.drug | 0.00948 0.00 1 0.9583
1b.drug | . . 1 .
------------+---------------------------------------------------
global test | 0.00 1 0.9583
----------------------------------------------------------------P = 0.9582ですので、比例ハザード性が成立していないわけではないと解釈できます。
次は、年齢で補正した多変量Cox比例ハザードモデルを作成してみましょう。stcoxに年齢を追加します。
. stcox drug age
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -83.551879
Iteration 2: log likelihood = -83.324009
Iteration 3: log likelihood = -83.323546
Refining estimates:
Iteration 0: log likelihood = -83.323546
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(2) = 33.18
Log likelihood = -83.323546 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
drug | .1048772 .0477017 -4.96 0.000 .0430057 .2557622
age | 1.120325 .0417711 3.05 0.002 1.041375 1.20526
------------------------------------------------------------------------------年齢＋1歳あたりのハザード比は1.120325 [1.041375 – 1.20526]であり、P = 0.002です。比例ハザード性も確認してみましょう。
. estat phtest, detail
Test of proportional-hazards assumption
Time: Time
----------------------------------------------------------------
| rho chi2 df Prob>chi2
------------+---------------------------------------------------
drug | 0.00949 0.00 1 0.9603
age | -0.11758 0.42 1 0.5168
------------+---------------------------------------------------
global test | 0.43 2 0.8064
----------------------------------------------------------------いずれの変数もP > 0.05ですので、比例ハザードが成立していないわけではないということなります。
もし、P < 0.05となってしまって、比例ハザード性が成立していない可能性がでてきたときは、graphicalに評価してみましょう。大規模コホート研究では、検出力が高すぎて、しばしばP < 0.05になってしまいます。estat phtestコマンドに、plot(変数名)オプションを追加します。
estat phtest, plot(age) xsize(3.7) ysize(3.7)
lowessカーブが、時間に依存せずに横一直線になっていれば、比例ハザード性は成立していると考えます。対象数が少ないので、なんとも言えない感じですが、これくらいならばフラットと考えていいでしょう、多分。
このモデルでは、年齢1歳あたりハザード比が1.12 [1.04-1.20]だけ上昇することがわかりました。しかしながら、年齢1歳あたりでは幅が小さすぎて、臨床的感覚にマッチしません。そこで、年齢10歳あたりハザード比がどれくらい上がるかを算出してみます。そのためには、generateコマンドを使って、ageを10で割った変数age10を作成し、age10をモデルに投入します。
. generate age10 = age/10
. stcox drug age10
failure _d: died
analysis time _t: studytime
Iteration 0: log likelihood = -99.911448
Iteration 1: log likelihood = -83.551877
Iteration 2: log likelihood = -83.324008
Iteration 3: log likelihood = -83.323544
Refining estimates:
Iteration 0: log likelihood = -83.323544
Cox regression -- Breslow method for ties
No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(2) = 33.18
Log likelihood = -83.323544 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
drug | .1048772 .0477017 -4.96 0.000 .0430057 .2557622
age10 | 3.114865 1.161372 3.05 0.002 1.499927 6.46857
------------------------------------------------------------------------------年齢10歳あたりハザード比が3.11 [1.16-1.50]だけ上昇することがわかりました。かなりのリスク上昇ですね。
追跡期間が、連続変数ではなく、カテゴリ変数（1-2ヶ月、3-6ヶ月、7-12ヶ月等）である場合、通常のCox比例ハザードモデルではなく、discrete time hazard modelを用いた方がよい場合があります。
Discrete time hazard modelの概念については、下記の論文がよくまとまっています。
link関数をlogitにする場合とcomplementary log-logにする場合がありますが、Cox比例ハザードモデルに近いcomplementary log-log関数を利用した方がよさそうです。
2010年JAMAにdiscrete time hazard model with a complementary log-log linkを利用した論文が実際に発表されていました。
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)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 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.001Stataの特徴は比較的簡単に美しいグラフの作成できることです。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)
Y軸を、累積生存率ではなく、累積アウトカム発症率に変更し、durgカテゴリ別のKaplan-Meier曲線を描きます。
sts graph, xsize(3.7) ysize(3.7) by(drug) failure
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)
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)
グラフのタイトル”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)) 
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))
凡例を、生存曲線のプロット領域内（ring(0)）の右下（5時方向）に、縦１列（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")))
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)))
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))
テキストの配置（placement）には、sw以外にも下記を設定できます。テキストの表示位置に応じて適宜変更して下さい。
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))
大阪大学キャンパスライフ支援センターの山本陵平です。最近、いろいろなところで疫学研究の基本を教える機会が増えてきたので、ホームページにまとめて公開することにしました。統計、研究デザイン、バイアス、プログラミング等の疫学研究を行う上で役に立ちそうな内容を、書き連ねていこうと思います。