Peristiwa regresi logistik yang langka: bagaimana mensimulasikan p yang diremehkan dengan contoh minimal?

CrossValidated memiliki beberapa pertanyaan tentang kapan dan bagaimana menerapkan koreksi bias peristiwa langka oleh King dan Zeng (2001) . Saya mencari sesuatu yang berbeda: demonstrasi berbasis simulasi minimal bahwa bias ada.

Secara khusus, Raja dan negara Zeng

"... dalam data kejadian yang jarang, bias dalam probabilitas dapat secara substantif bermakna dengan ukuran sampel dalam ribuan dan berada dalam arah yang dapat diprediksi: estimasi probabilitas kejadian terlalu kecil."

Ini adalah usaha saya untuk mensimulasikan bias pada R:

# FUNCTIONS
do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)  
    # Extract the fitted probability.
    #    If p is constant, glm does y ~ 1, the intercept-only model.
    #    If p is not constant, assume its smallest value is p[1]:
    glm(y ~ p, family = 'binomial')$fitted[1]
}
mean.of.K.estimates = function(p, K){
    mean(replicate(K, do.one.sim(p) ))
}

# MONTE CARLO
N = 100
p = rep(0.01, N)
reps = 100
# The following line may take about 30 seconds
sim = replicate(reps, mean.of.K.estimates(p, K=100))
# Z-score:
abs(p[1]-mean(sim))/(sd(sim)/sqrt(reps))
# Distribution of average probability estimates:
hist(sim)

Ketika saya menjalankan ini, saya cenderung mendapatkan skor-z yang sangat kecil, dan histogram perkiraan sangat dekat dengan terpusat pada kebenaran p = 0,01.

Apa yang saya lewatkan? Apakah simulasi saya tidak cukup besar menunjukkan bias yang benar (dan ternyata sangat kecil)? Apakah bias membutuhkan semacam kovariat (lebih dari intersepsi) untuk dimasukkan?

Pembaruan 1: King dan Zeng menyertakan perkiraan kasar untuk bias dalam persamaan 12 dari makalah mereka. Memperhatikan penyebutnya, saya secara drastis mengurangi menjadi dan menjalankan kembali simulasi, tetapi masih tidak ada bias dalam probabilitas kejadian yang diperkirakan. (Saya menggunakan ini hanya sebagai inspirasi. Perhatikan bahwa pertanyaan saya di atas adalah tentang perkiraan probabilitas acara, bukan .)$\beta_0$NN5$\hat \beta_0$

Pembaruan 2: Mengikuti saran dalam komentar, saya menyertakan variabel independen dalam regresi, yang mengarah ke hasil yang setara:

p.small = 0.01
p.large = 0.2
p = c(rep(p.small, round(N/2) ), rep(p.large, N- round(N/2) ) )
sim = replicate(reps, mean.of.K.estimates(p, K=100))

Penjelasan: Saya menggunakan pdirinya sebagai variabel independen, di mana padalah vektor dengan pengulangan nilai kecil (0,01) dan nilai lebih besar (0,2). Pada akhirnya, simmenyimpan hanya probabilitas yang diperkirakan sesuai dengan dan tidak ada tanda bias.$p = 0.01$

Pembaruan 3 (5 Mei 2016): Ini tidak terasa mengubah hasil, tetapi fungsi simulasi batin baru saya adalah

do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(0)
    }else{
        # Extract the fitted probability.
        #    If p is constant, glm does y ~ 1, the intercept only model.
        #    If p is not constant, assume its smallest value is p[1]:
        return(glm(y ~ p, family = 'binomial')$fitted[1])
    }
}

Penjelasan: MLE ketika y secara identik nol tidak ada ( terima kasih untuk komentar di sini untuk pengingat ). R gagal memberikan peringatan karena " toleransi konvergensi positif " -nya benar-benar puas. Lebih bebas berbicara, MLE ada dan minus tanpa batas, yang sesuai dengan ; karenanya pembaruan fungsi saya. Satu-satunya hal yang koheren lain yang dapat saya pikirkan lakukan adalah membuang menjalankan simulasi di mana y adalah identik nol, tapi itu jelas akan mengarah pada hasil yang lebih bertentangan dengan klaim awal bahwa "probabilitas kejadian diperkirakan terlalu kecil".$p=0$

