# AIC model selection is keeping a variable with p = 0.47

I am modeling migration departure timing for swallows to try and figure out which of the predictor variables that I have data for influenced departure timing. All of the predictors are variables that I have reason to expect might affect departure timing: sex, wing chord (proxy for body size), mass (proxy for body condition), breeding latitude, and year (I have two years of data). I also have a random effect for tagging location, so the full model looks like this:

lmer(departureDate ~ sex + wingChord + mass + latitude + year + (1|tagSite))

I checked for correlation between my numeric predictors and the highest correlation was $$r = 0.3$$ between latitude and mass.

mass wingChord latitude
mass 1.00000000 0.10356130 0.30809965
wingChord 0.10356130 1.00000000 0.08280734
latitude 0.30809965 0.08280734 1.00000000

The model summary for the full model shown above is this (with REML = TRUE, I use REML = FALSE for AICc model selection as described below; continuous predictors have been scaled):

AIC BIC logLik deviance df.resid
4042.6 4077.2 -2013.3 4026.6 549

Scaled residuals:

Min 1Q Median 3Q Max
-2.5559 -0.5642 -0.0359 0.4780 6.0517

Random effects: | Groups | Name| Variance | Std.Dev | | --- | --- | --- | --- |--- |--- | | tagSite | (Intercept) | 14.32 | 3.784 | | Residual | 76.47 | 8.745 |

Number of obs: 557, groups: tagSite, 19

Fixed effects:

Estimate Std.Error df t value Pr(>abs(t))
(Intercept) 199.9260 1.2695 35.6941 157.485 < 2e-16 ***
latitude 3.6840 1.0170 16.6520 3.622 0.002165 **
wingChord -1.2529 0.4201 555.6471 -2.982 0.002985 **
mass 0.2959 0.4180 552.4835 0.708 0.479298
sexM -1.3458 0.7851 549.1789 -1.714 0.087046
year2023 -3.7353 1.0155 475.3548 -3.678 0.000262 ***

Correlation of Fixed Effects:

(Intr) latitude wingChord mass sexM
latitude 0.157
wingChord -0.016 -0.074
mass 0.047 -0.099 -0.070
sexM -0.263 -0.036 0.016 0.170
year2023 -0.547 -0.111 0.064 -0.152 0.077

I am currently trying to use AICc model selection to come up with the most parsimonious model for departure timing. I am not using it to predict. I am just trying to understand what are the most important factors (from the ones that I have data for) that drive migratory departure timing. I have checked the residuals for the full model above and they look good, so the fit of the model itself does not seem to be an issue from what I can tell.

I've tried running the full set of possible models for all 32 combinations of my fixed effects (including a null model with just the random effect) and selecting the model with the lowest AIC. Ive also tried calculating the ratio of the estimate/SE (from the full model with REML = TRUE) and sequentially eliminating predictors in increasing order of this ratio (so, starting with the predictor with the lowest number) and stopping when the AIC starts to increase.

In both cases, I am left with the full model as having the lowest AICc. What I am confused about is why mass is left in the model when the $$p$$-value is so high. I understand (based on other posts on this issue and based on Sutherland et al. 2023 (https://doi.org/10.1098/rspb.2023.1261) that AICc model selection will retain variables with a $$p$$-value $$< 0.157$$. However, the $$p$$-value for mass is $$0.48$$, so why is the AICc for the model with mass almost $$6$$ points lower than the model without mass? (AICc for the full model is $$4042.861$$, AICc for the model with only mass removed is $$4048.136$$).

I also understand that there are differing opinions on model selection and AIC. I am unsure if this is the final route I will go for this analysis, however I am concerned that the fact that the model selection process using AICc is keeping mass, despite it's high $$p$$-value, is indicative of a bigger issue with my data/analysis.

I am not able to post my raw data but I would be able to post other statistical test results or figures if required to answer this question.

• @purplebubbles93 do any of the predictors have missing values? AIC only works when the exact same data set us used in all models (and I know I've lost a day or two of my life dealing with this issue before ....) overall this does seem wonky Commented Apr 4 at 0:01
• @NBrouwer I did have one NA in there for mass, which I hadn't realized, but I still have the same issue now that I've removed that individual. Commented Apr 4 at 15:26
• Weird. Since I have other good suggestion I'll just toss things out there that may be useful: 1) how are you calculating AICc? Try two different packages and see what happens (I like bbmle and it's AICtab). 2) check the raw likelihood and df for each model and make sure they check out (this helps esp when there are NAs lurking) 3) look at BIC or other ICs to see if they are also influenced by mass Commented Apr 4 at 22:21
• 4) calculate p values more than 1 one way 5) use a lasso on the full model 6) try w/ and w/o rand effect 7) run model on different portions of the data Commented Apr 4 at 22:31
• AIC (AICc) doesn't work well for small $n$. AIC is only asymptotically correct for large $n$ if the number of variables is kept small, and $n>400$ has been reported for this. If the number of variables approaches the number of sample realizations all bets are off. One problem with AIC is that there is no probability of a correct model selection available, so something else is needed, e.g., cross-validation, simulation studies, or whatever.
– Carl
Commented Apr 4 at 22:59

