Отбросьте year из вашей модели, поскольку она не имеет изменений, заново установите модель, а затем передайте flights в качестве аргумента newdata методу predict() модели.
Пример, используятермины и сокращения со страницы Википедии на ROC :
library(nycflights13)
late_arrival<- flights$arr_delay>50
late_arrival.lr <- glm(late_arrival~carrier+dep_delay+month, data=flights, family='binomial')
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
fit <- predict(late_arrival.lr, newdata = flights, type = "response")
d <- data.frame(late_arrival, fit)
# "Confusion matrix" of actual vs predicted outcomes
# for a cutpoint of 0.5:
xtabs(~ late_arrival + I(fit > 0.5), data = d)
#> I(fit > 0.5)
#> late_arrival FALSE TRUE
#> FALSE 290637 3091
#> TRUE 7386 26232
# Now do this for a range of cutpoints.
# Sensitivity = true positive rate = TPR
# Specificity = true negative rate = TNR
# 1 - Specificity = false positive rate = FPR = 1 - TNR
# The ROC plot is
# x = 1 - Specificity = FPR
# y = Sensitivity = TPR
fun <- function(cutpoint) {
pred <- d$fit > cutpoint
# cm = "confusion matrix"
cm <- xtabs(~ late_arrival + I(fit > cutpoint), data = d)
cm <- as.list(cm)
names(cm) <- c("TN", "FN", "FP", "TP")
sens <- with(cm, TP / (TP + FN))
spec <- with(cm, TN / (TN + FP))
return(data.frame(cutpoint, sens, spec))
}
# Example output:
fun(0.5)
#> cutpoint sens spec
#> 1 0.5 0.7802963 0.9894767
cutpoints <- seq(0.02, 0.98, by = 0.02)
# This does
# rbind(fun(cutpoints[1]), fun(cutpoints[2], ...)
roc <- do.call(rbind, lapply(cutpoints, fun))
plot(1 - roc$spec, roc$sens, type = "b",
xlab = "False positive rate (1 - specificity)",
ylab = "True positive rate (sensitivity)",
xlim = c(0, 1),
ylim = c(0, 1))
Создано в 2019-04-07 по представпакет (v0.2.1.9000)
Обратите внимание, что перед тем, как ответить на ваш главный вопрос, необходимо решить несколько вопросов:
Эффект year в вашем примере оценивается в NA, потому что в этой переменной нет изменений , поэтому оценить ее влияние невозможно.
> unique(flights$year)
[1] 2013
Если вы отбросите этот предиктор и заново подгонитевыходные данные имеют смысл (имеется в виду, что нет NA или огромных стандартных ошибок):
> late_arrival.lr <- glm(late_arrival~carrier+dep_delay+month, data=flights, family='binomial')
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
> coef(summary(late_arrival.lr))
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.325540101 0.0564526220 -94.3364527 0.000000e+00
carrierAA 0.335139676 0.0622536491 5.3834543 7.306979e-08
carrierAS -0.980666348 0.3701250164 -2.6495544 8.059801e-03
carrierB6 0.524971196 0.0542918253 9.6694335 4.066226e-22
carrierDL 0.406813418 0.0576767561 7.0533339 1.746810e-12
carrierEV 0.350366432 0.0535144496 6.5471370 5.865056e-11
carrierF9 0.776012126 0.2084826127 3.7221911 1.975015e-04
carrierFL 0.773647203 0.1077982499 7.1768067 7.135846e-13
carrierHA -2.225896541 0.8684691013 -2.5630118 1.037685e-02
carrierMQ 0.847415433 0.0601677914 14.0842037 4.749822e-45
carrierOO 0.232324503 1.3043323784 0.1781176 8.586307e-01
carrierUA 0.157191477 0.0549977051 2.8581461 4.261241e-03
carrierUS 0.649304471 0.0697493204 9.3091154 1.289014e-20
carrierVX 0.237994726 0.1131585684 2.1031967 3.544858e-02
carrierWN 0.032542799 0.0736491439 0.4418626 6.585887e-01
carrierYV 0.861814625 0.2373042135 3.6316870 2.815745e-04
dep_delay 0.089655081 0.0004428296 202.4595603 0.000000e+00
month 0.005089147 0.0032449949 1.5683066 1.168096e-01
Предупреждение fitted probabilities numerically 0 or 1 occurred часто означает, что результат идеально предсказывается одним из ваших непрерывно оцениваемых значенийпредсказатели .Например:
> x <- c(1, 2, 3)
> y <- c(0, 0, 1)
> coef(summary(glm(y ~ x, family="binomial")))
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Estimate Std. Error z value Pr(>|z|)
(Intercept) -115.57626 226884.08 -0.0005094067 0.9995936
x 46.34447 94156.73 0.0004922056 0.9996073
Здесь наилучшей оценкой будет
P (y = 1) = (0, если x <порог), иначе 1 </p>
, но это поднимает двачисленные задачи:
- Обычно сигмовидная кривая P (y = 1) против x теперь должна быть шаговой функцией .Это требует бесконечно крутой сигмоидальной формы, поэтому «наклон» относительно x стремится к бесконечности.
- Любой порог между 2 и 3 будет работать одинаково хорошо, поэтому невозможно определить одну лучшую оценку для перехвата.
В случае flights, однако, я думаю, что предупреждение просто означает, что оно говорит : Некоторые прогнозы настолько уверены, что любые нюансы будут потеряны при ошибке округления.
При проверке, действительно ли late_arrival действительно может быть точно предсказано одной переменной x, я использовал следующий код:
# Make warnings print as they appear.
# options() returns the previous settings, and we store it
warn <- options(warn = 1)$warn
for (i in c("carrier", "dep_delay", "month", "year")) {
print(i)
glm(late_arrival~flights[[i]], family='binomial')
}
# Restore the previous warning setting
options(warn = warn)
, который печатает
[1] "carrier"
[1] "dep_delay"
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
[1] "month"
[1] "year"
Но plot(flights$dep_delay, late_arrival) (занимает несколько секунд) показывает, что на самом деле нет полного разделения, когда все late_arrival происходят для dep_delay> некоторого порога.