Iバイナリロジスティック回帰モデルの曲線下面積(AUC)またはc統計量を手動で計算することに関心があります。

たとえば、検証データセットには、従属変数の真の値、保持(1 =保持、0 =保持なし)、およびトレーニングを使用して構築されたモデルを使用した回帰分析によって生成された各観測値の予測保持ステータスがあります。セット(これは0から1の範囲になります)。

私の最初の考えは、モデル分類の「正しい」数を特定し、「正しい」観測数を単純に合計観測数で割って計算することでした。 c-統計。 「正しい」とは、観測値の真の保持ステータス= 1で、予測された保持ステータスが> 0.5の場合、それは「正しい」分類です。さらに、観測値の真の保持ステータス= 0で、予測された保持ステータスが< 0.5の場合、これも「正しい」分類です。予測値= 0.5のときに「同点」が発生すると想定しますが、その現象は検証データセットでは発生しません。一方、「誤った」分類は、観測値の真の保持ステータス= 1で、予測された保持ステータスが< 0.5の場合、または結果の真の保持ステータス= 0であり、予測される保持ステータスは> 0.5です。 TP、FP、FN、TNについては知っていますが、この情報からc統計量を計算する方法については知りません。

回答

ハンリーの&マクニールの1982年の論文 受信者動作特性(ROC)の下の領域の意味と使用法をお勧めします)曲線

次の病状と検査結果の表があります(たとえば、ロジスティックモデルからの推定リスクに対応します)。右側の最初の数字は、 true の病状が「正常」である患者の数です。2番目の数字は true の病状が「異常」である患者の数です。

(1)間違いなく正常:33/3
(2)おそらく正常:6/2
(3)疑わしい:6/2
(4)おそらく異常: 11/11
(5)明らかに異常:2/33

つまり、合計58人の「正常な」患者と「51人の」異常な患者がいます。予測子が1の場合、「完全に正常」の場合、患者は通常正常であり(36人の患者のうち33人に当てはまります)、5の場合、「完全に異常」の場合、患者は通常異常です(33人の患者に当てはまります)。 35人の患者)なので、予測子は理にかなっています。しかし、スコアが2、3、または4の患者をどのように判断する必要がありますか?患者を異常または正常と判断するためのカットオフを設定して、結果のテストの感度と特異度を決定します。

感度と特異度

推定を計算できます。 / em>さまざまなカットオフの感度と特異性。 (これからは「感度」と「特異度」を記述し、値の推定された性質を暗黙的にします。)

カットオフを選択してすべてを分類する場合患者が異常であるとテスト結果が何を示していても(つまり、カットオフ1+を選択)、感度は51/51 = 1になります。特異度は0/58 = 0になります。そうではありません。とても良い音です。

OK、それではそれほど厳密ではないカットオフを選択しましょう。検査結果が2以上の場合のみ異常と分類します。次に、3人の異常な患者を見逃し、感度は48/51 = 0.94になります。しかし、33/58 = 0.57と、はるかに高い特異性があります。

これを継続して、さまざまなカットオフ(3、4、5、> 5)を選択できます。 (最後のケースでは、テストスコアが5である場合でも、 患者を異常として分類しません。)

ROC曲線

すべての可能なカットオフに対してこれを行い、1から特異度を引いたものに対する感度をプロットすると、ROC曲線が得られます。次のRコードを使用できます。

 # Data norm = rep(1:5, times=c(33,6,6,11,2)) abnorm = rep(1:5, times=c(3,2,2,11,33)) testres = c(abnorm,norm) truestat = c(rep(1,length(abnorm)), rep(0,length(norm))) # Summary table (Table I in the paper) ( tab=as.matrix(table(truestat, testres)) )  

出力は次のとおりです。

  testres truestat 1 2 3 4 5 0 33 6 6 11 2 1 3 2 2 11 33  

さまざまな統計を計算できます:

 ( tot=colSums(tab) ) # Number of patients w/ each test result ( truepos=unname(rev(cumsum(rev(tab[2,])))) ) # Number of true positives ( falsepos=unname(rev(cumsum(rev(tab[1,])))) ) # Number of false positives ( totpos=sum(tab[2,]) ) # The total number of positives (one number) ( totneg=sum(tab[1,]) ) # The total number of negatives (one number) (sens=truepos/totpos) # Sensitivity (fraction true positives) (omspec=falsepos/totneg) # 1 − specificity (false positives) sens=c(sens,0); omspec=c(omspec,0) # Numbers when we classify all as normal  

これを使用して、(推定)ROC曲線をプロットできます。

 plot(omspec, sens, type="b", xlim=c(0,1), ylim=c(0,1), lwd=2, xlab="1 − specificity", ylab="Sensitivity") # perhaps with xaxs="i" grid() abline(0,1, col="red", lty=2)  

