---
title: "What model to use for rare events"
date: "2017-10-26"
---
## Introduction
In empirical studies, people are worried about rare event situation. That is, when you have, for example, lots of 0's and only a few 1's, or vice versa. Do you run a logit model, or do you use a "rare event logit"? When should you use either approach? Or there is a third approach?
Paul Allison said in his blog (https://statisticalhorizons.com/logistic-regression-for-rare-events):
"Prompted by a 2001 article by King and Zeng, many researchers worry about whether they can legitimately use conventional logistic regression for data in which events are rare. Although King and Zeng accurately described the problem and proposed an appropriate solution, there are still a lot of misconceptions about this issue.
The problem is not specifically the rarity of events, but rather the possibility of a small number of cases on the rarer of the two outcomes. If you have a sample size of 1000 but only 20 events, you have a problem. If you have a sample size of 10,000 with 200 events, you may be OK. If your sample has 100,000 cases with 2000 events, you're golden."
In general I agree with him. However, when exactly should we use King and Zeng's "relogit"?
Allison also mentioned two other methods. One is called the Firth method, a penalized likelihood approach. The other one is the exact logistic regression, which is for small samples.
In this simulation exercise, I compare logit against bias-reduced logistic
regression (`brglm2` with `method="brglmFit"`, `type="AS_mean"`) and the Firth
model (`logistf`).
One thing to be clear about before reading any comparison of the latter two: for a
logistic regression **they are the same estimator.** Firth's penalized likelihood
(penalizing by the Jeffreys prior) coincides with mean bias reduction for GLMs
with a canonical link, and logit is canonical for the binomial. So `brglm2`'s
`AS_mean` and `logistf` solve the same estimating equations and differ only in
numerical tolerance and in how each package builds standard errors and intervals
(`logistf` uses penalized profile likelihood). A demonstration, on a
deliberately rare-event DGP: $n = 500$, $x \sim N(0,1)$, and
$y \sim \text{Bernoulli}(\Lambda(-3.5 + 0.8x))$, which puts the event rate
around 3% — roughly fifteen events in the whole sample.
```{r}
#| label: firth-equivalence
#| message: false
#| warning: false
library(brglm2); library(logistf)
set.seed(11)
n <- 500; x <- rnorm(n); y <- rbinom(n, 1, plogis(-3.5 + 0.8 * x))
d <- data.frame(y, x)
rbind(
`plain logit` = coef(glm(y ~ x, data = d, family = binomial)),
`brglm2 AS_mean` = coef(glm(y ~ x, data = d, family = binomial,
method = "brglmFit", type = "AS_mean")),
`logistf (Firth)` = coef(logistf(y ~ x, data = d))
)
```
The two bias-reduced fits give $-3.134453$ and $0.7277674$ against
$-3.134453$ and $0.7277675$: they agree to about $10^{-7}$, which is numerical
tolerance rather than any statistical difference. Plain logit gives $-3.176$ and
$0.738$, so on this draw the bias correction moves the slope by about 0.010 and
the intercept by 0.042 — both estimators sitting below the true $0.8$ and above
the true $-3.5$, which is what fifteen events buys you. This matters for reading the findings below: what
the simulation can tell us is how **logit** compares with **Firth-type bias
reduction**, not which of two implementations of the same estimator is better. If
a genuinely distinct third method is wanted, `brglm2`'s `type = "AS_median"`
(median bias reduction) is one, and King and Zeng's `relogit` correction is
another.
## Simulation
Here I have some code for using multiple cores to run these three models. The bias-reduced logistic regression is implemented in R via the `brglm2` package, the Firth method is implemented in `logistf`.
```{r}
#| label: rare-events-sim
#| eval: false
library(brglm2)
library(logistf)
require(snowfall)
set.seed(666)
# initialize parallel cores.
sfInit(parallel=TRUE, cpus=16)
gen.sim <- function(df){
x <- rnorm(df['nobs'], 0, 1)
# generate binary data
p <- df['p']
# NB: alpha = qlogis(p) makes p the event probability at x = 0, not the
# marginal event probability. Because the slope on x is 2 and x ~ N(0,1),
# the marginal rate is much higher than p: integrating plogis(alpha + 2x)
# over x gives 0.044, 0.134 and 0.203 for p = 0.01, 0.05 and 0.10.
alpha <- -log((1-p)/p)
z = alpha + 2*x
pr = 1/(1+exp(-z))
y = rbinom(df['nobs'], 1, pr)
df = data.frame(y=y, x=x)
# With small nobs and small p (e.g. nobs=10, p=.01), y can easily come out
# all-zero (>90% chance in that cell), which makes glm/logistf unable to
# estimate a coefficient on x. Guard against that rather than letting a
# single degenerate draw crash the whole parallel simulation loop.
if (var(y) == 0) {
return(c(logit=NA, brglm2=NA, logistf=NA))
}
# logit
m1 <- tryCatch(glm(y ~ x, family='binomial'), error = function(e) NULL)
m1.x <- if (!is.null(m1)) summary(m1)$coefficients['x','Estimate'] - 2 else NA
# bias-reduced logistic regression (brglm2)
m2 <- tryCatch(
glm(y ~ x, family = binomial(link = "logit"), data = df,
method = "brglmFit", type = "AS_mean"),
error = function(e) NULL
)
m2.x <- if (!is.null(m2)) coef(m2)['x'] - 2 else NA
# logistf (Firth)
m3 <- tryCatch(logistf(y ~ x, data=df), error = function(e) NULL)
m3.x <- if (!is.null(m3)) coef(m3)['x'] - 2 else NA
return(c(logit=m1.x, brglm2=m2.x, logistf=m3.x))
}
# set parameter space
sim.grid = seq(1, 100, 1)
p.grid = c(.01, .05, .1)
nobs.grid = c(10, 30, 50, 100, 200, 500, 1000, 10000)
data.grid <- expand.grid(nobs.grid, sim.grid, p.grid)
names(data.grid) <- c('nobs', 'nsim', 'p')
# export functions and libraries to parallel workers
sfExport(list=list("gen.sim"))
sfLibrary(brglm2)
sfLibrary(logistf)
results <- data.frame(t(sfApply(data.grid, 1, gen.sim)))
# stop the cluster
sfStop()
forshiny <- cbind(data.grid, results)
write.csv(forshiny, 'results.csv')
```
We simulate 100 times with sample size from 10 to 10000, event probability .01, .05, and .1.
Since there are many simulations, we used the `snowfall` library to speed things up.
(The original post plotted bias and MSE by sample size and event probability from the `results.csv` output above; those figures are not reproduced here — the findings are summarized in prose below.)
A caveat on how to read the grid before the findings. The parameter `p` is the
event probability at $x=0$, not the marginal event probability: with a slope of 2
on $x \sim N(0,1)$, the marginal rate is roughly 0.044, 0.134 and 0.203 for
`p` = 0.01, 0.05 and 0.10. So the expected number of events in a cell is closer
to $n$ times those numbers than to $n p$, by a factor of about four at the rare
end. What follows is stated in terms of the realized event count, which is the
quantity that actually governs the behaviour.
When the rarer of the two outcomes has only a handful of observations -- fewer
than about five events in the sample -- none of these three models performs well.
This is understandable: there is almost nothing in the data to identify the slope.
Once the rarer group has more than about fifty events, there is not much
difference between the three. In between, we found that Firth-type bias reduction
--- whether computed by `brglm2` or by `logistf` --- performs better than plain
logit. Any apparent ranking *between* those two should be discounted, since as shown above they are the same estimator; differences at that level reflect numerical tolerance or differing standard-error conventions, not a substantive advantage.
In the small sample situation, maybe it's better to use the exact logistic regression.