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 中注意到分母中的 ，我大幅减少并重新运行了模拟，但估计的事件概率仍然没有明显的偏差。（我仅将此用作灵感。请注意，我上面的问题是关于估计的事件概率，而不是。）NN5
更新 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 相同为零的模拟运行，但这显然会导致结果与最初声称的“估计的事件概率太小”更加相反。
