CrossValidated has several questions on when and how to apply the rare event bias correction by King and Zeng (2001). I am looking for something different: a minimal simulation-based demonstration that the bias exists.
In particular, King and Zeng state
"... in rare events data the biases in probabilities can be substantively meaningful with sample sizes in the thousands and are in a predictable direction: estimated event probabilities are too small."
Here is my attempt to simulate such a bias in 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)
When I run this, I tend to get very small z-scores, and the histogram of estimates is very close to centered over the truth p = 0.01.
What am I missing? Is it that my simulation is not large enough show the true (and evidently very small) bias? Does the bias require some kind of covariate (more than the intercept) to be included?
Update 1: King and Zeng include a rough approximation for the bias of $\beta_0$ in equation 12 of their paper. Noting the N in the denominator, I drastically reduced N to be 5 and re-ran the simulation, but still no bias in the estimated event probabilities is evident. (I used this only as inspiration. Note that my question above is about estimated event probabilities, not $\hat \beta_0$.)
Update 2: Following a suggestion in the comments, I included an independent variable in the regression, leading to equivalent results:
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))
Explanation: I used p itself as the independent variable, where p is a vector with a repetitions of a small value (0.01) and a larger value (0.2). In the end, sim stores only the estimated probabilities corresponding to $p = 0.01$ and there is no sign of bias.
Update 3 (May 5, 2016): This does not noticeably change the results, but my new inner simulation function is
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])
}
}
Explanation: The MLE when y is identically zero does not exist (thanks to comments here for the reminder). R fails to throw a warning because its "positive convergence tolerance" actually gets satisfied. More liberally speaking, the MLE exists and is minus infinity, which corresponds to $p=0$; hence my function update. The only other coherent thing I can think of doing is discarding those runs of the simulation where y is identically zero, but that would clearly lead to results even more counter to the initial claim that "estimated event probabilities are too small".