Ini adalah pertanyaan yang menarik - saya telah melakukan beberapa simulasi yang saya posting di bawah ini dengan harapan ini merangsang diskusi lebih lanjut.

Pertama-tama, beberapa komentar umum:

• Makalah yang Anda kutip adalah tentang bias peristiwa langka. Apa yang tidak jelas bagi saya sebelumnya (juga sehubungan dengan komentar yang dibuat di atas) adalah jika ada sesuatu yang istimewa tentang kasus di mana Anda memiliki 10/10000 yang bertentangan dengan pengamatan 10/30. Namun, setelah beberapa simulasi, saya setuju ada.

• Masalah yang ada dalam pikiran saya (saya sering menemukan ini, dan baru-baru ini ada sebuah makalah dalam Metode dalam Ekologi dan Evolusi tentang itu, saya tidak dapat menemukan rujukannya) adalah bahwa Anda bisa mendapatkan kasus degenerasi dengan GLM dalam data kecil situasi, di mana MLE adalah FAAAR jauh dari kebenaran, atau bahkan pada - / + infinity (karena tautan nonlinier saya kira). Tidak jelas bagi saya bagaimana seseorang harus memperlakukan kasus-kasus ini dalam estimasi bias, tetapi dari simulasi saya, saya akan mengatakan mereka tampaknya kunci untuk bias kejadian langka. Intuisi saya adalah menghapusnya, tetapi tidak jelas seberapa jauh mereka harus dihapus. Mungkin ada sesuatu yang perlu diingat untuk koreksi-bias.

• Juga, kasus-kasus degenerasi ini tampaknya cenderung menyebabkan masalah numerik (oleh karena itu saya telah meningkatkan maxit dalam fungsi glm, tetapi orang dapat berpikir tentang peningkatan epsilon juga untuk memastikan seseorang benar-benar melaporkan MLE yang sebenarnya).

Bagaimanapun, di sini beberapa kode yang menghitung perbedaan antara perkiraan dan kebenaran untuk intersep, kemiringan dan prediksi dalam regresi logistik, pertama untuk ukuran sampel rendah / situasi kejadian sedang:

set.seed(123)
replicates = 1000
N= 40
slope = 2 # slope (linear scale)
intercept = - 1 # intercept (linear scale)

bias <- matrix(NA, nrow = replicates, ncol = 3)
incidencePredBias <- rep(NA, replicates)

for (i in 1:replicates){
  pred = runif(N,min=-1,max=1) 
  linearResponse = intercept + slope*pred
  data = rbinom(N, 1, plogis(linearResponse))  
  fit <- glm(data ~ pred, family = 'binomial', control = list(maxit = 300))
  bias[i,1:2] = fit$coefficients - c(intercept, slope)
  bias[i,3] = mean(predict(fit,type = "response")) - mean(plogis(linearResponse))
}

par(mfrow = c(1,3))
text = c("Bias intercept", "Bias slope", "Bias prediction")

for (i in 1:3){
  hist(bias[,i], breaks = 100, main = text[i])
  abline(v=mean(bias[,i]), col = "red", lwd = 3)  
}

apply(bias, 2, mean)
apply(bias, 2, sd) / sqrt(replicates)

Bias yang dihasilkan dan kesalahan standar untuk intersep, kemiringan dan prediksi adalah

-0.120429315  0.296453122 -0.001619793
 0.016105833  0.032835468  0.002040664

Saya akan menyimpulkan bahwa ada bukti yang cukup bagus untuk sedikit bias negatif dalam intersepsi, dan sedikit bias positif di lereng, meskipun melihat hasil yang diplot menunjukkan bahwa bias tersebut kecil dibandingkan dengan varians dari nilai estimasi.

Jika saya mengatur parameter ke situasi peristiwa langka

N= 4000
slope = 2 # slope (linear scale)
intercept = - 10 # intercept (linear scale)

Saya mendapatkan bias yang lebih besar untuk intersep, tetapi masih TIDAK ada prediksi

   -1.716144e+01  4.271145e-01 -3.793141e-06
    5.039331e-01  4.806615e-01  4.356062e-06

Dalam histogram nilai estimasi, kita melihat fenomena perkiraan parameter degenerasi (jika kita harus memanggilnya seperti itu)

Mari kita hapus semua baris yang perkiraan intersepnya <20