AUC曲線

手動で計算AUC

台形の面積の式を使用して、ROC曲線の下の面積を非常に簡単に計算できます。

 height = (sens[-1]+sens[-length(sens)])/2 width = -diff(omspec) # = diff(rev(omspec)) sum(height*width)  

結果は0.8931711です。

一致度測定値

AUCは一致度測定値と見なすこともできます。一方が正常でもう一方が異常である可能性のあるすべてのペアの患者を取得すると、最も高い(最も「異常に見える」)テスト結果を持つのが異常な患者である頻度を計算できます(それらは同じ値を持っているので、これを「半分の勝利」と見なします):

 o = outer(abnorm, norm, "-") mean((o>0) + .5*(o==0))  

答えは、ROC曲線の下の領域である0.8931711です。これは常に当てはまります。

一致のグラフィカルビュー

ハレルが回答で指摘したように、これにもグラフィカルな解釈があります。 y 軸にテストスコア(リスク推定値)をプロットし、 x 軸に実際の病状をプロットしてみましょう(ここでは、重複するポイントを示すために、多少のジッターがあります):

 plot(jitter(truestat,.2), jitter(testres,.8), las=1, xlab="True disease status", ylab="Test score")  

真の疾患に対するリスクスコアの散布図ステータス。

次に、左側の各ポイント(「正常な」患者)と右側の各ポイント(「異常な」患者)の間に線を引きます。正の傾きを持つ線の割合(つまり、一致ペアの割合)は一致指数です(平らな線は「50%一致」としてカウントされます)。

同点の数(リスクスコアが等しい)のため、この例の実際の線を視覚化するのは少し難しいですが、ある程度のジッターと透明性があれば、妥当なプロットを得ることができます:

 d = cbind(x_norm=0, x_abnorm=1, expand.grid(y_norm=norm, y_abnorm=abnorm)) library(ggplot2) ggplot(d, aes(x=x_norm, xend=x_abnorm, y=y_norm, yend=y_abnorm)) + geom_segment(colour="#ff000006", position=position_jitter(width=0, height=.1)) + xlab("True disease status") + ylab("Test\nscore") + theme_light() + theme(axis.title.y=element_text(angle=0))  

真の病状に対するリスクスコアの散布図、可能なすべての観測ペア間の線。

ほとんどの線が上向きに傾斜しているため、一致指数が高くなります。また、各タイプの観測ペアからのインデックスへの寄与も確認できます。そのほとんどは、リスクスコアが1の正常な患者とリスクスコアが5の異常な患者(1〜5ペア)からのものですが、1〜4ペアおよび4〜5ペアからもかなり多くのものがあります。また、勾配の定義に基づいて実際の一致指数を計算するのは非常に簡単です。

 d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm)) mean((d$slope > 0) + .5*(d$slope==0))  

答えは再び0.8931711、つまりAUCです。

ウィルコクソン-マン-ホイットニー検定

一致度とウィルコクソン-マン-ホイットニーの間には密接な関係があります。テスト。実際、後者は、一致の確率(つまり、最も「異常に見える」テスト結果を持つランダムの正常と異常のペアの異常な患者)が正確に0.5であるかどうかをテストします。そして、その検定統計量は、推定された一致確率の単純な変換です。

 > ( wi = wilcox.test(abnorm,norm) ) Wilcoxon rank sum test with continuity correction data: abnorm and norm W = 2642, p-value = 1.944e-13 alternative hypothesis: true location shift is not equal to 0  

検定統計量(W = 2642)は、一致するペアの数をカウントします。可能なペアの数で割ると、おなじみの数になります。

 w = wi$statistic w/(length(abnorm)*length(norm))  

はい、ROC曲線の下の面積は0.8931711です。

AUCを計算する簡単な方法(R)

しかし、私たち自身の生活を楽にしてみましょう。 AUCを自動的に計算するさまざまなパッケージがあります。

Epiパッケージ

Epiパッケージは、さまざまなROC曲線を作成します。埋め込まれた統計(AUCを含む):

 library(Epi) ROC(testres, truestat) # also try adding plot="sp"  

EpiパッケージのROC曲線

pROCパッケージ

pROCパッケージも気に入っています。 ROC推定値を平滑化する(および平滑化されたROCに基づいてAUC推定値を計算する):

pROCパッケージからのROC曲線(平滑化および平滑化)

(赤い線は元のROCで、黒い線は平滑化されたROCです。デフォルトの1:1のアスペクト比にも注意してください。感度と特異性の両方が0–1であるため、これを使用するのが理にかなっています。範囲。)

