使用AIC逐步回归拟合含交互项三阶多项式模型的共线性问题咨询

Hey there, I totally get how frustrating that wonky β coefficient issue from multicollinearity can be when building polynomial regression models—those correlated raw polynomial terms are definitely the culprit here. Let’s walk through two reliable fixes tailored to your use case:

1. Use Orthogonal Polynomials (Quickest Model Selection Fix) #

The core problem is that raw terms like a, I(a^2), and I(a^3) are highly correlated. R’s built-in poly() function creates orthogonal polynomial terms by default, which eliminates this collinearity entirely without altering the model’s predictive power or AIC values. Perfect for stepwise selection!

Here’s your adjusted code:

library(MASS)  # Required for stepAIC

# Build initial model with orthogonal 3rd-order polynomials
lm_poly_ortho <- lm(y ~ poly(a, 3) + poly(b, 3) + c + a:b, data = my_data)

# Run stepwise AIC regression to find the optimal model
optimal_model <- stepAIC(lm_poly_ortho, direction = "both")

A quick heads-up: Orthogonal polynomial coefficients aren’t directly interpretable as "a’s linear effect" or "b’s cubic effect," but that doesn’t impact model selection—AIC depends only on fit quality, so the stepwise process will still pick the best term combination. If you later need raw polynomial coefficients for interpretation, you can reverse-engineer the orthogonal transformation or use predict() with type = "terms".

2. Center Predictors (For Interpretable Raw Coefficients) #

If you need to keep raw polynomial terms to interpret coefficients in your original variable context, centering predictors first will drastically reduce collinearity. Centering subtracts the mean from each variable, breaking the strong correlation between low-order and high-order terms (e.g., a and a^2).

Try this workflow:

library(MASS)
library(car)  # For VIF checks to verify collinearity is fixed

# Center your predictors (handle NA values with na.rm = TRUE if needed)
my_data <- transform(my_data,
                     a_cent = a - mean(a, na.rm = TRUE),
                     b_cent = b - mean(b, na.rm = TRUE))

# Build model with centered polynomial terms
lm_poly_cent <- lm(y ~ a_cent + I(a_cent^2) + I(a_cent^3) + b_cent + I(b_cent^2) + I(b_cent^3) + c + a_cent:b_cent, data = my_data)

# Verify collinearity is reduced (VIF < 5 is ideal)
vif(lm_poly_cent)

# Run stepwise AIC regression
optimal_model_cent <- stepAIC(lm_poly_cent, direction = "both")

This approach keeps your β coefficients stable and interpretable, while letting stepAIC do its job properly.

Bonus: Diagnose Collinearity First #

Before and after making changes, use the vif() function from the car package to check for multicollinearity. A VIF (Variance Inflation Factor) above 5 (or 10, depending on convention) signals significant collinearity that’s distorting your coefficients.

内容的提问来源于stack exchange，提问作者Jans.K.