Berechnen von Vorhersageintervallen für die logistische Regression

20

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')])
Carbocation
quelle
Eine grundlegende Frage, warum wird sqrt (sum (model.vcov * square.student)) als Standardfehler angenommen? Ist es nicht die Standardabweichung und muss durch sqrt (n) geteilt werden? Wenn ja, welches n sollte verwendet werden, n sollte zur Anpassung an das Modell verwendet werden, oder n des neuen Datenrahmens, der zur Vorhersage verwendet wird?
Rafael

Antworten:

6

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

Greg Snow
quelle
6
Ich suche nach einem Vorhersageintervall von 95% für eine Vorhersage, die sich im Log-Odds-Bereich befindet. Später transformiere ich das in Wahrscheinlichkeitsraum. Ein Vorhersageintervall von 100% wäre für kein Verfahren interessant, oder? Ein Vorhersageintervall von 100% für die lineare Regression würde beispielsweise -Inf to Inf ... enthalten. Wie Sie in meinem Code sehen können, wird das Vorhersageintervall jedenfalls im Protokollquotenraum berechnet, der später in den Wahrscheinlichkeitsraum transformiert wird . Daher halte ich meine Frage nicht für sinnlos.
Carbo
2
Die Log-Odds können in eine Wahrscheinlichkeit umgewandelt werden und Sie können ein Konfidenzintervall für die Wahrscheinlichkeit (oder die Log-Odds) berechnen. Für die Antwortvariable gibt es jedoch ein Vorhersageintervall von 0 oder 1. Wenn Ihr Ergebnis das Überleben mit 0 = tot und 1 = lebendig ist, können Sie die Wahrscheinlichkeit voraussagen, dass Sie für eine bestimmte Menge von Kovariaten am Leben sind, und ein Konfidenzintervall für berechnen diese Wahrscheinlichkeit. Aber das Ergebnis ist 0/1. Sie können keinen Patienten haben, der zu 62% am Leben ist. Es muss 0 oder 1 sein. Die einzig möglichen Vorhersageintervalle sind 0-0, 0-1 und 1-1 (das ist warum sich die meisten Menschen an Vertrauensbereiche halten).
Greg Snow
8
Wenn Sie eine Situation haben, in der die Antwort binomisch ist (was unter den gleichen Bedingungen ein Aggregat von 0-1s sein kann), kann ein Vorhersageintervall sinnvoll sein.
Glen_b
7
Logistische Regression ist die Regression einer Wahrscheinlichkeit, bei der versucht wird, die Wahrscheinlichkeit eines Ereignisses als Funktion von Regressorvariablen zu modellieren. Vorhersageintervalle in dieser Einstellung werden als Intervalle auf der Wahrscheinlichkeitsskala oder der logarithmischen Wahrscheinlichkeitsskala herangezogen, um perfekte Werte zu erzielen.
kjetil b halvorsen
2
@Cesar, die Vorhersageintervallformel wird abgeleitet, indem angenommen wird, dass Y normal über die Linie verteilt ist, aber in der logistischen Regression haben wir keine Normalverteilung, sondern ein Bernoulli oder Binomial. Das Anwenden der Formeln auf dieser Seite würde entweder zu einem Konfidenzintervall führen (kann dies bereits tun) oder zu einem künstlich erweiterten Konfidenzintervall, das nicht der Definition eines Vorhersageintervalls entspricht (Vorhersage der tatsächlichen Ergebnisse auf der ursprünglichen Ergebnisskala). Wie bereits in Glen_b erwähnt, kann ein Vorhersageintervall sinnvoll sein, wenn das Ergebnis wirklich binomisch ist.
Greg Snow