Kontribusi setiap kovariat untuk prediksi tunggal dalam model regresi logistik

8

Katakanlah, misalnya, bahwa kita memiliki model regresi logistik yang menghasilkan probabilitas bahwa seorang pasien akan mengembangkan penyakit tertentu berdasarkan banyak kovariat.

Kita bisa mendapatkan gambaran tentang besarnya dan arah pengaruh masing-masing kovariat secara umum dengan memeriksa koefisien model dan mempertimbangkan perubahan dalam odds-rasio.

Bagaimana jika kita ingin tahu untuk seorang pasien tunggal apa faktor risiko terbesarnya / faktor terbesar yang menguntungkannya. Saya sangat tertarik pada hal-hal yang sebenarnya bisa dilakukan oleh pasien.

Apa cara terbaik untuk melakukan ini?

Cara saya saat ini mempertimbangkan ditangkap dalam kode R berikut (diambil dari utas ini ):

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

Saya sedang mempertimbangkan melihat tambahan di

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

dan mencoba untuk mendapatkan informasi dari masing-masing penambahan jumlah yang merupakan perkiraan probabilitas, tapi saya tidak yakin bagaimana melakukannya.

Saya bisa melihat

  • Variabel mana yang memberikan kontribusi absolut terbesar pada estimasi probabilitas dan menjadikannya sebagai faktor risiko terbesar.
  • Variabel mana yang berbeda dengan jumlah terbesar dari proporsi rata-rata mereka, yaitu melihat proporsi masing-masing variabel yang berkontribusi terhadap estimasi probabilitas rata-rata dan melihat variabel mana yang berbeda dari proporsi ini dengan jumlah terbesar dalam pengamatan khusus ini.
  • Kombinasi keduanya: bobot perbedaan absolut antara proporsi rata-rata dan proporsi yang diamati dengan proporsi rata-rata dan ambil variabel-variabel tersebut dengan nilai tertimbang terbesar.

Manakah dari ini yang paling masuk akal? Akankah salah satu dari pendekatan ini menjadi cara yang masuk akal untuk menjawab pertanyaan?

Selain itu, saya ingin tahu bagaimana saya bisa mendapatkan interval kepercayaan untuk kontribusi aditif kovariat individu ke estimasi probabilitas.

dave
sumber

Jawaban:

10

Anda dapat menggunakan predictfungsi dalam R. Sebut dengan type='terms'dan itu akan memberi Anda kontribusi setiap istilah dalam model (koefisien dikalikan nilai variabel). Ini akan berada pada skala log-odds.

Pilihan lain adalah menggunakan TkPredictfungsi dari paket TeachingDemos. Ini akan menampilkan grafik dari nilai yang diprediksi vs salah satu dari prediktor, lalu biarkan pengguna secara interaktif mengubah nilai dari berbagai prediktor untuk melihat bagaimana hal itu mempengaruhi prediksi.

Greg Snow
sumber
1
Prediksi 'istilah', saya kumpulkan, berpusat. Apakah Anda tahu bagaimana ini dilakukan?
dave
4
The predict.glmfungsi memanggil predict.lmfungsi, yang memiliki bagian di dalamnya yang jika ada intercept maka setiap kolom dari matriks model memiliki mean dikurangkan dari itu sebelum dikalikan dengan vektor koefisien.
Greg Snow