Wykreślanie przedziałów ufności dla przewidywanych prawdopodobieństw z regresji logistycznej

Ok, mam regresję logistyczną i wykorzystałem tę predict()funkcję do opracowania krzywej prawdopodobieństwa na podstawie moich oszacowań.

## LOGIT MODEL:
library(car)
mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))

## PROBABILITY CURVE:
all.x <- expand.grid(won=unique(won), bid=unique(bid))
y.hat.new <- predict(mod1, newdata=all.x, type="response")
plot(bid<-000:1000,predict(mod1,newdata=data.frame(bid<-c(000:1000)),type="response"), lwd=5, col="blue", type="l")

To świetnie, ale jestem ciekawy, jak wykreślić przedziały ufności dla prawdopodobieństw. Próbowałem, plot.ci()ale nie miałem szczęścia. Czy ktoś może wskazać mi kilka sposobów, aby to zrobić, najlepiej z carpakietem lub bazą R.

Użyty kod szacuje model regresji logistycznej za pomocą glmfunkcji. Nie podałeś danych, więc po prostu je uzupełnię.

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

Model regresji logistycznej modeluje związek między binarną zmienną odpowiedzi a, w tym przypadku, jednym ciągłym predyktorem. Wynikiem jest prawdopodobieństwo przekształcone logitem jako liniowa zależność od predyktora. W twoim przypadku wynik jest binarną odpowiedzią na wygraną lub nie na wygraną w grach hazardowych i jest przewidziany na podstawie wartości zakładu. Współczynniki z mod1są podane w zarejestrowanych szansach (trudnych do interpretacji), zgodnie z:

Aby przekonwertować zarejestrowane kursy na prawdopodobieństwa, możemy przetłumaczyć powyższe na

Możesz użyć tych informacji do skonfigurowania fabuły. Po pierwsze, potrzebujesz zakresu zmiennej predykcyjnej:

plotdat <- data.frame(bid=(0:1000))

Następnie za pomocą predictmożesz uzyskać prognozy na podstawie swojego modelu

preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)

Należy pamiętać, że dopasowane wartości można również uzyskać poprzez

mod1$fitted

Określając se.fit=TRUE, otrzymujesz również błąd standardowy związany z każdą dopasowaną wartością. Wynikiem data.framejest macierz z następującymi składnikami: dopasowane predykcje ( fit), szacowane błędy standardowe ( se.fit) i skalar dający pierwiastek kwadratowy z dyspersji użytej do obliczenia błędów standardowych ( residual.scale). W przypadku dwumianowego logit, wartość będzie wynosić 1 (który można zobaczyć wpisując preddat$residual.scalew R). Jeśli chcesz zobaczyć przykład tego, co dotychczas obliczyłeś, możesz wpisać head(data.frame(preddat)).

Następnym krokiem jest skonfigurowanie fabuły. Najpierw chcę ustawić pusty obszar kreślenia z parametrami:

with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))

Teraz możesz zobaczyć, gdzie ważne jest, aby wiedzieć, jak obliczyć dopasowane prawdopodobieństwa. Możesz narysować linię odpowiadającą dopasowanym prawdopodobieństwom zgodnie z drugim wzorem powyżej. Za pomocą preddat data.framemożesz przekonwertować dopasowane wartości na prawdopodobieństwa i użyć tego do wykreślenia linii względem wartości zmiennej predykcyjnej.

with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))

Na koniec odpowiedz na pytanie, przedziały ufności można dodać do wykresu, obliczając prawdopodobieństwo dopasowanych wartości +/- 1.96razy błąd standardowy:

with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

Powstały wykres (z losowo wygenerowanych danych) powinien wyglądać mniej więcej tak:

Na wszelki wypadek, oto cały kod w jednym kawałku:

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))
plotdat <- data.frame(bid=(0:1000))
preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)
with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))
with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))
with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

(Uwaga: jest to mocno zredagowana odpowiedź, która ma na celu uczynienie jej bardziej adekwatną do stats.stackexchange.)

Oto modyfikacja rozwiązania @ smillig. Używam tutaj narzędzi tidyverse, a także linkinvfunkcji, która jest częścią obiektu modelu GLM mod1. W ten sposób nie musisz ręcznie odwracać funkcji logistycznej, a to podejście będzie działać bez względu na to, jaki konkretny GLM pasuje.

library(tidyverse)
library(magrittr)


set.seed(1234)

# create fake data on gambling. Does prob win depend on bid size? 
mydat <- data.frame(
  won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
  bid=runif(250, min=0, max=1000)
)

# logistic regression model: 
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

# new predictor values to use for prediction: 
plotdat <- data.frame(bid=(0:1000))

# df with predictions, lower and upper limits of CIs: 
preddat <- predict(mod1,
               type = "link",
               newdata=plotdat,
               se.fit=TRUE) %>% 
  as.data.frame() %>% 
  mutate(bid = (0:1000), 

         # model object mod1 has a component called linkinv that 
         # is a function that inverts the link function of the GLM:
         lower = mod1$family$linkinv(fit - 1.96*se.fit), 
         point.estimate = mod1$family$linkinv(fit), 
         upper = mod1$family$linkinv(fit + 1.96*se.fit)) 


# plotting with ggplot: 
preddat %>% ggplot(aes(x = bid, 
                   y = point.estimate)) + 
  geom_line(colour = "blue") + 
  geom_ribbon(aes(ymin = lower,
                  ymax = upper), 
              alpha = 0.5) + 
  scale_y_continuous(limits = c(0,1))