Меня интересует вычисление площади под кривой (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-статистику, учитывая эту информацию.

Я бы рекомендовал статью Hanley & McNeil 1982 года «Значение и использование области под кривой рабочих характеристик приемника (ROC) ».

пример

Они имеют следующую таблицу статуса заболевания и результата теста (например, в соответствии с предполагаемым риском из логистической модели). Первое число справа - это количество пациентов с истинным статусом заболевания «нормальный», а второе число - количество пациентов с истинным статусом заболевания «аномальный»:

(1) Определенно нормально: 33/3
(2) Вероятно нормально: 6/2
(3) Сомнительно: 6/2
(4) Вероятно ненормально: 11/11
(5) Определенно ненормально: 2/33

Таким образом, в общей сложности 58 «нормальных» пациентов и «51» ненормальных. Мы видим, что когда предиктор равен 1, «Определенно нормальный», пациент обычно нормальный (верно для 33 из 36 пациентов), а когда он равен 5, «Определенно ненормальный», пациент обычно ненормальный (истинный для 33 из 35 пациентов), поэтому предиктор имеет смысл. Но как мы должны судить пациента с оценкой 2, 3 или 4? То, что мы устанавливаем в качестве предела для оценки пациентов как ненормальных или нормальных, определяет чувствительность и специфичность полученного теста.

Чувствительность и специфичность

Мы можем рассчитать расчетную чувствительность и специфичность для разных срезов. (Теперь я просто напишу «чувствительность» и «специфичность», чтобы приблизительная природа значений была неявной.)

Если мы выберем нашу отсечку так, чтобы мы классифицировали всех пациентов как ненормальных, независимо от того, что говорят результаты их тестов (то есть мы выбрали отсечку 1+), мы получим чувствительность 51/51 = 1. Специфичность будет 0 / 58 = 0. Звучит не очень хорошо.

Хорошо, давайте выберем менее строгую отсечку. Мы классифицируем пациентов как ненормальных, только если у них результат теста 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 вручную

Мы можем очень легко вычислить площадь под кривой 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))

Ответ снова 0,8931711, площадь под кривой ROC. Это всегда будет так.

Графическое представление соответствия

Как указал Харрелл в своем ответе, это также имеет графическую интерпретацию. Давайте нарисуем тестовую оценку (оценку риска) по оси 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))

Да, это 0.8931711, площадь под кривой ROC.

Более простые способы расчета AUC (в R)

Но давайте сделаем жизнь проще для себя. Существуют различные пакеты, которые автоматически рассчитывают AUC для нас.

Пакет Эпи

EpiПакет создает хороший ROC кривой с различными статистическими данными (включая АУК) встроенные:

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

Пакет PROC

Мне также нравится pROCпакет, так как он может сгладить оценку ROC (и вычислить оценку AUC на основе сглаженной ROC):

(Красная линия - это исходная ROC, а черная линия - сглаженная ROC. Также обратите внимание на стандартное соотношение сторон 1: 1. Это имеет смысл использовать, поскольку чувствительность и специфичность имеют диапазон 0–1.)

Расчетное значение AUC для сглаженного ROC составляет 0,9107, аналогично, но немного больше, чем AUC для сглаженного ROC (если вы посмотрите на рисунок, вы легко поймете, почему он больше). (Хотя у нас действительно слишком мало возможных различных значений результатов теста, чтобы рассчитать плавный AUC).

Пакет RMS

rmsПакет Harrell может вычислять различные связанные статистические данные соответствия, используя 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

Заключительные слова

Многие люди думают, что AUC говорит нам, насколько «хорош» тест. И некоторые люди думают, что AUC - это вероятность того, что тест правильно классифицирует пациента. Это не так . Как видно из приведенного выше примера и расчетов, AUC говорит нам кое-что о семействе тестов, по одному тесту для каждого возможного отсечения.

А AUC рассчитывается на основе отсечек, которые никогда бы не использовались на практике. Почему мы должны заботиться о чувствительности и специфичности «бессмысленных» предельных значений? Тем не менее, именно на этом (частично) основывается AUC. (Конечно, если AUC очень близко к 1, почти каждый возможный тест будет иметь большую дискриминационную силу, и мы все были бы очень рады.)

Интерпретация AUC «случайная нормальная аномальная» пара хороша (и может быть распространена, например, на модели выживания, где мы видим, является ли человек с самой высокой (относительной) опасностью, умирающим самым ранним). Но никто бы никогда не использовал это на практике. Это редкий случай, когда кто-то знает, что у него один здоровый и один больной человек, не знает, какой человек болен, и должен решить, кого из них лечить. (В любом случае, решение легко; отнеситесь к тому, у кого самый высокий предполагаемый риск.)

Поэтому я думаю, что изучение фактической кривой ROC будет более полезным, чем просто рассмотрение сводного показателя AUC. И если вы используете ROC вместе с (оценками) издержек ложных срабатываний и ложных отрицаний, а также с базовыми показателями того, что вы изучаете, вы можете получить что-нибудь.

Также обратите внимание, что AUC измеряет только дискриминацию , а не калибровку. Таким образом, он измеряет, можете ли вы различать двух человек (один болен и один здоров), основываясь на оценке риска. Для этого, он смотрит только на относительных значений риска (или ряды, если вы будете, сравни интерпретации тест Вилкоксона-Манна-Уитни), а не абсолютные, которые вы должны быть заинтересованы. Например, если вы разделите каждый риск По оценкам вашей логистической модели на 2, вы получите точно такой же AUC (и ROC).

При оценке модели риска калибровка также очень важна. Чтобы проверить это, вы посмотрите на всех пациентов с оценкой риска около, например, 0,7, и посмотрите, действительно ли примерно 70% из них были больны. Делайте это для каждой возможной оценки риска (возможно, с использованием некоторого сглаживания / локальной регрессии). Составьте график результатов, и вы получите графическую меру калибровки .

Если есть модель с обеими хорошей калибровкой и хорошей дискриминацией, то вы начинаете иметь хорошую модель. :)

Посмотрите на этот вопрос: Понимание кривой ROC

Вот как построить кривую ROC (из этого вопроса):

Рисование ROC кривой

учитывая набор данных, обработанный вашим классификатором рейтинга

• ранговые тестовые примеры по убыванию баллов
• $(0, 0)$
• $x$
  • $x$$1/\text{pos}$
  • $x$$1/\text{neg}$

$\text{pos}$$\text{neg}$

Вы можете использовать эту идею для вычисления AOC ROC вручную, используя следующий алгоритм:

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

Эта хорошая анимированная картинка должна проиллюстрировать этот процесс

В сообщении Карла много отличной информации. Но я еще не видел за последние 20 лет пример кривой ROC, которая изменила бы чье-либо мышление в правильном направлении. Единственное значение кривой ROC, по моему скромному мнению, состоит в том, что ее площадь равна очень полезной вероятности согласования. Сама кривая ROC заставляет читателя использовать срезы, что является плохой статистической практикой.

$c$$Y=0,1$$x$$Y=1$$y$$Y=0$$Y=1$

$n$

Для функции Hmiscпакета R rcorr.censвыведите весь результат, чтобы увидеть больше информации, особенно стандартную ошибку.

Вот альтернатива естественному способу вычисления 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{Sum of true ranks}-0.5PN(PN+1)}{PN(N-PN)}$$