平滑化 ROCからの推定AUCは0.9107であり、平滑化されていないROCからのAUCと似ていますが、わずかに大きくなっています(図を見ると、なぜ大きいのかが簡単にわかります)。 (ただし、実際には、スムーズなAUCを計算するには、可能な個別のテスト結果値が少なすぎます。)

rmsパッケージ

Harrellのrmsパッケージrcorr.cens()関数を使用して、関連するさまざまな一致統計を計算できます。出力のC IndexはAUCです:

 > library(rms) > rcorr.cens(testres,truestat)[1] C Index 0.8931711  

caToolsパッケージ

最後に、caToolsパッケージとそのcolAUC()関数があります。他のパッケージに比べていくつかの利点があり(主に速度と多次元データを処理する機能– ?colAUCを参照)、時々役立つことがあります。しかしもちろん、何度も計算したのと同じ答えが得られます:

 library(caTools) colAUC(testres, truestat, plotROC=TRUE) [,1] 0 vs. 1 0.8931711  

caToolsパッケージのROC曲線

最後の言葉

多くの人は、AUCが「良い」方法を教えてくれると考えているようです。テストはです。また、AUCは、テストによって患者が正しく分類される確率であると考える人もいます。 ではありません。上記の例と計算からわかるように、AUCはファミリのテストについて何かを教えてくれます。可能なカットオフごとに1つのテストです。

そして、AUCはに基づいて計算されます。実際には決して使用しないカットオフ。 「無意味な」カットオフ値の感度と特異性を気にする必要があるのはなぜですか?それでも、それはAUCが(部分的に)基づいているものです。 (もちろん、AUCが非常に 1に近い場合、ほぼすべての可能なテストに大きな識別力があり、私たちは皆非常に満足しています。)

ランダム正規分布–AUCの異常なペアの解釈は優れています(たとえば、最も早く死亡する(相対的な)ハザードが最も高い人であるかどうかを確認する生存モデルに拡張できます)。しかし、実際には決して使用しません。健康な人と病気の人がいることを知っている 、病気の人が誰であるかわからない、そしてしなければならないというのはまれなケースです。それらのどれを扱うかを決定します。 (いずれの場合も、決定は簡単です。推定リスクが最も高いものを処理してください。)

したがって、実際の ROC曲線を調べる方が、単に見るよりも役立つと思います。 AUC要約測定。また、ROCを、誤検知と誤検知の(推定)コストと、学習しているものの基本レートとともに使用すると、どこかに到達できます。

また、AUCは識別のみを測定し、キャリブレーションは測定しないことに注意してください。つまり、リスクスコアに基づいて、2人(病気の人と健康な人)を区別できるかどうかを測定します。このため、相対リスク値(または、必要に応じてランク、ウィルコクソン-マン-ホイットニー検定の解釈を参照)のみを調べ、絶対値は調べません。 em>興味があります。たとえば、ロジスティックモデルの各リスク推定値を2で割ると、まったく同じAUC(およびROC)が得られます。

リスクモデルを評価する場合、キャリブレーションも非常に重要です。これを調べるために、リスクスコアが約0.7(0.7など)のすべての患者を調べ、これらの約70%が実際に病気であるかどうかを確認します。考えられるリスクスコアごとにこれを実行します(おそらく、ある種の平滑化/局所回帰を使用します)。結果をプロットすると、キャリブレーションのグラフィカルな測定値が得られます。

優れたキャリブレーションと優れた識別の両方を備えたモデルがある場合は、 良いモデルを持ち始めます。 🙂

コメント

  • ありがとう、@ Karl Ove Hufthammer、これは私が今まで受け取った中で最も徹底的な答えです。特に、”最後の言葉”セクションに感謝します。素晴らしい仕事です!もう一度ありがとう!
  • この詳細な回答をありがとうございました。私はEpi :: ROC()v2.2.6がAUCが1.62であると確信しているデータセットを使用しています(メンタリストの研究ではありません)が、ROCによると、上記のコードの結果は0.56であると信じていますin。
  • sens=c(sens,0); omspec=c(omspec,0)に小さなエラーがあると思います。’これはsens=c(0, sens); omspec=c(0, omspec)?先頭の0で正しくプロットされますが、現在の回答の方法ではありません。
  • いいえ、現在の定義は、AFAICS、正しい、@ steveb、正しいプロットになります。おそらく紛らわしいのは、ROC曲線が左から左へではなく、右から左へ(つまり、右上隅から左下隅へ)に描かれていることだと思います。右、ほとんどのプロットがそうであるように。これは、変数を定義した結果です。左から右にプロットすることもできます(プロットする前にsensomspecの両方を逆にすることによって)

回答

次の質問をご覧ください: ROC曲線を理解する

ROC曲線を作成する方法は次のとおりです(その質問から):

