Как получить граф ggplot для имитации логистической регрессии в примере из Википедии в R?

Question

(добавлен воспроизводимый пример)

Я пытался имитировать пример логистической регрессии из Википедии "Вероятность сдачи экзамена в зависимости от количества часов обучения" здесь :

Мне не удалось получить такой же график ggplot на этой странице, и я не мог понять, почему.

df <- data.frame(hour=c(0.50,0.75,1.00,1.25,1.50,1.75,1.75,2.00,2.25,2.50,2.75,3.00,3.25,3.50,4.00,4.25,4.50,4.75,5.00,5.50), pass=c(0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,1,1,1))

df
#   hour pass
#1   0.50    0
#2   0.75    0
#3   1.00    0
#4   1.25    0
#5   1.50    0
#6   1.75    0
#7   1.75    1
#8   2.00    0
#9   2.25    1
#10  2.50    0
#11  2.75    1
#12  3.00    0
#13  3.25    1
#14  3.50    0
#15  4.00    1
#16  4.25    1
#17  4.50    1
#18  4.75    1
#19  5.00    1
#20  5.50    1

df$pass <- as.factor(df$pass)
my_fit <- glm(df$pass ~ df$hour, data=df, na.action=na.exclude, family="binomial")
summary(my_fit)

NON GGPLOT PLOT РАБОТАЕТ ОТЛИЧНО:

my_table <- summary(my_fit)     
my_table$coefficients[,1] <- invlogit(coef(my_fit))
my_table

anova(my_fit)
library(pscl); pR2(my_fit) # for McFadden rho^2

plot(df$hour, df$pass, xlab="x", ylab="logit values")

LinearPredictions <- predict(my_fit); LinearPredictions
# LinearPredictions is NOT equal to 0.01666 + 0.81827*(1:20)
# LinearPredictions is NOT equal to -4.0777+1.5046*(1:20)
# LinearPredictions are equal to what (I couldn't solve)?

EstimatedProbability.hat <- exp(LinearPredictions)/(1 + exp(LinearPredictions))
EstimatedProbability.hat

EstimatedProbability <- c(0.25, 0.50, 0.75) # Estimated probabilities for which their x levels are wanted to be found

HoursStudied <- (log(EstimatedProbability/(1- EstimatedProbability)) - my_fit$coefficients[1])/ my_fit$coefficients[2]
HoursStudied.summary <- data.frame(EstimatedProbability, HoursStudied)
HoursStudied.summary

plot(df$hour, EstimatedProbability.hat, xlab="studying hours", ylab="estimated probability (pass)") # , xlim=c(0,6), ylim=c(0,1)
# Add red curve
lines(df$hour, EstimatedProbability.hat, lty=1, col="red")
# Vertical dashes
segments(x0=HoursStudied.summary$HoursStudied, y0=0, x1=HoursStudied.summary$HoursStudied, y1=HoursStudied.summary$EstimatedProbability,
         lty=2, col=c("darkblue","darkred","darkgreen"))
# Horizontal dashes
segments(x0=0, y0=HoursStudied.summary$EstimatedProbability, x1=HoursStudied.summary$HoursStudied, 
y1=HoursStudied.summary$EstimatedProbability, lty=2, col=c("darkblue","darkred","darkgreen"))

legend("bottomright", legend=c("HS0.25", "HS0.50", "HS0.75"), lty=2, col=c("darkblue","darkred","darkgreen"), bty="n", cex=0.75)

СБОЙ УЧАСТКА В GGPLOT:
Я пытался сделать то же самое в ggplot, но не получилось:

df$EstimatedProbabilities <- EstimatedProbability.hat; df
HoursStudied.summary$group <- c('HS0.25','HS0.50','HS0.75')

library(ggplot2)
ggplot(df, aes(x=hour, y=df$pass)) + 
geom_point() + 
geom_line(aes(y=EstimatedProbabilities), colour="black") + 
geom_segment(data=HoursStudied.summary, aes(y=EstimatedProbability,
xend=HoursStudied, yend=EstimatedProbability, col=group), x=-Inf, linetype="dashed") + 
geom_segment(data=HoursStudied.summary, aes(x=HoursStudied,
xend=HoursStudied, yend=EstimatedProbability, col=group), y=-Inf, linetype="dashed")

Проблема: Кривая ggplot такая же, как и у plot, однако, она рисует всю функцию ниже линии y = 0. Почему?

Answer 1

Вопрос усложняет то, что можно сделать простым с geom_smooth. Обратите внимание, что прогнозы на type = "response", следующие этот пост до CrossValidated .

my_fit <- glm(pass ~ hour, data = df, na.action = na.exclude,
              family = "binomial")
pred <- predict(my_fit, type = "response")
pred_df <- data.frame(hour = df$hour, pred)

library(ggplot2)

ggplot(df, aes(x = hour, y = pass)) +
  geom_point() +
  geom_smooth(method = "glm", 
              method.args = list(family = "binomial"),
              se = FALSE) +
  geom_point(data = pred_df, aes(x = hour, y = pred), colour = "blue") +
  geom_hline(data = data.frame(c(0.25, 0.50, 0.75)),
             aes(yintercept = c(0.25, 0.50, 0.75)),
             colour = "darkgrey", linetype = "dashed")

Answer 2

edit: вам нужно, чтобы ваш df$pass был числовым, а не множителем. Я также не стал бы отображать эстетику при первоначальном вызове ggplot, а просто передавал бы их при вызовах geom_point и geom_line.

df$pass <- as.numeric(df$pass) - 1

ggplot(df) +
    geom_point(aes(x=hour,y=pass)) +
    geom_line(aes(x=hour,y=EstimatedProbabilities)) +
    geom_segment(data=HoursStudied.summary, aes(y=EstimatedProbability, xend=HoursStudied, yend=EstimatedProbability, col=group), x=-Inf, linetype="dashed") + 
    geom_segment(data=HoursStudied.summary, aes(x=HoursStudied, xend=HoursStudied, yend=EstimatedProbability, col=group), y=-Inf, linetype="dashed")