罕見事件邏輯回歸偏差：如何用最小示例模擬被低估的 p？
CrossValidated 有幾個關於何時以及如何應用King 和 Zeng (2001)的罕見事件偏差校正的問題。我正在尋找不同的東西：基於最小模擬的證明存在偏差。
特別是國王和曾國
“……在罕見事件數據中，概率偏差在數千個樣本量下可能具有實質性意義，並且處於可預測的方向：估計的事件概率太小。”
這是我在 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)當我運行這個時，我傾向於得到非常小的 z 分數，並且估計的直方圖非常接近於以真值 p = 0.01 為中心。
我錯過了什麼？是不是我的模擬不夠大，無法顯示真實的（顯然非常小的）偏差？偏差是否需要包含某種協變量（超過截距）？
更新 1： King 和 Zeng 對在他們論文的等式12中。注意到
N分母中的 ，我大幅減少並N重新5運行了模擬，但估計的事件概率仍然沒有明顯的偏差。（我僅將此用作靈感。請注意，我上面的問題是關於估計的事件概率，而不是.)**更新 2：**根據評論中的建議，我在回歸中包含了一個自變量，得到了相同的結果：
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))說明：我用
p自己作為自變量，這裡p是一個向量，有一個小值（0.01）和一個大值（0.2）的重複。最後，sim只存儲對應的估計概率並且沒有偏見的跡象。更新 3（2016 年 5 月 5 日）： 這並沒有明顯改變結果，但我的新內部模擬函數是
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]) } }說明：當 y 為零時的 MLE 不存在（感謝此處的評論提醒）。R 未能發出警告，因為它的“正收斂容差”實際上得到了滿足。更自由地說，MLE 存在並且是負無窮大，對應於; 因此我的功能更新。我能想到的唯一其他連貫的事情是丟棄那些 y 相同為零的模擬運行，但這顯然會導致結果與最初聲稱的“估計的事件概率太小”更加相反。
這是一個有趣的問題——我做了一些模擬，我在下面發布，希望這能激發進一步的討論。
首先，一些一般性評論：
- 你引用的論文是關於罕見事件偏差的。我之前不清楚（也關於上面的評論）是如果你有 10/10000 而不是 10/30 觀察的情況有什麼特別之處。但是，經過一些模擬，我同意存在。
- 我想到的一個問題（我經常遇到這種情況，最近在生態學和進化的方法中有一篇論文，但我找不到參考資料）是你可以在小數據中得到 GLM 的退化案例在這種情況下，MLE 與事實相距甚遠，甚至在 - / + 無窮大（我想是由於非線性鏈接）。我不清楚在偏差估計中應該如何處理這些情況，但從我的模擬來看，我會說它們似乎是罕見事件偏差的關鍵。我的直覺是刪除它們，但是還不清楚它們必須被刪除多遠。對於偏差校正，也許要記住一些事情。
- 此外，這些退化的情況似乎容易導致數值問題（因此我增加了 glm 函數中的 maxit，但人們也可以考慮增加 epsilon 以確保實際報告真正的 MLE）。
無論如何，這裡有一些代碼計算邏輯回歸中截距、斜率和預測的估計值和真值之間的差異，首先是針對低樣本量/中等發生率的情況：
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)截距、斜率和預測的偏差和標準誤差為
-0.120429315 0.296453122 -0.001619793 0.016105833 0.032835468 0.002040664我會得出結論，有很好的證據表明截距有輕微的負偏差，斜率有輕微的正偏差，儘管看一下繪製的結果表明，與估計值的方差相比，偏差很小。
如果我將參數設置為罕見情況
N= 4000 slope = 2 # slope (linear scale) intercept = - 10 # intercept (linear scale)我對截距有更大的偏差，但預測仍然沒有
-1.716144e+01 4.271145e-01 -3.793141e-06 5.039331e-01 4.806615e-01 4.356062e-06在估計值的直方圖中，我們看到了退化參數估計的現象（如果我們應該這樣稱呼的話）
讓我們刪除所有截距估計值 <20 的行
apply(bias[bias[,1] > -20,], 2, mean) apply(bias[bias[,1] > -20,], 2, sd) / sqrt(length(bias[,1] > -10))偏差減少了，圖中的情況變得更加清晰——參數估計顯然不是正態分佈的。我想知道這意味著報告的 CI 的有效性。
-0.6694874106 1.9740437782 0.0002079945 1.329322e-01 1.619451e-01 3.242677e-06
我會得出結論，截距上的罕見事件偏差是由罕見事件本身驅動的，即那些罕見的、極小的估計。不確定我們是否要刪除它們，不確定截止日期是什麼。
但需要注意的重要一點是，無論哪種方式，響應規模的預測似乎都沒有偏差——鏈接函數只是吸收了這些極小的值。
