Contribution de chaque covariable à une seule prédiction dans un modèle de régression logistique

8

Supposons, par exemple, que nous avons un modèle de régression logistique qui génère la probabilité qu'un patient développe une maladie particulière sur la base de nombreuses covariables.

Nous pouvons avoir une idée de l'ampleur et de la direction de l'effet de chaque covariable en général en examinant les coefficients du modèle et en considérant le changement du rapport de cotes.

Et si nous voulons savoir pour un seul patient quels sont ses plus grands facteurs de risque / les plus grands facteurs en sa faveur. Je suis particulièrement intéressé par ceux pour lesquels le patient pourrait réellement faire quelque chose.

Quelle est la meilleure façon de procéder?

La façon dont je considère actuellement est capturée dans le code R suivant (extrait de ce fil ):

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

J'envisage de regarder en plus

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

et essayer d'obtenir les informations des additifs individuels de la somme qui est l'estimation de la probabilité, mais je ne sais pas comment le faire.

Je pourrais regarder

  • Quelles variables apportent la plus grande contribution absolue à l'estimation de probabilité et prennent celles-ci comme les facteurs de risque les plus importants.
  • Quelles variables diffèrent le plus par rapport à leur proportion moyenne, c'est-à-dire voir quelle proportion chaque variable contribue à l'estimation de probabilité en moyenne et voir quelles variables diffèrent de cette proportion par le plus grand nombre dans cette observation particulière
  • Une combinaison de ceux-ci: pondérer la différence absolue entre la proportion moyenne et la proportion observée par la proportion moyenne et prendre les variables avec les plus grandes valeurs pondérées

Lesquelles ont le plus de sens? L'une de ces approches serait-elle une façon raisonnable de répondre à la question?

De plus, j'aimerais savoir comment obtenir des intervalles de confiance pour les contributions additives de covariables individuelles à l'estimation de probabilité.

Dave
la source

Réponses:

10

Vous pouvez utiliser la predictfonction dans R. Appelez-la avec type='terms'et elle vous donnera la contribution de chaque terme dans le modèle (le coefficient multiplié par la valeur variable). Ce sera sur l'échelle log-odds.

Une autre option consiste à utiliser la TkPredictfonction du package TeachingDemos. Cela affichera un graphique de la valeur prédite par rapport à l'un des prédicteurs, puis permettra à l'utilisateur de modifier interactivement la valeur des différents prédicteurs pour voir comment cela affecte la prédiction.

Greg Snow
la source
1
Je suppose que les prédictions des «termes» sont centrées. Savez-vous comment cela se fait?
dave
4
La predict.glmfonction appelle la predict.lmfonction, qui a une section en elle que s'il y a une interception alors chaque colonne de la matrice du modèle a sa moyenne soustraite avant d'être multipliée par le vecteur de coefficient.
Greg Snow