ループとROC曲線作成の練習①前編 | 医療統計とStataプログラミングの部屋

このブログでは、統計解析ソフトStataのプログラミングのTipsや便利コマンドを紹介しています．

Facebook groupでは、ちょっとした疑問や気づいたことなどを共有して貰うフォーラムになっています．ブログと合わせて個人の学習に役立てて貰えれば幸いです．

さて、今回はROC（Receiver operating characteristic）曲線とAUC（Area under curve）を求める方法については以前概説しました．

今回は、いくつかのサブグループに分けて解析した結果をまとめる作業を通じて、

ループコマンドに強くなる
cutptコマンドのおさらい
“addplot” オプションに強くなる

事を目標としたいと思います．作っていたら結構長くなってしまったので前編と後編に分割しますね．今回は前編です．

1．サンプルデータのダウンロード
2．グラフを描いてみる

1．サンプルデータのダウンロード

ダミーデータを置いておきます．これを元にして解析をしてみましょう．

sampledata_cutpt ダウンロード

2．グラフを描いてみる

import delimited "C:\...\sampledata_cutpt.csv"
*** 変数設定
recode age (min/55=1) (55/65=2) (65/max=3), gen(agecat)
label define agecat 1 "Age < 55" 2 "Age 55 to 65" 3 "Age >= 65"
label values agecat agecat
*** Outcomeを従属変数、marker_a, marker_bの性能を調べる．まずはmarker_aから
roctab outcome marker_a, graph

これでROC曲線とAUCが得られますが、殺風景なのでいろいろと加えていってみましょう．

まずは前回の復習ということで、cutptを使ってYouden indexから求めたカットオフをグラフ上に表示させます．マーカーの色は赤く目立たせてみます．

**1 addplot
でグラフの一部を目立たせる
cutpt outcome marker_a, noadjust youden   
roctab outcome marker_a, graph addplot(scatteri `e(sens)' `=1 - e(spec)', mcolor(red))

legendの位置がイマイチなのでグラフエリアに縦に格納してみます．legend内はlabel(# “labname”) label(# “labname”) と並べていきます．

**2 legend
を整える
cutpt outcome marker_a, noadjust youden   
roctab outcome marker_a, graph addplot(scatteri `e(sens)' `=1 - e(spec)', mcolor(red)) legend(label(1 "ROC") label(3 "Cutpoint") col(1) ring(0) pos(5))

グラフエリアを白背景にしつつ、グリッドを消すのですが、twowayと異なり、グリッドを消すときはそれぞれの軸ラベル内のオプションで処理することが必要です．

**3 graph areaの設定
cutpt outcome marker_a, noadjust youden   
roctab outcome marker_a, graph addplot(scatteri `e(sens)' `=1 - e(spec)', mcolor(red)) legend(label(1 "ROC") label(3 "Cutpoint") col(1) ring(0) pos(5)) graphregion(color(white) lcolor(white)) xlabel(, nogrid) ylabel(, nogrid)

さて、ここまで行くと感度特異度をグラフ中に埋め込みたくなりますね！ということで

**4 感度・特異度テキストで表示
cutpt outcome marker_a, noadjust youden   
roctab outcome marker_a, graph addplot(scatteri `e(sens)' `=1 - e(spec)', mcolor(red)) legend(label(1 "ROC") label(3 "Cutpoint") col(1) ring(0) pos(5)) graphregion(color(white) lcolor(white)) xlabel(, nogrid) ylabel(, nogrid) text(0.8 0.4 "Sensitivity: `e(sens)'" "Specificity: `e(spec)'")

おや、これはどうしたことでしょうか．桁数が増えてしまっています．こんな時は少し前の記事を参照してもらって、local macroに桁数を絞って格納してしまいましょう．

**5 感度・特異度の桁数調整
cutpt outcome marker_a, noadjust youden   
  local sens: dis %3.2f `e(sens)'
  local spec: dis %3.2f `e(spec)'
roctab outcome marker_a, graph addplot(scatteri `e(sens)' `=1 - e(spec)', mcolor(red)) legend(label(1 "ROC") label(3 "Cutpoint") col(1) ring(0) pos(5)) graphregion(color(white) lcolor(white)) xlabel(, nogrid) ylabel(, nogrid) text(0.8 0.4 "Sensitivity: `sens'" "Specificity: `spec'")