Darstellen von Konfidenzintervallen für die vorhergesagten Wahrscheinlichkeiten aus einer logistischen Regression

20

Ok, ich habe eine logistische Regression und habe die predict()Funktion verwendet, um eine Wahrscheinlichkeitskurve basierend auf meinen Schätzungen zu entwickeln.

## 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")

Das ist großartig, aber ich bin gespannt darauf, die Konfidenzintervalle für die Wahrscheinlichkeiten zu zeichnen. Ich habe es versucht, plot.ci()aber kein Glück gehabt. Kann jemand mich auf einige Möglichkeiten hinweisen, um dies zu erreichen, vorzugsweise mit dem carPaket oder der Basis R.

ATMathew
quelle
4
(+1) Als Antwort auf die Stimmen, um das Thema abzuschließen: Anscheinend ist die Basis für diese Stimmen, dass die Frage eine rein softwarebezogene Frage zu stellen scheint ("wie man so und so in R zeichnet"), a Frage, die in der Tat auf SO erscheinen sollte. Beachten Sie jedoch, dass in der aktuellen Antwort statistische Formeln zum Erstellen der Darstellungspunkte enthalten sind. Dies deutet darauf hin, dass ein statistisches Interesse an der Frage besteht, und ich zögere daher, für die Migration zu stimmen. Eine gute Antwort hier würde diesen statistischen Punkt hervorheben und erläutern.
whuber

Antworten:

26

Der von Ihnen verwendete Code schätzt mithilfe der glmFunktion ein logistisches Regressionsmodell . Sie haben keine Daten angegeben, deshalb werde ich nur einige erfinden.

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"))

Ein logistisches Regressionsmodell modelliert die Beziehung zwischen einer binären Antwortvariablen und in diesem Fall einem kontinuierlichen Prädiktor. Das Ergebnis ist eine logittransformierte Wahrscheinlichkeit als lineare Beziehung zum Prädiktor. In Ihrem Fall ist das Ergebnis eine binäre Antwort, die dem Gewinnen oder Nichtgewinnen beim Spielen entspricht, und sie wird durch den Wert des Einsatzes vorhergesagt. Die Koeffizienten von mod1werden in protokollierten Quoten angegeben (die schwer zu interpretieren sind), gemäß:

logit(p)=Log(p(1-p))=β0+β1x1

Um protokollierte Gewinnchancen in Wahrscheinlichkeiten umzuwandeln, können wir das Obige in übersetzen

p=exp(β0+β1x1)(1+exp(β0+β1x1))

Mit diesen Informationen können Sie den Plot einrichten. Zunächst benötigen Sie einen Bereich der Prädiktorvariablen:

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

Anschließend predictkönnen Sie anhand Ihres Modells Vorhersagen abrufen

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

Beachten Sie, dass die angepassten Werte auch über abgerufen werden können

mod1$fitted

Durch die Angabe erhalten se.fit=TRUESie auch den Standardfehler für jeden angepassten Wert. Das Ergebnis data.frameist eine Matrix mit den folgenden Komponenten: den angepassten Vorhersagen ( fit), den geschätzten Standardfehlern ( se.fit) und einem Skalar, der die Quadratwurzel der Dispersion angibt, die zur Berechnung der Standardfehler ( residual.scale) verwendet wird. Bei einem Binomial-Logit ist der Wert 1 (was Sie durch Eingabe preddat$residual.scalevon sehen können R). Wenn Sie ein Beispiel dessen sehen möchten, was Sie bisher berechnet haben, können Sie Folgendes eingeben head(data.frame(preddat)).

Der nächste Schritt ist das Einrichten des Plots. Ich möchte zuerst einen leeren Zeichenbereich mit den Parametern einrichten:

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

Jetzt können Sie sehen, wo es wichtig ist, die angepassten Wahrscheinlichkeiten zu berechnen. Sie können die den angepassten Wahrscheinlichkeiten entsprechende Linie gemäß der obigen zweiten Formel zeichnen. Mit preddat data.framekönnen Sie die angepassten Werte in Wahrscheinlichkeiten umwandeln und damit eine Linie gegen die Werte Ihrer Prädiktorvariablen zeichnen.

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

Beantworten Sie schließlich Ihre Frage. Die Konfidenzintervalle können dem Diagramm hinzugefügt werden, indem die Wahrscheinlichkeit für die angepassten Werte +/- 1.96multipliziert mit dem Standardfehler berechnet wird:

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))

Das resultierende Diagramm (aus den zufällig generierten Daten) sollte ungefähr so ​​aussehen:

Bildbeschreibung hier eingeben

Aus Gründen der Zweckmäßigkeit ist hier der gesamte Code in einem Block:

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))

(Hinweis: Dies ist eine stark bearbeitete Antwort, um sie für stats.stackexchange relevanter zu machen.)

smillig
quelle
Wo ist die Variable se.fitdefiniert?
Makro
In predict(..., se.fit=TRUE).
smillig
(-1) Diese CIs sind für jeden Einzelfall? In diesem Fall ist für ein binäres Ergebnis der einzig sinnvolle CI für eine vorhergesagte Wahrscheinlichkeit [0,1]. Auch wenn dies eine technisch kompetente Antwort sein könnte.
Rolando2
Per @ whubers Kommentar denke ich, dass eine gute Antwort eine Formel enthalten sollte, wie die SE berechnet wird. Könnte vielleicht jemand die Antwort bearbeiten und verbessern?
Heisenberg
1
Ihre Antwort scheint nur das "mittlere Vorhersageintervall" zu geben. Wie würde ich das "Punktvorhersageintervall" hinzufügen?
Bob Hopez
0

Hier ist eine Modifikation der @ smillig-Lösung. Ich verwende hier Tidyverse-Tools und verwende auch die linkinvFunktion, die Teil des GLM- Modellobjekts ist mod1. Auf diese Weise müssen Sie die Logistikfunktion nicht manuell invertieren, und dieser Ansatz funktioniert unabhängig davon, welches spezifische GLM Sie anpassen.

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))
Nayef
quelle
3
Obwohl die Implementierung bei Fragen häufig mit inhaltlichen Inhalten vermischt wird, soll es sich bei uns um eine Website handeln, die Informationen zu Statistiken, maschinellem Lernen usw. und nicht um Code enthält. Es kann auch gut sein, Code bereitzustellen, aber bitte erarbeiten Sie Ihre inhaltliche Antwort in Textform für Personen, die diese Sprache nicht gut genug lesen, um die Antwort aus dem Code zu erkennen und zu extrahieren.
gung - Wiedereinsetzung von Monica