AIC model selection is keeping a variable with p = 0.47

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

The other answers make good comments about stepwise regression, but they are not answers to the posed question.

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:

• Replicate t or F test from regression using regression likelihoods
• Relationship between the t-statistic and F-statistic in simple linear regression

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


### with your values
f = 0.708^2
df = 552.4835
df/2*log(1+f*(1)/(df))
### the result is 0.2505184

The example above has a model because I was trying to reproduce the error by getting some weird situations like parameters that correlate with the random effect, and small correlatio s but high collinearity. But I couldn't get a situation where the t-value can be so small for such large AIC difference.