I really don't like the approach you have proposed here for model selection. If you have a set of a-priori defined models, test those directly. Then calculate AIC and see which fits look better (regardless of their $$p$$ values). From there you have tested a theoretical model with a mathematical model, and the AIC can help validate that selection in some way.

Fitting every combination of the predictors and then trying to find the lowest AIC is odd. You are potentially capitalizing on complete chance effects in doing so, and this masks the study as confirmatory when it is slowly becoming exploratory instead. More importantly, if your goal is theory testing, then having weak predictors isn't a problem. In fact, if you find that some important factor actually has no effect, people need to learn from that and this would be interesting to report. But AIC will typically favor models with strong effects, so naturally omitting what should be important tests of your theory may be completely lost using the method you are seeking here.

Finally, the $$p$$ values should be much lower on the list of things to consider here. These are just probabilities and say nothing about the practical or theoretical implications of your model.

• +1. Re "p-values ... are just probabilities:" they aren't even that. They are hypothetical probabilities under the assumed null hypotheses.
– whuber
Commented Apr 3 at 21:28
• Yes, that is an important qualifier. Commented Apr 3 at 21:30
• +1, because your answer is useful. However, it does not really answer the OP's question, does it? We have a predictor whose coefficient has $p=0.47$, so its effect should be anything but "strong" (whatever our thoughts about p values and NHST), so if AIC(c) prefers models with strong effects, why is it lower when we include this weak predictor? Commented Apr 4 at 7:08
• @ShawnHemelstrand thank you for the response, but I'm not sure I understand it. I don't think that I have a set of a-priori defined models, I have a set of predictors that I think may affect departure timing (based on other studies with small numbers of the same species as mine or on different bird species), and I am wondering which (if any, it's fine if the answer is none) may have affected departure timing for the individuals in my study. I actually do see my study as more exploratory than confirmatory, as there is really limited research in this area for my study species. Commented Apr 4 at 15:21
• „Finally, the p values should be much lower on the list of things to consider here“ This is missing the point of the question. Whether or not the p-values are on the list of things to consider, the result tables from the lmerTest function indicate a p-value that is conflicting with the AICc values. Setting p-values aside as being statistically not relevant in this exploratory approach, they still are relevant as a way of quality control. The selected model via AIC should not have parameters with p-values a lot much above 0.15 Commented Apr 25 at 11:49

Shawn's idea of avoiding stepwise variable selection should be heeded. As has been discussed countless times on this site, variable selection is a disaster. The first thing it does is to ruin standard errors of remaining variables. Also it has almost no chance of finding the right variables. I would try to pose one model and compute confidence intervals for the importance of each variable in that model if you want to learn about what's important. The confidence intervals will often tell you that the data don't have sufficient information to tell you what's important about the data. See this.

• But why is the difference in AIC such high while the p-value is so much insignificant? That's the question. Commented Apr 25 at 11:51

The question is why it is possible that a feature with a low significant parameter is retained in a stepwice AIC selection procedure. t-values typically relate to f-statistics, likelihood ratios, and eventually AIC scores. Question about this equivalence are:

What is happening here?

It is difficult to say whether the below situation is the case here without much further information. It would be helpful if a reproducible example can be given.

If we compute the likelihood ratio based on the t-test and degrees of freedom,

$$\log(LR) = \frac{df}{2} \log\left(1+\frac{t^2}{df}\right)$$

and if we compute the likelihood based on AIC

$$\log(LR) = \frac{AIC_0-AIC_1}{2} - 1$$

then we should get reasonably close results. An example computation below shows that likelihood ratio, F-test and t-test should match closely.

If I compute this for your case then the $$\log(LR)$$ is approximately 0.25, and it is weird that the full model ends up with an AICc score 6 points higher.

There are maybe some discrepancies in how the different functions compute the degrees of freedom and deviance, but it seems more like a computation error. So youay need to dig into the code and for the question here provide a reproducible example.

library(lmerTest)
n = 100
m = 10
set.seed(1)

### random effect
z = sample(1:m, n, replace = T)
r_eff = rnorm(m)[z]

### random error
error = rnorm(n)

### fixed effects
x1 = rnorm(n,0,1)
x2 = rnorm(n,0,1)
x3 = rnorm(n,0,1)
x4 = rnorm(n,0,1)
mass = x1+x2+x3+x4 + rnorm(n,0,1)

z = as.factor(z)
y = r_eff + x1 + x2 + x3 + x4 + 0.05 * mass + error

### modeling
mod1 = lmer(y ~ x1 + x2 + x3 + x4 + mass + (1|z), REML = 0)
mod0 = lmer(y ~ x1 + x2 + x3 + x4 + (1|z), REML = 0)

summary(mod1)
anova(mod1)

### log likelihood ratio 0.09421622
- (AIC(mod1)-(AIC(mod0)+2))/2

### log likelihood computed from t-statistic and degrees of freedom
### the value is 0.09422343
t = summary(mod1)$$coefficients[6,4] f = t^2 df = summary(mod1)$$coefficients[6,3]
df/2*log(1+f*(1)/(df))