Does automatic model selection via AIC bias the p-values of the selected model? [Looking for simulation-based evidence]

Question

Let's say I run a procedure where I fit every possible model given some set of covariates and I select the model with the minimum AIC. I know that if my selection criteria was based on minimizing p-values, the p-values of the selected model would be misleading. But what if my selection criteria was AIC alone? To what extent would this bias the p-values?

I had assumed the effect on p-values would be real, but negligible, but came across this paper, which proves the following:

P values are intimately linked to confidence intervals and to differences in Akaike's information criterion (ΔAIC), two metrics that have been advocated as replacements for the P value.

If this is true, does it imply that p-values are misleading after automatic selection based on AIC? To what extent will they be biased, and what determines this?

Answer 1

Yes, it does bias the p-values

Your intuition on this is correct --- generally speaking, whenever we select a model via optimisation we bias the resulting p-value for any tests that fail to account for that optimisation step. This is true regardless of whether we optimise over p-values, or the maximum log-likelihood, or AIC, BIC, etc.

By minimising AIC you are maximising the log-likelihood (and therefore the likelihood) over each model with the same number of terms. To see this, suppose we let $\mathcal{M}$ denote an individual model and we let $\mathscr{M}$ denote the class of all models under consideration. If we let $\mathscr{M}_k$ be the class of models with $k$ parameters then we can write the maximum log-likelihood over this class as:

$$R(k) \equiv \max_{\mathcal{M} \in \mathscr{M}_k} \hat{\ell}_\mathcal{M}$$

The minimum AIC over the class of all models is related to this quantity by:

$$\begin{align} \min_{\mathcal{M} \in \mathscr{M}} \text{AIC}_\mathcal{M} &= \min_{\mathcal{M} \in \mathscr{M}} (2 k_\mathcal{M} - 2 \hat{\ell}_\mathcal{M}) \\[6pt] &= \min_{k} \min_{\mathcal{M} \in \mathscr{M}_k} (2 k - 2 \hat{\ell}_\mathcal{M}) \\[6pt] &= \min_{k} \Big( 2 k - 2 \max_{\mathcal{M} \in \mathscr{M}_k} \hat{\ell}_\mathcal{M} \Big) \\[6pt] &= \min_{k} \Big( 2 k - 2 R(k) \Big). \\[6pt] \end{align}$$

As you can see, minimising AIC is equivalent to a two step process: first we find all the models that maximise the log-likelihood compared to other models with the same number of parameters; then we minimise the quantity $(2k-R(k))$ using the models from the first step (which each have different values for $k$). This means that when you select a model by minimising AIC, you are implicitly maximising log-likelihood over subclasses of models. By optimising the log-likelihood, you are inducing a bias towards parameter estimates that give a highly "peaked" likelihood function (giving a higher maximum), which is going to tend to happen when you use parameter values that are far from the "null" hypotheses in standard model tests. You then optimise over $k$ which exaccerbates this further by selecting the number of model terms that gives a highly "peaked" likelihood function relative to the number of model terms. This means that your optimisation procedure is going to tend to give you a model with "significant" model terms, biasing your p-values downward in these tests.

In regard to the size of the bias, this is going to depend largely on the size of the model class over which you are optimising. If you optimise over a small class of models then the bias may be modest, but if you optimise over a large class of models then the resulting bias will be large. Since you have stated that you are using an "all possible models" approach, that means that you are optimising over $2^m$ possible models where you have $m$ model terms to consider. The size of the model space grows exponentially in $m$ so you rapidly get to a large model space which will lead to large bias when optimising over this space.

Correcting for this kind of bias is complicated, and it generally entails running simulations where you perform the same (AIC) optimisation over a set of null models where all model obey the null hypothesis of interest. If you can form suitable pivotal quantities, this kind of simulation can give you an estimated null distribution for your test statistics under the optimisation method and you can then use this to estimate the true (unbiased) p-value of the test. This is quite a complicated exercise, as it generallly requires you to custom-program a computation that loops over model fitting and optimisation, and computation and extraction of resulting test statistics.

Answer 2

Since simulations were requested as a followup to Ben's excellent answer, here's a simple example.

Model

Consider a polynomial regression model of the form:

$$y_i = \beta_0 + \sum_{d=1}^D \beta_d x_i^d + \epsilon_i \quad \quad \epsilon_i \underset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)$$

The intercept $\beta_0$, coefficients $\beta_1, \dots, \beta_D$, and noise variance $\sigma^2$ will be fit by maximum likelihood for each choice of degree $D$. The degree will be chosen to minimize the AIC.

Hypothesis testing

Suppose we want to test the null hypothesis that all coefficients are zero (except the intercept). Ordinarily, we could use an F test. But, as Ben described, the resulting p values will be biased here because the test doesn't account for the model selection step. We can confirm this by examining the distribution of p values in an example where the null hypothesis is known to be true.

Simulation

Repeat $10^4$ times:

1. Generate a dataset containing $n=20$ points, where explanatory variables $x_i$ and responses $y_i$ are sampled independently from the standard normal distribution. Since they're independent, the null hypothesis is true.

2. Fit polynomial regression models with degree $D$ ranging from 1 to $D_{max}=10$, using maximum likelihood.

3. Select the model with minimum AIC.

4. Run an F test for the selected model, as above. Record the resulting p value.

The parameters here are unrealistic (e.g. nobody would fit a 10th degree polynomial to 20 data points). But, a proper hypothesis test should be able to handle it. I've set things up this way to make the failure/bias more obvious.

Biased p values

Under the null hypothesis, proper p values should be uniformly distributed between 0 and 1. We can check whether this is true, since p values from the simulations are samples from the distribution under the null hypothesis. The histogram below shows that the distribution is decidedly non-uniform, with low p values much more probable than they should be.

This indicates that the p values are optimistically biased, and we'd run into trouble using them for hypothesis testing. A test is deemed significant if the p value falls below a threshold $\alpha$, which specifies the maximum acceptable type I error rate---the probability of wrongly obtaining a significant result, assuming the null hypothesis is true. Since the null hypothesis is true in the simulations, all significant results are type I errors. So, for any nominal type I error rate $\alpha$, we can estimate the actual type I error rate as the fraction of p values that fall below $\alpha$. Given proper p values, the nominal and actual rates should be equal. However, as shown in the plot above, the actual type I error rate exceeds nominal rate for all choices of $\alpha$.

Dependence on model size

As Ben also mentioned, comparing a greater number of models during the model selection step will increase the amount of bias. To show this effect, I repeated the simulation above for different model sizes. In particular, I varied the maximum permitted degree $D_{max}$ of the polynomial regression model. The plot below shows the actual vs. nominal type I error rates for each choice of $D_{max}$.

Notice that the curve is the identity line for $D_{max}=1$, indicating that p values behave as theyy should. This is because no model selection happens in this case--a degree 1 polynomial is the only choice. However, p values are increasingly biased for greater choices of $D_{max}$ because an increasing number of models are compared in each case.

Code

Matlab code implementing the simulations described above:

% parameters
n = 20; % how many data points
Dmax = 10; % maxmimum polynomial degree
nreps = 1e4; % how many simulations

% compute p value for each simulation
p = zeros(1, nreps);
for rep = 1 : nreps
    fprintf('%d/%d\n', rep, nreps);
    
    % generate data
    x = randn(n, 1);
    y = randn(n, 1);
    
    % design matrix for all polynomial degrees
    A = x .^ (1 : Dmax);
    
    % fit polynomial regression models
    % choose degree that minimizes aic
    best_mdl = [];
    best_aic = inf;
    for D = 1 : Dmax
        mdl = fitlm(A(:, 1:D), y);
        aic = mdl.ModelCriterion.AIC;
        if aic < best_aic
            best_mdl = mdl;
            best_aic = aic;
        end
    end
    
    % run F test on model with best aic
    p(rep) = best_mdl.coefTest();
end

% nominal type I error rates
alpha = linspace(0, 1, 100);

% actual type I error rate for each alpha level
fpr = sum(p < alpha(:), 2) / nreps;