The King and Zeng Estimators

King and Zeng’s article initially focuses on data gathering and notes that, for rare events, substantial cost savings can be had by undersampling the non-events. So long as a weighted version of the logistic estimator is used (with weights based on the proportion of events in the sample and in the population), unbiased parameter estimates can still be had. Even if undersampling of non-events is not used, however, there are consequences to proceeding simply with the usual logit model. Specifically, the following two points are made:

• Parameters for logistic regression are well known to be biased in small samples, but the same bias can exist in large samples if the event is rare.
• Even a bias-corrected estimator for the model parameters does not necessarily lead to optimal predicted probabilities.

If we can determine the amount of bias in a parameter estimate, we can simply subtract it out.

\[\tilde{\boldsymbol{\beta}} = \hat{\boldsymbol{\beta}} – bias(\hat{\boldsymbol{\beta}})\] Based on McCullagh and Nelder’s (1989) foundational work on generalized linear models, the bias for any GLM is:

\[\text{bias}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^{T}\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{W}\xi\] where \(\xi\) is a function of the first and second derivatives of the inverse link function and \(\mathbf{W}\) is, in the logit case, a function of the observed rate of events. The innovation of the King and Zeng estimator given undersampling is to respecify the usual logistic likelihood to be weighted based on the known proportion of events in the population and the proportion of events in the sample. Hence, the estimator adjusts for non-random case selection introduced by the case-control sampling, if used. In the absence of down-sampling non-events, all cases are weighted equally (\(w_i = 1\)), and the bias can be calculated without the weights.

King and Zeng demonstrate how the bias operates on the intercept term in a simple model with one predictor. The bias can be shown (their Appendix D) to be

\[E(\hat{\beta}_0 – \beta_0) \approx \frac{\bar{\pi} – 0.5}{n\bar{\pi}(1 – \bar{\pi})}\]

where \(\bar{\pi}\) is the proportion of events in the data. For rare events, \(\bar{\pi}\) will be less than .5, making the bias negative. The denominator shows that bias is driven by two factors, sample size and the proportion of events. As \(n\) increases, or as \(\bar{\pi}\) increases towards .5 (an even distribution of events and non-events), the bias gets smaller because the denominator gets bigger. Predicted probabilities calculated from the biased estimate will then tend to be too small, as moving the intercept up increases predicted probabilities.

However, even this adjustment may still yield probability estimates that are too low. This becomes clear when one recalls that logistic regression is a linear model of a hypothetical unobserved latent variable, \(y^*\) (see a prior blog post), which takes on a value of zero below a threshold \(\tau\) and one above it. Given a model-based distribution of possible values on the \(y^*\) scale for case i, the predicted probability is the area under the curve above the threshold.

Assuming a logistic distribution, this leads to the usual estimator for the predicted probabilities, \(\frac{1}{1 + e^{-\mathbf{x}_i\boldsymbol{\beta}}}\). However, the \(\boldsymbol{\beta}\) are estimated with some amount of uncertainty, whereas the usual calculation assumes that the \(\boldsymbol{\beta}\) are known perfectly. An alternative estimator proposed by King and Zeng averages over the uncertainty in \(\boldsymbol{\beta}\), which turns out to be equivalent to widening the prediction interval on \(y^*\).

To illustrate, the model makes a prediction for case \(i\) on the scale for \(y^*\). The probability is the area under the distribution that exceeds the threshold.

By averaging over the uncertainty in \(\boldsymbol{\beta}\), the distribution becomes wider, and a larger amount of the area under the curve exceeds the threshold.

In the first figure, the area under the curve above the threshold (i.e. the probability of an event) was 0.023. In the second figure, the variance increased, and the area under the curve above the threshold (or probability of an event) grew to 0.091. In other words, upping the variance for \(y^*\) increased the predicted probability that an event occurs.

Probabilities estimated in the context of rare events tend to be low, and it may even be the case that no probabilities exceed .5. By making an adjustment to the predicted probability estimates, in this case one that is based on the covariance matrix of the estimated parameters, we end up (correctly) with more predicted cases. King and Zeng’s adjusted probability estimator is

\[ p(Y_i = 1) \approx \tilde{\pi}_i + C_i\]

where

