Berechnen von Vorhersageintervallen für die logistische Regression

Ich möchte verstehen, wie man Vorhersageintervalle für logistische Regressionsschätzungen erzeugt.

Mir wurde geraten, die Verfahren in Colletts Modeling Binary Data , 2nd Ed S.98-99, zu befolgen. Nachdem predict.glmich dieses Verfahren implementiert und mit Rs verglichen habe, denke ich, dass dieses Buch das Verfahren zum Berechnen von Konfidenzintervallen und nicht von Vorhersageintervallen zeigt.

Die Implementierung der Prozedur von Collett mit einem Vergleich zu predict.glmwird unten gezeigt.

Ich möchte wissen, wie ich von hier aus ein Vorhersageintervall anstelle eines Konfidenzintervalls erzeuge.

#Derived from Collett 'Modelling Binary Data' 2nd Edition p.98-99
#Need reproducible "random" numbers.
seed <- 67

num.students <- 1000
which.student <- 1

#Generate data frame with made-up data from students:
set.seed(seed) #reset seed
v1 <- rbinom(num.students,1,0.7)
v2 <- rnorm(length(v1),0.7,0.3)
v3 <- rpois(length(v1),1)

#Create df representing students
students <- data.frame(
    intercept = rep(1,length(v1)),
    outcome = v1,
    score1 = v2,
    score2 = v3
)
print(head(students))

predict.and.append <- function(input){
    #Create a vanilla logistic model as a function of score1 and score2
    data.model <- glm(outcome ~ score1 + score2, data=input, family=binomial)

    #Calculate predictions and SE.fit with the R package's internal method
    # These are in logits.
    predictions <- as.data.frame(predict(data.model, se.fit=TRUE, type='link'))

    predictions$actual <- input$outcome
    predictions$lower <- plogis(predictions$fit - 1.96 * predictions$se.fit)
    predictions$prediction <- plogis(predictions$fit)
    predictions$upper <- plogis(predictions$fit + 1.96 * predictions$se.fit)


    return (list(data.model, predictions))
}

output <- predict.and.append(students)

data.model <- output[[1]]

#summary(data.model)

#Export vcov matrix 
model.vcov <- vcov(data.model)

# Now our goal is to reproduce 'predictions' and the se.fit manually using the vcov matrix
this.student.predictors <- as.matrix(students[which.student,c(1,3,4)])

#Prediction:
this.student.prediction <- sum(this.student.predictors * coef(data.model))
square.student <- t(this.student.predictors) %*% this.student.predictors
se.student <- sqrt(sum(model.vcov * square.student))

manual.prediction <- data.frame(lower = plogis(this.student.prediction - 1.96*se.student), 
    prediction = plogis(this.student.prediction), 
    upper = plogis(this.student.prediction + 1.96*se.student))

print("Data preview:")
print(head(students))
print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by Collett's procedure:"))
manual.prediction
print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by R's predict.glm:"))    
print(output[[2]][which.student,c('lower','prediction','upper')])

Vorhersageintervalle sagen voraus, wo die tatsächlichen Antwortdatenwerte mit einer gegebenen Wahrscheinlichkeit abfallen werden. Da die möglichen Werte der Antwort eines Logistikmodells auf 0 und 1 beschränkt sind, ist das 100% -Vorhersageintervall daher . Keine anderen Intervalle sind für eine Vorhersage mit logistischer Regression wirklich sinnvoll. Da es immer das gleiche Intervall ist, ist es im Allgemeinen nicht interessant genug, um es zu generieren oder zu diskutieren.$0 <= y <= 1$