Beitrag jeder Kovariate zu einer einzelnen Vorhersage in einem logistischen Regressionsmodell

Nehmen wir zum Beispiel an, wir haben ein logistisches Regressionsmodell, das die Wahrscheinlichkeit ausgibt, dass ein Patient eine bestimmte Krankheit entwickelt, die auf vielen Kovariaten basiert.

Wir können uns ein Bild von der Größe und Richtung des Effekts jeder Kovariate im Allgemeinen machen, indem wir die Koeffizienten des Modells untersuchen und die Änderung des Odds-Ratio berücksichtigen.

Was ist, wenn wir für einen einzelnen Patienten wissen möchten, welche seiner größten Risikofaktoren / welche die größten Faktoren zu seinen Gunsten sind? Ich interessiere mich besonders für diejenigen, gegen die der Patient tatsächlich etwas tun könnte.

Was ist der beste Weg, dies zu tun?

Die Art und Weise, wie ich derzeit überlege, wird im folgenden R-Code (aus diesem Thread entnommen ) erfasst :

#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')])

Ich denke darüber nach, zusätzlich zu schauen

this.student.prediction.list <- this.student.predictors * coef(data.model)

und versuchen, die Informationen aus den einzelnen Addenden der Summe herauszuholen, die die Wahrscheinlichkeitsschätzung darstellt, aber ich bin mir nicht sicher, wie ich das machen soll.

Ich könnte schauen

• Welche Variablen leisten den größten absoluten Beitrag zur Wahrscheinlichkeitsschätzung und nehmen diese als die größten Risikofaktoren an.
• Welche Variablen unterscheiden sich am stärksten von ihrem mittleren Anteil, dh sehen Sie, welchen Anteil jede Variable im Durchschnitt zur Wahrscheinlichkeitsschätzung beiträgt, und sehen Sie, welche Variablen sich in dieser speziellen Beobachtung um den größten Betrag von diesem Anteil unterscheiden
• Eine Kombination davon: Gewichten Sie die absolute Differenz zwischen dem mittleren Anteil und dem beobachteten Anteil mit dem mittleren Anteil und nehmen Sie die Variablen mit den größten gewichteten Werten

Welche davon sind am sinnvollsten? Wäre einer dieser Ansätze ein vernünftiger Weg, um die Frage zu beantworten?

Außerdem möchte ich wissen, wie ich Konfidenzintervalle für die additiven Beiträge einzelner Kovariaten zur Wahrscheinlichkeitsschätzung erhalten kann.

Sie können die predictFunktion in R verwenden. Rufen Sie sie mit auf type='terms'und Sie erhalten den Beitrag jedes Terms im Modell (der Koeffizient multipliziert mit dem Variablenwert). Dies wird auf der Log-Odds-Skala sein.

Eine weitere Option ist die Verwendung der TkPredictFunktion aus dem TeachingDemos-Paket. Dies zeigt eine grafische Darstellung des vorhergesagten Werts gegenüber einem der Prädiktoren. Anschließend kann der Benutzer den Wert der verschiedenen Prädiktoren interaktiv ändern, um zu sehen, wie sich dies auf die Vorhersage auswirkt.