\[C_i = (0.5 – \tilde{\pi_i})\tilde{\pi_i}(1 – \tilde{\pi_i})\mathbf{X}V(\tilde{\boldsymbol{\beta}})\mathbf{X}^T \]

The tilde ~ is used to represent the fact that the probabilities come from the bias-adjusted logit model described above.

A Demonstration in R

Say we want to predict an event that occurs 3% of the time in the population. How well do we do with out-of-sample predictions using the regular logit model versus the bias-corrected version?

To answer this question, we’ll generate fake data from a simple model. First, generate a random independent variable \(x_1 \sim \mathcal{N}(0, 1)\). \(Y^*\), the latent continuous value of the observed outcome, will be a linear function of this variable plus some noise drawn from a standard normal distribution, \(Y^* = \beta_0 + \beta_1 x_1 + \mathcal{N}(0, 1)\) (yes, the data generating process is more probit than logit, but variances of 1 are easier to work with than \(\frac{\pi^2}{3}\)). Finally, \(y_i\) is a discrete representation of \(Y_i^*\) such that, if \(Y_i^* >= 0\) then \(y_i = 1\). Otherwise \(Y_i = 0\).

The following code generates the data.

set.seed(12345)

signal_to_noise <- .2
n               <- 750
prop_ones       <- .03

beta_1 <- sqrt(signal_to_noise)
  
dat <- data_frame(x1 = rnorm(n),
                  intercept = qnorm(prop_ones, sd = sqrt(beta_1 ^ 2 + 1)),
                  y_signal = beta_1 * x1,
                  y_star = intercept + y_signal + rnorm(n),
                  y = y_star >= 0,
                  true_prob = pnorm(y_signal + intercept))

signal_to_noise, n, and prop_ones are inputs we could theoretically vary.

• n is population size, here set to 750. We’ll assess the out-of-sample accuracy by randomly setting aside 1/3 of the cases as a test set.
• signal_to_noise is the ratio of variation explainable by the predictor compared to the noise term, here set to .2.
  
• prop_ones is the proportion of events in the data, here set as 3%.

dat %>%
  count(y) %>%
  mutate(prop = round(n/sum(n),2))
## # A tibble: 2 x 3
##   y         n  prop
##   <lgl> <int> <dbl>
## 1 FALSE   726  0.97
## 2 TRUE     24  0.03

Next, split the data into 66% for model fitting and 34% for out of sample prediction, then fit the traditional logistic GLM.

train <- dat %>%
  sample_frac(.66)

test <- dat %>%
  anti_join(train, by = c("x1", "y_star", "y", "true_prob"))

X_train <- train %>%
  select(x1)

X_test <- test %>%
  select(x1)

fit <- glm(as.factor(y) ~ x1, data = train, family = "binomial")

Extract coefficients and the coefficient covariance matrix.

coefs <- fit %>% 
    coef() %>% 
    as.numeric()
  
V <- vcov(fit)

Although the predict.glm method will return predicted probabilities, we’ll create our own function that can be fed the bias-corrected estimates.

logisticPred <- function(X, coef) {
  
  X %>%
    na.omit() %>%
    mutate(int = 1) %>%
    select(int, everything()) %>%
    as.matrix(.) %*% coef %>%
    as.vector() %>%
    (function(x) 1 / (1 + exp(-x)))
  
}

Using the training and test data:

pred <- X_train %>%
  logisticPred(coefs)

test_pred <- X_test %>%
  logisticPred(coefs)

The following figure shows the distribution of probabilities calculated for the test set based on the glm fit to the training data. The xs along the horizontal axis represent cases where an event actually occurred.

None of the probabilities come close to the .5 threshold, so no events are predicted in the test data.

What happens if we make the recommended adjustments to the parameter estimates? The formulas from the King and Zeng paper require linear algebra, so convert the data frame to a matrix. bias is a column in the design matrix representing the intercept.

X_matrix <- X_train %>%
  mutate(bias = 1) %>%
  select(bias, everything()) %>%
  as.matrix()

X_test_matrix <- X_test %>%
  mutate(bias = 1) %>%
  select(bias, everything()) %>%
  as.matrix()

Use formulas to calculate the amount of bias in the coefficients.

W <- diag(pred * (1 - pred))
  
