Vorhersage- und Konfidenzintervalle für die logistische Regression

Question

Im Folgenden finden Sie eine Reihe von fiktiven Wahrscheinlichkeitsdaten, die ich mit einer threshold of 0.5 in binomial konvertiert habe. Ich führte ein glm() Modell auf den diskreten Daten durch, um zu testen, ob die von glm() zurückgegebenen Intervalle "mittlere Vorhersageintervalle" ("Confidence Interval") oder "Punktvorhersageintervalle" ("Prediction Interval") waren. Aus dem folgenden Diagramm geht hervor, dass die zurückgegebenen Intervalle die letzten sind - "Point Prediction Intervals"; Beachten Sie, dass bei einer Stichprobe von 95% 2/20 Punkte außerhalb der Linie liegen.

Wenn dies tatsächlich der Fall ist, wie generiere ich das 'mittlere Vorhersage-Intervall' (d. H. "Confidence Intervals") in R für einen Binomial-Datensatz mit 0 und 1 mit Glm()? Bitte zeigen Sie Ihren Code und meinen ähnlichen Plan mit der Fit-Linie, mit den gegebenen Wahrscheinlichkeiten, Konfidenzintervallen und Prädiktionsintervallen.

# Fictitious data 
xVal <- c(15,15,17,18,32,33,41,42,47,50, 
     53,55,62,63,64,65,66,68,70,79, 
     94,94,94,95,98) 
randRatio <- c(.01,.03,.05,.04,.01,.2,.1,.08,.88,.2, 
       .2,.99,.49,.88,.2,.88,.66,.87,.66,.90, 
       .98,.88,.95,.95,.95) 
# Converted to binomial 
randBinom <- ifelse(randRatio < .5, 0, 1) 

# Data frame for model 
binomData <- data.frame(
    randBinom = randBinom, 
    xVal = xVal 
) 

# Model 
mode1 <- glm(randBinom~ xVal, data = binomData, family = binomial(link = "logit")) 

# Predict all points in xVal range 
frame <- data.frame(xVal=(0:100)) 
predAll <- predict(mode1, newdata = frame,type = "link", se.fit=TRUE) 

# Params for intervals and plot 
confidence <- .95 
score <- qnorm((confidence/2) + .5) 
frame <- data.frame(xVal=(0:100)) 

#Plot 
with(binomData, plot(xVal, randBinom, type="n", ylim=c(0, 1), 
       ylab = "Probability", xlab="xVal")) 
lines(frame$xVal, plogis(predAll$fit), col = "red", lty = 1) 
lines(frame$xVal, plogis(predAll$fit + score * predAll$se.fit), col = "red", lty = 3) 
lines(frame$xVal, plogis(predAll$fit - score * predAll$se.fit), col = "red", lty = 3) 
points(xVal, randRatio, col = "red") # Original probabilities 
points(xVal, randBinom, col = "black", lwd = 3) # Binomial Points used in glm 

Hier ist die Handlung, die vermutlich mit '-Punkt Prädiktionsintervalle' (d.h. "Prediction Intervals") in gestrichelten roten und dem mittleren Sitz in festen rot. Schwarze Punkte stellen die diskrete binomische Daten von der ursprünglichen Wahrscheinlichkeiten in randRatio:

Answer 1

Ich bin nicht sicher, ob Sie für die gerade nach oben Prädiktionsintervalls fragen, aber wenn Sie sind, können Sie es einfach berechnen.

Sie können als solche für das Modell ein traditionelles Konfidenzintervall extrahieren:

confint(model) 

Und dann, wenn Sie eine Prognose ausführen, können Sie ein Prognoseintervall berechnen basierend auf der Vorhersage, wie so:

upper = predAll$fit + 1.96 * predAll$se.fit 
lower = predAll$fit - 1.96 * predAll$se.fit 

Sie nehmen einfach die Vorhersage (an jedem beliebigen Punkt, wenn Sie einen einzelnen Satz von Prädiktorvariablen verwenden) und addieren und subtrahieren 1,96 * absoluten Wert des Standardfehlers. (1,96 Se enthält 97,5% der Normalverteilung und stellt das 95% Intervall wie für die Standardabweichung in der Normalverteilung dar)

Dies ist die gleiche Formel, die Sie für ein traditionelles Konfidenzintervall verwenden würden, außer dass die Verwendung der Standardfehler (im Gegensatz zur Standardabweichung) machen das Intervall breiter, um die Unsicherheit in der Vorhersage selbst zu berücksichtigen.

Update:

Method for plotting prediction invervals courtesy of Rstudio!

Wie gewünscht ... wenn auch nicht von mir gemacht!