Прогнозируемые значения логистической регрессии из glm и stat_smooth в ggplot2 различаются

Я пытаюсь сделать этот график логистической регрессии в ggplot2.

df <- structure(list(y = c(2L, 7L, 776L, 19L, 12L, 26L, 7L, 12L, 8L,
24L, 20L, 16L, 12L, 10L, 23L, 20L, 16L, 12L, 18L, 22L, 23L, 22L,
13L, 7L, 20L, 12L, 13L, 11L, 11L, 14L, 10L, 8L, 10L, 11L, 5L,
5L, 1L, 2L, 1L, 1L, 0L, 0L, 0L), n = c(3L, 7L, 789L, 20L, 14L,
27L, 7L, 13L, 9L, 29L, 22L, 17L, 14L, 11L, 30L, 21L, 19L, 14L,
22L, 29L, 28L, 28L, 19L, 10L, 27L, 22L, 18L, 18L, 14L, 23L, 18L,
12L, 19L, 15L, 13L, 9L, 7L, 3L, 1L, 1L, 1L, 1L, 1L), x = c(18L,
19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L,
32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L,
45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 59L,
62L, 63L, 66L)), .Names = c("y", "n", "x"), class = "data.frame", row.names = c(NA,
-43L))


mod.fit <- glm(formula = y/n ~ x, data = df, weight=n, family = binomial(link = logit),
        na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T))
summary(mod.fit)

Pi <- c(0.25, 0.5, 0.75)
LD <- (log(Pi /(1-Pi))-mod.fit$coefficients[1])/mod.fit$coefficients[2]
LD.summary <- data.frame(Pi , LD)
LD.summary


plot(df$x, df$y/df$n, xlab = "x", ylab = "Estimated probability")

lin.pred <- predict(mod.fit)
pi.hat <- exp(lin.pred)/(1 + exp(lin.pred))
lines(df$x, pi.hat, lty = 1, col = "red")


segments(x0 = LD.summary$LD, y0 = -0.1, x1 = LD.summary$LD, y1 = LD.summary$Pi,
         lty=2, col=c("darkblue","darkred","darkgreen"))
segments(x0 = 15, y0 = LD.summary$Pi, x1 = LD.summary$LD, y1 = LD.summary$Pi,
         lty=2, col=c("darkblue","darkred","darkgreen"))
legend("bottomleft", legend=c("LD25", "LD50", "LD75"), lty=2, col=c("darkblue","darkred","darkgreen"), bty="n", cex=0.75)

Вот моя попытка с ggplot2

library(ggplot2)

p <- ggplot(data = df, aes(x = x, y = y/n)) +
            geom_point() +
            stat_smooth(method = "glm", family = "binomial")

p <- p + geom_segment(aes(
                            x = LD.summary$LD
                          , y = 0
                          , xend = LD.summary$LD
                          , yend = LD.summary$Pi
                         )
                         , colour="red"
                       )

p <- p + geom_segment(aes(
                            x = 0
                          , y = LD.summary$Pi
                          , xend = LD.summary$LD
                          , yend = LD.summary$Pi
                         )
                         , colour="red"
                       )

print(p)

Вопросы

Предсказанные значения для glm и stat_smooth выглядят иначе. Эти два метода дают разные результаты, или я что-то пропускаю здесь.
Мой график ggplot2 не является таким же базовым графиком R.
Как использовать разные цвета для сегментов линии в ggplot2?
А как поставить легенду в ggplot2?

Заранее благодарим за помощь и время. Благодаря

Ответ 1

Просто несколько незначительных дополнений к ответу @mathetmatical.coffee. Как правило, geom_smooth не должен заменять фактическое моделирование, поэтому он может казаться неудобным в тех случаях, когда вы хотите использовать определенный вывод, который вы получите от glm и тому подобное. Но на самом деле все, что нам нужно сделать, это добавить установленные значения в наш кадр данных:

df$pred <- pi.hat
LD.summary$group <- c('LD25','LD50','LD75')

ggplot(df,aes(x = x, y = y/n)) + 
    geom_point() + 
    geom_line(aes(y = pred),colour = "black") + 
    geom_segment(data=LD.summary, aes(y = Pi,
                                      xend = LD,
                                      yend = Pi,
                                      col = group),x = -Inf,linetype = "dashed") + 
    geom_segment(data=LD.summary,aes(x = LD,
                                     xend = LD,
                                     yend = Pi,
                                     col = group),y = -Inf,linetype = "dashed")

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

Урок здесь состоит в том, что если все, что вы хотите сделать, это добавить плавный график, и от него ничего не зависит, используйте geom_smooth. Если вы хотите обратиться к выходному сигналу с установленной модели, то, как правило, она легче подходит для модели вне ggplot, а затем для построения графика.

Ответ 2

Измените свой LD.summary, чтобы включить новый столбец с group (или соответствующей меткой).

LD.summary$group <- c('LD25','LD50','LD75')

Затем измените свои команды geom_segment, чтобы иметь col=LD.summary$group в нем (и удалите colour="red"), который отображает каждый сегмент другого цвета и добавляет легенду:

geom_segment( aes(...,col=LD.summary$group) )

Кроме того, чтобы избежать необходимости делать LD.summary$xxx все время, отправьте data=LD.summary на geom_segment:

geom_segment(data=LD.summary, aes(x=0, y=Pi,xend=LD, yend=Pi, colour=group) )

Что касается того, почему графики не совсем то же самое, в базовом R-графике ось x переходит от ~ 20 вперед, тогда как в ggplot она идет от нуля вперед. Это связано с тем, что ваш второй geom_segment начинается с x=0. Чтобы исправить это, вы можете изменить x=0 на x=min(df$x).

Чтобы ваша метка оси y использовала + scale_y_continuous('Estimated probability').

Вкратце:

LD.summary$group <- c('LD25','LD50','LD75')
p <- ggplot(data = df, aes(x = x, y = y/n)) +
            geom_point() +
            stat_smooth(method = "glm", family = "binomial") +
            scale_y_continuous('Estimated probability')    # <-- add y label
p <- p + geom_segment(data=LD.summary, aes( # <-- data=Ld.summary
                            x = LD
                          , y = 0
                          , xend = LD
                          , yend = Pi
                          , col = group     # <- colours
                         )
                       )    
p <- p + geom_segment(data=LD.summary, aes( # <-- data=Ld.summary
                            x = min(df$x)   # <-- don't plot all the way to x=0
                          , y = Pi
                          , xend = LD
                          , yend = Pi
                          , col = group     # <- colours
                         )
                       )
print(p)

который дает: