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)
```

 A: 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.