ROC曲線の描画

ランク付け分類子

  • スコアの減少に関するテスト例のランク付け
  • $(0、0)$から開始
  • 各例$ x $(降順)
    • $ x $が正の場合、$ 1 / \ text {pos} $を上に移動します
    • $ x $が負の場合、$ 1 / \ text {neg} $を右に移動します

ここで、$ \ text {pos} $と$ \ text {neg} $は、それぞれ正と負の例の割合です。

このアイデアを使用して、次のアルゴリズムを使用してAUCROCを手動で計算できます。

auc = 0.0 height = 0.0 for each training example x_i, y_i if y_i = 1.0: height = height + tpr else auc = auc + height * fpr return auc 

この素敵なgifアニメーションの画像は、これを示しているはずです。より明確なプロセス

曲線の作成

コメント

  • @Alexeyに感謝グリゴレフ、これは素晴らしいビジュアルであり、将来的に役立つ可能性があります! +1
  • “正と負の例の割合”について少し説明してください。 2軸の最小単位値?
  • @Allan Ruin:posここで、正のデータの数を意味します。 20個のデータポイントがあり、そのうち11個のポイントが1であるとします。したがって、グラフを描画すると、11×9(高さx幅)の長方形ができます。 Alexey Grigorevはスケーリングを行いましたが、必要に応じて’のままにします。さて、各ステップでチャート上で1を移動するだけです。

回答

Karlの投稿にはたくさんありますしかし、過去20年間、誰かの考え方を良い方向に変えたROC曲線の例はまだ見ていません。私の謙虚な意見におけるROC曲線の唯一の価値は、その面積がたまたま非常に有用な一致確率に等しいということです。 ROC曲線自体は、読者にカットオフを使用するように誘惑しますが、これは悪い統計手法です。

$ c $インデックスを手動で計算する限り、$ Y = 0,1 $でプロットを作成します。 x $軸と$ y $軸上の$ Y = 1 $の連続予測子または予測確率。 $ Y = 0 $のすべてのポイントを$ Y = 1 $のすべてのポイントに接続する場合、正の勾配を持つ線の比率が一致確率になります。

分母がこの設定の$ n $は不適切な精度のスコアリング規則であるため、回避する必要があります。これには、正しく分類された比率、感度、および特異性が含まれます。

R Hmiscパッケージrcorr.cens関数の場合、結果全体で、より多くの情報、特に標準エラーを確認できます。

コメント

  • ありがとう、@ Frank Harell、あなたの視点に感謝します。 ‘カットオフが好きではないので、一致確率としてc統計量を使用します。もう一度ありがとう!

回答

これは、台形を使用するだけでAUCを計算する自然な方法の代替手段です。 ROC曲線の下の領域を取得するルール。

AUCは、ランダムにサンプリングされた正の観測値が、ランダムにサンプリングされた負の観測値よりも大きい(正の)予測確率を持つ確率に等しくなります。これを使用して、正と負の観測値のすべてのペアワイズの組み合わせを調べることにより、任意のプログラミング言語でAUCを非常に簡単に計算できます。サンプルサイズが大きすぎる場合は、観測値をランダムにサンプリングすることもできます。ペンと紙を使用してAUCを計算する場合、サンプルが非常に少ない/時間が長い場合を除いて、これは最善のアプローチではない可能性があります。たとえば、Rの場合:

n <- 100L x1 <- rnorm(n, 2.0, 0.5) x2 <- rnorm(n, -1.0, 2) y <- rbinom(n, 1L, plogis(-0.4 + 0.5 * x1 + 0.1 * x2)) mod <- glm(y ~ x1 + x2, "binomial") probs <- predict(mod, type = "response") combinations <- expand.grid(positiveProbs = probs[y == 1L], negativeProbs = probs[y == 0L]) mean(combinations$positiveProbs > combinations$negativeProbs) [1] 0.628723 

pROCパッケージを使用して確認できます:

library(pROC) auc(y, probs) Area under the curve: 0.6287 

ランダムサンプリングの使用:

mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE)) [1] 0.62896 

回答

  1. あなたには観察の真の価値があります。
  2. 事後確率を計算し、この確率で観測値をランク付けします。
  3. $ P $のカットオフ確率と観測数$ N $を想定:
    $$ \ frac {\ text {真のランクの合計} -0.5PN(PN + 1)} { PN(N-PN)} $$

コメント

  • @ user73455 … 1)はい、私には真の価値があります観察のため。 2)事後確率は、各観測値の予測確率と同義ですか? 3)理解した;ただし、”真のランクの合計”とは何ですか?また、この値をどのように計算しますか?おそらく、例はこの答えをより完全に説明するのに役立つでしょうか?ありがとうございます!

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です