8
$\begingroup$

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.

$\endgroup$
10
+50
$\begingroup$

You can use the predict function in R. Call it with type='terms' and it will give you the contribution of each term in the model (the coefficient times the variable value). This will be on the log-odds scale.

Another option is to use the TkPredict function from the TeachingDemos package. This will show a graph of the predicted value vs. one of the predictors, then let the user interactively change the value of the various predictors to see how that affects the prediction.

$\endgroup$
  • 1
    $\begingroup$ The 'terms' predictions, I gather, are centered. Do you know how this is done? $\endgroup$ – dave Sep 7 '12 at 2:25
  • 4
    $\begingroup$ The predict.glm function calls the predict.lm function, which has a section in it that if there is an intercept then each column of the model matrix has its mean subtracted from it before being multiplied by the coefficient vector. $\endgroup$ – Greg Snow Sep 7 '12 at 15:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.