Merencanakan interval kepercayaan untuk probabilitas yang diprediksi dari regresi logistik

Ok, saya memiliki regresi logistik dan telah menggunakan predict()fungsi untuk mengembangkan kurva probabilitas berdasarkan perkiraan saya.

## LOGIT MODEL:
library(car)
mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))

## PROBABILITY CURVE:
all.x <- expand.grid(won=unique(won), bid=unique(bid))
y.hat.new <- predict(mod1, newdata=all.x, type="response")
plot(bid<-000:1000,predict(mod1,newdata=data.frame(bid<-c(000:1000)),type="response"), lwd=5, col="blue", type="l")

Ini bagus, tetapi saya ingin tahu tentang merencanakan interval kepercayaan untuk probabilitas. Saya sudah mencoba plot.ci()tetapi tidak berhasil. Adakah yang bisa mengarahkan saya ke beberapa cara untuk menyelesaikan ini, sebaiknya dengan carpaket atau pangkalan R.

Kode yang Anda gunakan memperkirakan model regresi logistik menggunakan glmfungsi. Anda tidak memasukkan data, jadi saya hanya akan menebusnya.

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

Model regresi logistik memodelkan hubungan antara variabel respons biner dan, dalam hal ini, satu prediktor kontinu. Hasilnya adalah probabilitas logit-transformed sebagai hubungan linier dengan prediktor. Dalam kasus Anda, hasilnya adalah respons biner yang terkait dengan menang atau tidaknya menang dalam judi dan diprediksi oleh nilai taruhan. Koefisien dari mod1diberikan dalam odds log (yang sulit diinterpretasikan), sesuai dengan:

Untuk mengonversi peluang yang dicatat menjadi probabilitas, kami dapat menerjemahkan di atas ke

Anda dapat menggunakan informasi ini untuk mengatur plot. Pertama, Anda memerlukan rentang variabel prediktor:

plotdat <- data.frame(bid=(0:1000))

Kemudian menggunakan predict, Anda dapat memperoleh prediksi berdasarkan model Anda

preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)

Perhatikan bahwa nilai yang dipasang juga dapat diperoleh melalui

mod1$fitted

Dengan menentukan se.fit=TRUE, Anda juga mendapatkan kesalahan standar yang terkait dengan setiap nilai yang dipasang. Hasilnya data.frameadalah matriks dengan komponen-komponen berikut: prediksi yang dipasang ( fit), estimasi kesalahan standar ( se.fit), dan skalar yang memberikan akar kuadrat dari dispersi yang digunakan untuk menghitung kesalahan standar ( residual.scale). Dalam kasus logit binomial, nilai akan menjadi 1 (yang dapat Anda lihat dengan memasukkan preddat$residual.scaledalam R). Jika Anda ingin melihat contoh dari apa yang Anda hitung sejauh ini, Anda bisa mengetik head(data.frame(preddat)).

Langkah selanjutnya adalah mengatur plot. Saya ingin mengatur area plot kosong dengan parameter terlebih dahulu:

with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))

Sekarang Anda dapat melihat di mana penting untuk mengetahui cara menghitung probabilitas yang dipasang. Anda dapat menggambar garis yang sesuai dengan probabilitas yang dipasang mengikuti rumus kedua di atas. Menggunakan preddat data.frameAnda dapat mengonversi nilai yang dipasang ke probabilitas dan menggunakannya untuk memplot garis terhadap nilai variabel prediktor Anda.

with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))

Akhirnya, jawab pertanyaan Anda, interval kepercayaan dapat ditambahkan ke plot dengan menghitung probabilitas untuk nilai yang dipasang +/- 1.96dikali kesalahan standar:

with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

Plot yang dihasilkan (dari data yang dihasilkan secara acak) akan terlihat seperti ini:

Demi kepentingan, inilah semua kode dalam satu potongan:

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))
plotdat <- data.frame(bid=(0:1000))
preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)
with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))
with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))
with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

(Catatan: Ini adalah jawaban yang diedit dalam upaya untuk membuatnya lebih relevan dengan stats.stackexchange.)

Ini adalah modifikasi dari solusi @ smillig. Saya menggunakan alat rapi di sini, dan juga menggunakan linkinvfungsi yang merupakan bagian dari objek model GLM mod1. Dengan begitu, Anda tidak perlu membalikkan fungsi logistik secara manual, dan pendekatan ini akan berfungsi terlepas dari GLM spesifik apa yang Anda muat.

library(tidyverse)
library(magrittr)


set.seed(1234)

# create fake data on gambling. Does prob win depend on bid size? 
mydat <- data.frame(
  won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
  bid=runif(250, min=0, max=1000)
)

# logistic regression model: 
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

# new predictor values to use for prediction: 
plotdat <- data.frame(bid=(0:1000))

# df with predictions, lower and upper limits of CIs: 
preddat <- predict(mod1,
               type = "link",
               newdata=plotdat,
               se.fit=TRUE) %>% 
  as.data.frame() %>% 
  mutate(bid = (0:1000), 

         # model object mod1 has a component called linkinv that 
         # is a function that inverts the link function of the GLM:
         lower = mod1$family$linkinv(fit - 1.96*se.fit), 
         point.estimate = mod1$family$linkinv(fit), 
         upper = mod1$family$linkinv(fit + 1.96*se.fit)) 


# plotting with ggplot: 
preddat %>% ggplot(aes(x = bid, 
                   y = point.estimate)) + 
  geom_line(colour = "blue") + 
  geom_ribbon(aes(ymin = lower,
                  ymax = upper), 
              alpha = 0.5) + 
  scale_y_continuous(limits = c(0,1))