apply(bias[bias[,1] > -20,], 2, mean)
apply(bias[bias[,1] > -20,], 2, sd) / sqrt(length(bias[,1] > -10))

Bias menurun, dan hal-hal menjadi sedikit lebih jelas dalam gambar - estimasi parameter jelas tidak terdistribusi normal. Saya ingin tahu bahwa itu berarti validitas CI yang dilaporkan.

-0.6694874106  1.9740437782  0.0002079945
1.329322e-01 1.619451e-01 3.242677e-06

Saya akan menyimpulkan bias peristiwa langka pada intersep didorong oleh peristiwa langka itu sendiri, yaitu perkiraan langka, sangat kecil. Tidak yakin apakah kami ingin menghapusnya atau tidak, tidak yakin akan seperti apa cutoffnya.

Satu hal yang penting untuk dicatat adalah bahwa, bagaimanapun juga, tampaknya tidak ada bias pada prediksi pada skala respons - fungsi tautan hanya menyerap nilai-nilai yang sangat kecil ini.

Peristiwa langka bias hanya terjadi ketika ada regressor. Itu tidak akan terjadi dalam model intersep-only seperti yang disimulasikan di sini. Lihat posting ini untuk detail: http://statisticalhorizons.com/linear-vs-logistic#comment-276108

Gambar 7 dalam makalah ini tampaknya paling langsung menjawab pertanyaan tentang bias dalam prediksi. Saya tidak sepenuhnya memahami angka tersebut (khususnya, interpretasi "perkiraan probabilitas kejadian terlalu kecil" tampaknya seperti penyederhanaan yang berlebihan) tetapi saya berhasil mereproduksi sesuatu yang mirip dengan itu berdasarkan deskripsi singkat mereka tentang simulasi mereka di Bagian 6.1:

n_grid = 40
x_grid = seq(0, 7, length.out = n_grid)
beta0 = -6
beta1 = 1

inverse_logit = function(x) 1/(1 + exp(-x))

do.one.sim = function(){
    N = 5000
    x = rnorm(N)
    p = inverse_logit(beta0 + beta1*x)
    # Draw fake data based on probabilities p
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(rep(0, n_grid))
    }else{
        # Extract the error
        mod = glm(y ~ x, family = 'binomial')
        truth = inverse_logit(beta0 + beta1*x_grid)
        pred = predict(mod, newdata = data.frame(x = x_grid),
            type = 'response')
        return(pred - truth)
    }
}
mean.of.K.estimates = function(K){
    rowMeans(replicate(K, do.one.sim()))
}

set.seed(1)
bias = replicate(10, mean.of.K.estimates(100))
maxes = as.numeric(apply(bias, 1, max))
mins = as.numeric(apply(bias, 1, min))

par(mfrow = c(3, 1), mar = c(4,4,2,2))
plot(x_grid, rowMeans(bias), type = 'l',
    ylim = c(min(bias), max(bias)),
    xlab = 'x', ylab = 'bias')
lines(x_grid, maxes, lty = 2)
lines(x_grid, mins, lty = 2)
plot(x_grid, dnorm(x_grid), type = 'l',
    xlab = 'x', ylab = 'standard normal density')
plot(x_grid, inverse_logit(beta0 + beta1*x_grid),
    xlab = 'x', ylab = 'true simulation P(Y = 1)',
    type = 'l')

Plot pertama adalah replikasi saya dari angka 7 mereka, dengan penambahan kurva putus-putus mewakili berbagai hasil selama 10 percobaan.

Sesuai makalah, xberikut adalah variabel prediktor dalam regresi yang diambil dari standar normal. Jadi, seperti yang diilustrasikan dalam plot kedua, frekuensi relatif pengamatan untuk x > 3(di mana bias paling terlihat terjadi di plot pertama) semakin kecil.

Plot ketiga menunjukkan probabilitas simulasi "benar" dalam proses menghasilkan sebagai fungsi dari x. Tampaknya bias terbesar terjadi di mana xlangka atau tidak ada.

Secara keseluruhan, ini menunjukkan bahwa OP sepenuhnya salah mengartikan klaim sentral dari makalah dengan membingungkan "peristiwa langka" (yaitu x > 3) dengan "acara yang P(Y = 1)sangat kecil". Agaknya makalah ini menyangkut yang pertama, bukan yang terakhir.