Memperkirakan rasio risiko yang disesuaikan dalam data biner menggunakan regresi Poisson

9

Saya tertarik untuk memperkirakan rasio risiko yang disesuaikan, analog dengan bagaimana seseorang memperkirakan rasio odds yang disesuaikan menggunakan regresi logistik. Beberapa literatur (misalnya, ini ) menunjukkan bahwa menggunakan regresi Poisson dengan kesalahan standar Huber-White adalah cara berbasis model untuk melakukan ini

Saya belum menemukan literatur tentang bagaimana penyesuaian untuk kovariat berkelanjutan mempengaruhi hal ini. Simulasi sederhana berikut menunjukkan bahwa masalah ini tidak mudah:

arr <- function(BLR,RR,p,n,nr,ce)
{
   B = rep(0,nr)
   for(i in 1:nr){
   b <- runif(n)<p 
   x <- rnorm(n)
   pr <- exp( log(BLR) + log(RR)*b + ce*x)
   y <- runif(n)<pr
   model <- glm(y ~ b + x, family=poisson)
   B[i] <- coef(model)[2]
   }
   return( mean( exp(B), na.rm=TRUE )  )
}

set.seed(1234)
arr(.3, 2, .5, 200, 100, 0)
[1] 1.992103
arr(.3, 2, .5, 200, 100, .1)
[1] 1.980366
arr(.3, 2, .5, 200, 100, 1)
[1] 1.566326

Dalam hal ini, rasio risiko sebenarnya adalah 2, yang pulih dengan andal ketika efek kovariat kecil. Tetapi, ketika efek kovariat besar, ini terdistorsi. Saya berasumsi ini muncul karena efek kovariat dapat mendorong ke atas terhadap batas atas (1) dan ini mencemari estimasi.

Saya telah melihat tetapi belum menemukan literatur tentang penyesuaian untuk kovariat berkesinambungan dalam estimasi rasio risiko yang disesuaikan. Saya mengetahui posting berikut di situs ini:

tetapi mereka tidak menjawab pertanyaan saya. Apakah ada makalah tentang ini? Apakah ada peringatan yang diketahui yang harus dilakukan?

r references biostatistics poisson-regression relative-risk kjetil b halvorsen
sumber
1
Mungkin menarik bagi Anda: aje.oxfordjournals.org/content/162/3/199.full
StatsStudent
Juga, T&J ini stats.stackexchange.com/questions/18595/… dapat membantu.
mdewey

Jawaban:

1

Saya tidak tahu apakah Anda masih membutuhkan jawaban untuk pertanyaan ini, tetapi saya memiliki masalah serupa di mana saya ingin menggunakan regresi Poisson. Dalam menjalankan kode Anda, saya menemukan bahwa jika saya mengatur model sebagai

model <- glm(y ~ b + x, family=binomial(logit)

alih-alih sebagai model regresi Poisson Anda, hasil yang sama terjadi: perkiraan OR adalah ~ 1,5 ketika ce mendekati 1. Jadi, saya tidak yakin bahwa contoh Anda memberikan informasi tentang kemungkinan masalah dengan penggunaan regresi Poisson untuk hasil biner.

David F
sumber
1
Masalah dengan pemasangan model logit, meskipun tidak mengarah pada risiko yang diprediksi lebih besar dari 1, adalah bahwa rasio odds adalah penaksir yang bias dari rasio risiko dan bahwa bias meningkat secara dramatis ketika hasilnya menjadi lebih lazim. Anda dapat menentukan binomial(link=log)untuk benar-benar cocok dengan model risiko relatif, tetapi jarang konvergen karena hasil yang terlalu tinggi.
AdamO
1

Saya menemukan bahwa menggunakan kemungkinan maksimum langsung dengan fungsi probabilitas yang tepat sangat meningkatkan estimasi risiko relatif. Anda dapat secara langsung menentukan fungsi risiko terpotong sebagai tingkat prediksi untuk proses tersebut.

masukkan deskripsi gambar di sini

Biasanya kami menggunakan Hessian untuk membuat CI untuk estimasi. Saya belum mengeksplorasi kemungkinan menggunakannya sebagai matriks "B" (daging) dalam kesalahan Huber White dan menggunakan risiko yang dipasang untuk mendapatkan matriks "A" (roti) ... tapi saya menduga itu bisa bekerja! Lebih layaknya Anda dapat menggunakan bootstrap untuk mendapatkan kesalahan model yang kuat untuk hubungan mean-variance yang salah ditentukan.

## the negative log likelihood for truncated risk function
negLogLik <- function(best, X, y) { 
  pest <- pmin(1, exp(X %*% best))
  -sum(dpois(x = y, lambda = pest, log=TRUE))
}

set.seed(100)

sim <- replicate(100, {
  n <- 200
  X <- cbind(1, 'b'=rbinom(n, 1, 0.5), 'x'=rnorm(n))
  btrue <- c(log(0.3), log(2), 1)
  ptrue <- pmin(1, exp(X %*% matrix(btrue)))
  y <- rbinom(n, 1, ptrue) ## or just take y=ptrue for immediate results
  nlm(f = logLik, p = c(log(mean(y)),0,0), X=X, y=y)$estimate
})

rowMeans(exp(sim))

Memberi:

> rowMeans(exp(sim))
[1] 0.3002813 2.0680780 3.0888280

Koefisien tengah memberi Anda apa yang Anda inginkan.

AdamO
sumber