Решающий порог для модели логистической регрессии glm в R

У меня есть некоторые данные с предикторами и двоичной целью. Например:

df <- data.frame(a=sort(sample(1:100,30)), b= sort(sample(1:100,30)), 
                 target=c(rep(0,11),rep(1,4),rep(0,4),rep(1,11)))

Я обучил модель логистической регрессии, используя glm()

model1 <- glm(formula= target ~ a + b, data=df, family=binomial)

Теперь я пытаюсь предсказать вывод (например, одних и тех же данных должно хватить)

predict(model1, newdata=df, type="response")

Это генерирует вектор вероятностных чисел. Но я хочу предсказать фактический класс. Я мог бы использовать round() для чисел вероятности, но это предполагает, что все, что меньше 0,5, является классом «0», а все, что выше, — классом «1». Это правильное предположение? Даже когда численность населения каждого класса может быть не равной (или близкой к равной)? Или есть способ оценить этот порог?


r logistic-regression predict glm
person user2175594    schedule 23.04.2014    source источник
comment
существуют разные критерии, например, точка, в которой сумма чувствительности и специфичности максимальна, см., например, этот вопрос: stackoverflow.com/questions/23131897/   -  person adibender    schedule 23.04.2014
comment
@adibender Спасибо! Но было бы, конечно, неправильно использовать порог как долю населения, верно? То есть, если в популяции 30 % случаев — это «0», а 70 % — «1», наивной оценкой будет использование 0,3 в качестве порога. Но это не было бы логичным подходом к этому?   -  person user2175594    schedule 23.04.2014
comment
Вы можете найти отличное руководство по этому вопросу здесь: hopstat.wordpress.com/2014/12/19/   -  person pbahr    schedule 08.01.2018

Ответы (6)


arrow_upward
6
arrow_downward

Наилучшая пороговая (или граничная) точка для использования в GL-моделях — это точка, которая максимизирует специфичность и чувствительность. Эта пороговая точка может не давать наивысшего прогноза в вашей модели, но она не будет смещена в сторону положительных или отрицательных результатов. Пакет ROCR содержит функции, которые могут помочь вам в этом. проверьте функцию performance() в этом пакете. Это даст вам то, что вы ищете. Вот картинка того, что вы ожидаете получить:

введите здесь описание изображения

Найдя точку отсечки, я обычно сам пишу функцию, чтобы найти количество точек данных, значение прогноза которых превышает отсечку, и сопоставить его с группой, к которой они принадлежат.

person Error404    schedule 23.04.2014
comment
не могли бы вы предоставить более конкретный код, который будет генерировать приведенный выше график? Кроме того, как значения отсечки могут находиться в диапазоне от 0 до 14 для вероятностей, принимающих значения от 0 до 1? - person Kasia Kulma; 29.09.2016
comment
Ниже я добавил подходы baseR/ggplot! - person user61871; 03.12.2018

arrow_upward
4
arrow_downward

Золотым стандартом для определения хороших параметров модели, в том числе "какое пороговое значение я должен установить" для логистической регрессии, является перекрестная проверка.

Общая идея состоит в том, чтобы удержать одну или несколько частей обучающей выборки и выбрать порог, который максимизирует количество правильных классификаций в этой удержанной выборке, но Википедия может дать вам гораздо больше подробностей.

person merlin2011    schedule 23.04.2014
comment
Так как мы будем настраивать пороговый параметр для данных перекрестной проверки, якобы, потребуется третий набор для оценки, чтобы сообщить о объективной ожидаемой ошибке? - person user2175594; 23.04.2014
comment
@ user2175594, Да, верно. Традиционно у вас будет как минимум три отдельных раздела ваших данных: обучение, проверка и тест (оценка). Однако, если вы делаете что-то вроде перекрестной проверки в k-кратном порядке, то обучение и проверка представляют собой, по сути, один и тот же набор, перераспределенный несколькими способами. - person merlin2011; 23.04.2014

arrow_upward
2
arrow_downward

Поковырялся, пытаясь воспроизвести первый график. Учитывая объект predictions <- prediction(pred,labels), тогда:

базовый подход

plot(unlist(performance(predictions, "sens")@x.values), unlist(performance(predictions, "sens")@y.values), 
     type="l", lwd=2, ylab="Specificity", xlab="Cutoff")
par(new=TRUE)
plot(unlist(performance(predictions, "spec")@x.values), unlist(performance(predictions, "spec")@y.values), 
     type="l", lwd=2, col='red', ylab="", xlab="")
axis(4, at=seq(0,1,0.2),labels=z)
mtext("Specificity",side=4, padj=-2, col='red')

введите здесь описание изображения

подход ggplot2

sens <- data.frame(x=unlist(performance(predictions, "sens")@x.values), 
                   y=unlist(performance(predictions, "sens")@y.values))
spec <- data.frame(x=unlist(performance(predictions, "spec")@x.values), 
                   y=unlist(performance(predictions, "spec")@y.values))

sens %>% ggplot(aes(x,y)) + 
  geom_line() + 
  geom_line(data=spec, aes(x,y,col="red")) +
  scale_y_continuous(sec.axis = sec_axis(~., name = "Specificity")) +
  labs(x='Cutoff', y="Sensitivity") +
  theme(axis.title.y.right = element_text(colour = "red"), legend.position="none")

введите здесь описание изображения

person user61871    schedule 03.12.2018

arrow_upward
1
arrow_downward

Чтобы получить порог в данных с ближайшими значениями чувствительности и специфичности (т. е. пересечение на графиках выше) программно, вы можете использовать этот код, который довольно близок:

predictions = prediction(PREDS, LABELS)

sens = cbind(unlist(performance(predictions, "sens")@x.values), unlist(performance(predictions, "sens")@y.values))
spec = cbind(unlist(performance(predictions, "spec")@x.values), unlist(performance(predictions, "spec")@y.values))
sens[which.min(apply(sens, 1, function(x) min(colSums(abs(t(spec) - x))))), 1]
person Adam Waring    schedule 16.04.2020

arrow_upward
0
arrow_downward

В функции PresenceAbsence::optimal.thresholds пакета PresenceAbsence реализовано 12 методов.

Это также рассматривается в Freeman, EA, & Moisen, GG (2008). Сравнение эффективности пороговых критериев для бинарной классификации с точки зрения прогнозируемой распространенности и каппа. Экологическое моделирование, 217(1-2), 48-58.

person irudnyts    schedule 13.11.2019

arrow_upward
-2
arrow_downward

Вы можете попробовать следующее:

perfspec <- performance(prediction.obj = pred, measure="spec", x.measure="cutoff")

plot(perfspec)

par(new=TRUE)

perfsens <- performance(prediction.obj = pred, measure="sens", x.measure="cutoff")

plot(perfsens)
person Dipayan Sarkar    schedule 07.01.2017