Why does AIC select the wrong logistic regression model?

I've noticed a phenomenon with logistic regression: when the probability of success is small and number of trials large the AIC consistently selects the wrong models. In fact the more wrong the model the lower the AIC.

Here's some example code to show what I mean. About $$96\%$$ of the time the model with the random probe gets chosen. $$100\%$$ of the time the model with $$3$$ random probes gets chosen. The models with probes have AIC between $$300$$ and $$2000$$ units lower.

If I replace the $$\lambda$$ in the Poisson distribution with something smaller, like $$1000$$, this doesn't happen anymore - the model with the random probe gets selected about $$9\%$$ of the time, the model with 3 probes $$6\%$$ of the time.

This seems very odd to me. As the number of trials goes up the number of successes does too, so why does the metric break down?

nreps <- 1000
registerDoMC(16)

N <- 40000
lambda = 1e7

result <- foreach(rep = 1:nreps) %dopar% {

# Number of trials
n_trials <- rpois(N, lambda = lambda)

# Probability of success is very small and depends on x
x <- rnorm(N)

rate <- rep(1e4, N)
rate[x > 0] <- 2e4

p <- rexp(N, rate)

# Convert probability of success into number of successes in n_trials
success <- floor(n_trials * p)
failure <- n_trials - success

# Set up target
y <- data.frame(success = success, failure = failure)
y <- as.matrix(y)

# Random numbers
probe0 <- rnorm(N)
probe1 <- rnorm(N)
probe2 <- rnorm(N)

# Fit models
fit  <- glm(y ~ x, family = 'binomial')
fit0 <- glm(y ~ x + probe0, family = 'binomial')
fit1 <- glm(y ~ x + probe0 + probe1, family = 'binomial')
fit2 <- glm(y ~ x + probe0 + probe1 + probe2, family = 'binomial')

res <- data.frame(fit = AIC(fit), fit0 = AIC(fit0), fit1 = AIC(fit1), fit2 = AIC(fit2))

return(res)

}

result <- Reduce(rbind, result)

mean(result$$fit > result$$fit0)
mean(result$$fit > result$$fit2)
$$$$

• Without looking into the details, could it be that AIC selects the best forecasting model (maximimizing the expected likelihood of a new, unseen observation) instead of the true model? That would be exactly as expected. – Richard Hardy Mar 3 at 21:27
• That doesn't explain why AIC gets uniformly better the more random variables are involved or why a smaller number of trials makes the problem go away. – badmax Mar 3 at 23:08

Your outcome variable $$y$$ is very overdispersed relative to a binomial, so the binomial model is misspecified. From one example of fit2, modified to quasibinomial:

> summary(glm(formula = y ~ x + probe0 + probe1 + probe2, family = "quasibinomial"))

Call:
glm(formula = y ~ x + probe0 + probe1 + probe2, family = "quasibinomial")

Deviance Residuals:
Min       1Q   Median       3Q      Max
-61.135  -22.809   -9.250    9.047  189.121

Coefficients:
Estimate Std. Error   t value Pr(>|t|)
(Intercept) -9.539830   0.005434 -1755.604   <2e-16 ***
x           -0.263607   0.005261   -50.109   <2e-16 ***
probe0       0.007280   0.005243     1.389    0.165
probe1      -0.001967   0.005278    -0.373    0.709
probe2      -0.004626   0.005269    -0.878    0.380
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 823.1188)

Null deviance: 28757637  on 39999  degrees of freedom
Residual deviance: 26686790  on 39995  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5
`

If you don't have any allowance for the overdispersion, larger models will always fit better. Basically, you are accepting a new model if it improves the binomial deviance by 2, and a random predictor will improve the binomial deviance by an average of over 800.

You can't just use AIC with quasilikelihood, you'd need something like Takeuchi's criterion, or you'd need to fit something like a beta-binomial model.