Say, for example, that we have a logistic regression model which outputs the probability that a patient will develop a particular disease based on many covariates.
We can get an idea of the magnitude and direction of the effect of each covariate in general by examining the coefficients of the model and considering the change in the odds-ratio.
What if we want to know for a single patient what his or her biggest risk factors / the biggest factors in his or her favor. I'm particularly interested in those which the patient could actually do something about.
What is the best way to do this?
The way I'm currently considering is captured in the following R code (taken from this thread):
#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')])
I'm considering looking additionally at
this.student.prediction.list <- this.student.predictors * coef(data.model)
and trying to get the information out of the individual addends of the sum which is the probability estimate, but I'm not sure how to do it.
I could look at
- Which variables make the largest absolute contribution to the probability estimate and take those to be the largest risk factors.
- Which variables differ by the greatest amount from their mean proportion, i.e. see what proportion each variable contributes to the probability estimate on average and see which variables differ from this proportion by the greatest amount in this particular observation
- A combination thereof: weight the absolute difference between mean proportion and observed proportion by the mean proportion and take those variables with the greatest weighted values
Which of these make the most sense? Would any of these approaches be a reasonable way to answer the question?
Additionally, I'd like to know how I could get confidence intervals for the additive contributions of individual covariates to the probability estimate.
