ROC 曲線を用いたカットオフ値の計算
SAS では LOGISTIC Procedure を用いて ROC 曲線の描画 (proc logistic plots = roc(id = prob) など) が簡単にできますが、カットオフ値については自前で計算する必要がありそうででした。 そこで、カットオフ値を求めるためのいくつかの指標を計算するマクロを作成しました。
カットオフ値の指標としては、Youden's index [1] (感度 + 特異度 – 1; 大きい方が良い) や、縦軸に感度を横軸に 1 – 特異度を取った ROC 曲線の図の左上からの距離 (小さい方が良い) などがあります。 Youden's index は推定法などに関する論文が多数ありますし、解釈も簡単なので、どちらかというと良く用いられているようです。 しかし、診断法のデザインと解析は奥深いものですし、様々な注意点がありますので、実際に利用する際には教科書など (例えば [2] など) をしっかりと勉強した方が良いように思います。
見ての通り Youden's index は感度と特異度の和に基づく指標です。 感度と特異度の和が高ければ、診断性能は悪くなさそうに見えますが、必要な感度や特異度の大きさは研究の目的によって変わってくるはずです。 したがって、機械的に Youden's index の高いものを最適カットオフだ、というのは少し乱暴だと思われます。 個人的には最適カットオフ値を計算に基づいて決めるという考え方ではなく、必要な感度や特異度を満たすようなカットオフ値がどのあたりになりそうなのかあたりを付けるために使うという見方の方がしっくりくるかなと考えています。
注意点 SASの出力の仕様上、感度や特異度を含んだROC曲線のデータと解析に使ったデータを結合する際に、ロジスティック回帰モデルを当てはめて得られた確率の予測値をキーとしてマージを行う必要があります。 しかし、SASの出力の仕様により上記の二つのデータ間に微妙な数値計算誤差の違いが出ることがあるため、マージがうまくいかない場合があります。 本マクロでは、確率予測値を 1e-12 で丸めてからマージをしますが、うまくいかない場合があるかもしれません。 出力データセットを確認してマージがうまくいっているかどうかを確認してからご利用下さい。
ページ下に作成したマクロ (%rocCut) を記載しました。このマクロを読み込んで、変数名などを適切に指定して実行すると、各指標が出力されます。
例
ここでは、疑似乱数を用いて適当なデータを生成し、%rocCut を用いて各指標を求めてみます。 このデータでは、二値結果変数 y に対して変数 x のよいカットオフ値を求めたいとします。
%macro rocCut(data, out, y, var, cov =, eventcode = 1, alpha = 0.05);
ods listing close;
ods results off;
proc logistic data = &data. ;
model &y.(event="&eventcode.") = &var. &cov. / outroc = _rocdat rocci roceps = 1e-12;
output out = _outp p = _PROB_;
ods output ROCAssociation = _auc(keep = Area LowerArea UpperArea);
run;
ods results on;
ods listing;
proc sort data = _outp; by _PROB_;
proc sort data = _rocdat; by _PROB_;
data _outp; set _outp; _PROB_ = round(_PROB_, 1e-12);
data _rocdat; set _rocdat; _PROB_ = round(_PROB_, 1e-12);
data _rocdat; merge _outp(where = (_PROB_ ^= .)) _rocdat; by _PROB_;
data &out.;
set _rocdat;
if _N_ = 1 then set _auc(obs = 1);
se = _SENSIT_;
_za_ = quantile("Normal", 1 - &alpha./2);
_nn_ = _POS_ + _FALNEG_; _x_ = _POS_; _p_ = _x_ / _nn_;
selcl = (_x_ + _za_**2*0.5)/(_nn_ + _za_**2) -
_za_*sqrt(_nn_)/(_nn_ + _za_**2)*sqrt(_p_*(1 - _p_) + _za_**2/_nn_*0.25);
seucl = (_x_ + _za_**2*0.5)/(_nn_ + _za_**2) +
_za_*sqrt(_nn_)/(_nn_ + _za_**2)*sqrt(_p_*(1 - _p_) + _za_**2/_nn_*0.25);
_SPEC_ = 1 - _1MSPEC_;
_nn_ = _NEG_ + _FALPOS_; _x_ = _NEG_; _p_ = _x_ / _nn_;
splcl = (_x_ + _za_**2*0.5)/(_nn_ + _za_**2) -
_za_*sqrt(_nn_)/(_nn_ + _za_**2)*sqrt(_p_*(1 - _p_) + _za_**2/_nn_*0.25);
spucl = (_x_ + _za_**2*0.5)/(_nn_ + _za_**2) +
_za_*sqrt(_nn_)/(_nn_ + _za_**2)*sqrt(_p_*(1 - _p_) + _za_**2/_nn_*0.25);
_nn_ = _POS_ + _FALPOS_ + _NEG_ + _FALNEG_; _x_ = _POS_ + _NEG_; _p_ = _x_ / _nn_;
acc = (_POS_ + _NEG_) / _nn_;
acclcl = (_x_ + _za_**2*0.5)/(_nn_ + _za_**2) -
_za_*sqrt(_nn_)/(_nn_ + _za_**2)*sqrt(_p_*(1 - _p_) + _za_**2/_nn_*0.25);
accucl = (_x_ + _za_**2*0.5)/(_nn_ + _za_**2) +
_za_*sqrt(_nn_)/(_nn_ + _za_**2)*sqrt(_p_*(1 - _p_) + _za_**2/_nn_*0.25);
youden = _SENSIT_ + _SPEC_ - 1;
distance = sqrt((1 - _SENSIT_)**2 + _1MSPEC_**2);
drop _SENSIT_ _1MSPEC_ _p_ _za_ _nn_ _x_ _SOURCE_;
rename
_SPEC_ = sp _POS_ = pos _NEG_ = neg
_FALPOS_ = fpos _FALNEG_ = fneg Area = auc
LowerArea = auclcl UpperArea = aucucl _LEVEL_ = level
_PROB_ = prob;
proc sort data = &out.; by descending youden distance;
proc datasets lib = work;
delete _rocdat _outp _auc;
run; quit;
%mend rocCut;
data sample;
call streaminit(130301);
beta = 2;
do i = 1 to 100;
x = rand('Normal', 0, 2);
logit = beta*x;
theta = 1 / (1 + exp(-logit));
y = rand('Bernoulli', theta);
output;
end;
drop beta logit theta;
run;
* %rocCut(入力データセット名, 出力データセット名, 結果変数, 連続値変数);
%rocCut(sample, outroc, y, x);
proc print data = outroc(obs = 5);
run;
実行結果です。
i x y level prob pos neg fpos fneg auc auclcl aucucl se
27 0.26754 1 1 0.54606 37 51 5 7 0.9456 0.9057 0.9855 0.84091
57 0.18006 0 1 0.50537 37 50 6 7 0.9456 0.9057 0.9855 0.84091
97 0.28887 1 1 0.55590 36 51 5 8 0.9456 0.9057 0.9855 0.81818
92 0.06122 1 1 0.45010 39 47 9 5 0.9456 0.9057 0.9855 0.88636
5 0.11846 1 1 0.47665 38 48 8 6 0.9456 0.9057 0.9855 0.86364
selcl seucl sp splcl spucl acc acclcl accucl youden distance
0.70634 0.92073 0.91071 0.80744 0.96126 0.88 0.80188 0.93001 0.75162 0.18243
0.70634 0.92073 0.89286 0.78532 0.94996 0.87 0.79020 0.92243 0.73377 0.19181
0.68039 0.90487 0.91071 0.80744 0.96126 0.87 0.79020 0.92243 0.72890 0.20256
0.76021 0.95047 0.83929 0.72194 0.91307 0.86 0.77863 0.91474 0.72565 0.19683
0.73291 0.93597 0.85714 0.74264 0.92579 0.86 0.77863 0.91474 0.72078 0.19749
pos は陽性的中数、neg は陰性的中数、fpos は偽陽性数、fneg は偽陰性数、auc は AUC/C-index (auclcl, aucucl は 95% CLs)、se は感度 (selcl, seucl は 95% CLs)、sp は特異度 (splcl, spucl は 95% CLs)、acc は正診率 (acclcl, accucl は 95% CLs)、youden が Youden's index [1]、distance が左上からの距離による指標です。
この結果では、x = 0.26754 をカットオフ値とすると、感度が 84.1 %、特異度が 91.1 % となりました。 このとき、Youden's index = 0.75162 となり、distance = 0.18243、どちらの指標を用いても同じカットオフ値が得られています (2 行目以降を見ると分かりますが、かならず同じカットオフ値が得られるわけではありません)。
参考文献
- Youden WJ. Index for rating diagnostic tests. Cancer 1950;3:32–35.
- Zhou X-H, McClish DK, Obuchowski NA. Statistical methods in diagnostic medicine. New York: John Wiley, 2002.
履歴
- 2016/12/21 プログラムのバグ修正 (Thanks to Dr. Kyoko Nomura)
- 2015/07/03 コメント追加
- 2015/01/22 プログラムのバグ修正
- 2014/06/04 プログラムのバグ修正
- 2013/03/01 公開