Q <- X_matrix %*% solve(t(X_matrix) %*% W %*% X_matrix) %*% t(X_matrix)

e <- 0.5 * diag(Q) * (2 * pred - 1)

bias <- (solve(t(X_matrix) %*% W %*% X_matrix) %*% t(X_matrix) %*% W %*% e) %>%
  as.numeric()

Now subtract the bias and get the corrected variance.

unbiased_coefs <- coefs - bias

updated_var <- (nrow(X_train) / (nrow(X_train) + ncol(X_train) + 1)) ^ 2 * V

To check our algebra, compare the results to what the Zelig package produces:

re_model <- zelig(y ~ x1, data = train, model = "relogit", cite = FALSE) 

unbiased_coefs
## [1] -3.909300  1.170486
coef(re_model)
## (Intercept)          x1 
##   -3.909300    1.170486
sqrt(diag(updated_var))
## (Intercept)          x1 
##   0.3788302   0.2743757
sqrt(diag(as.matrix(vcov(re_model)[[1]])))
## (Intercept)          x1 
##   0.3803608   0.2754843

The parameter estimates match, but the standard errors are different. Apparently, zelig reports bias corrected parameters but reports SEs from the regular GLM (!), as seen in the output from the next two code chunks.

zelig(y ~ x1, data = train, model = "relogit", cite = FALSE) %>% 
  summary
## Model: 
## 
## Call:
## z5$zelig(formula = y ~ x1, data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1729  -0.2897  -0.1987  -0.1343   3.0611  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -3.9093     0.3804 -10.278  < 2e-16
## x1            1.1705     0.2755   4.249 2.15e-05
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 154.65  on 494  degrees of freedom
## Residual deviance: 132.70  on 493  degrees of freedom
## AIC: 136.7
## 
## Number of Fisher Scoring iterations: 7
## 
## Next step: Use 'setx' method
glm(y ~ x1, data = train, family = "binomial") %>%
  summary
## 
## Call:
## glm(formula = y ~ x1, family = "binomial", data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1770  -0.2815  -0.1915  -0.1283   3.0905  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.9838     0.3804 -10.474  < 2e-16 ***
## x1            1.1958     0.2755   4.341 1.42e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 154.65  on 494  degrees of freedom
## Residual deviance: 132.70  on 493  degrees of freedom
## AIC: 136.7
## 
## Number of Fisher Scoring iterations: 7

Now we’ll make predictions on the test set without the additional adjustment for predicted probabilities but with the bias-corrected coefficients.

unbiased_coef_pred <- X_test %>%
  logisticPred(unbiased_coefs)

How many predicted events do we have?

We haven’t done much better. What if we now add the correction to the estimate of the predicted probability? To facilitate calculations, let \(\zeta = (.5 - \tilde{\pi}_i)\tilde{\pi}_i(1 - \tilde{\pi}_i)\) and \(\eta = \mathbf{x}V(\tilde{\boldsymbol{\beta}})\mathbf{x}^T\).

new_X_test <- test %>%
  as_data_frame %>%
  mutate(p = unbiased_coef_pred) %>%
  mutate(zeta = (.5-p)*p*(1-p),
         eta  = (1*updated_var[1,1] + x1*updated_var[1,2])*1 + 
           (1*updated_var[2,1] + x1*updated_var[2,2])*x1) %>%
  mutate(updated_p = p + zeta*eta) %>%
  mutate(event = updated_p > .5)

corrected_pred <- unlist(new_X_test$updated_p)

And the distribution of probabilities.

In this fake data set, we still failed to predict a single case, though we’ve slightly pushed up the probabilities for events.

King and Zeng provide results from a variety of Monte Carlo simulations that vary the proportion of events, the sample size, and the value of \(\beta_0\) and \(\beta_1\). The results find that very rare events (e.g. 0.15% of cases) are hard to predict even with the corrected method, whereas there is not much difference when the signal-to-noise ratio is very high. Somewhat less rare events predicted from noisier models seem to benefit the most from the adjustment.

It is well known that the logistic model does an inferior job at predicting out-of-sample events relative to machine learning methods like random forests and GBMs, and we will continue to use these when the our projects emphasize accurate forecasts. On the other hand, combining the bias-corrected estimator with propensity-score matched or weighted samples is an option worth considering in statistical projects.