R中多项式逻辑回归特征选择：如何保留高阶项对应的低阶协变量？

方法一：自定义约束的向后逐步选择 #

手动实现带层次检查的向后选择是最直接的方式，核心逻辑是禁止移除仍有高阶项在模型中的低阶变量，具体步骤：

1. 先构建包含所有多项式项的全模型（建议手动生成原始多项式项，比如x1_sq = x1^2，方便后续识别阶次；若用poly()生成正交多项式，需保证变量命名可区分阶次）
2. 循环筛选可移除变量：
   • 对每个候选变量，先检查同特征的高阶项是否在模型中：比如要移除主效应x1，若模型里存在x1_sq/x1_cu这类项，则跳过该变量；
   • 优先移除无高阶依赖、且移除后AIC下降最多（或性能损失最小）的变量；
3. 重复上述步骤，直到无法移除符合约束的变量，或达到预设简化目标。

示例代码片段：

library(MASS)
# 生成模拟数据
set.seed(123)
x1 <- rnorm(100)
x2 <- rnorm(100)
y <- factor(ifelse(0.5 + 2*x1 + 1.5*x1^2 + 1*x2 > rnorm(100), "1", "0"))
df <- data.frame(y, x1, x1_sq = x1^2, x2, x2_sq = x2^2)

# 全模型
full_model <- glm(y ~ x1 + x1_sq + x2 + x2_sq, family = binomial)

current_model <- full_model
while(TRUE) {
  # 提取当前模型的变量（截距除外）
  vars <- attr(terms(current_model), "term.labels")
  aic_values <- c()
  valid_vars <- c()
  
  for(var in vars) {
    # 提取基变量名称（移除阶次后缀）
    base_var <- gsub("_sq|_cu", "", var)
    # 查找同特征的高阶项
    higher_terms <- grep(paste0("^", base_var, "_"), vars, value = TRUE)
    
    # 若当前变量是低阶主效应且存在高阶项，跳过移除
    if(length(higher_terms) > 0 && var == base_var) next
    
    # 尝试移除变量并计算AIC
    temp_model <- update(current_model, paste0(". ~ . - ", var))
    aic_values <- c(aic_values, AIC(temp_model))
    valid_vars <- c(valid_vars, var)
  }
  
  if(length(valid_vars) == 0) break # 无符合条件的可移除变量
  
  # 找到移除后AIC最小的变量
  best_idx <- which.min(aic_values)
  best_var <- valid_vars[best_idx]
  new_aic <- aic_values[best_idx]
  
  # 若移除后AIC上升，停止选择
  if(new_aic >= AIC(current_model)) break
  
  current_model <- update(current_model, paste0(". ~ . - ", best_var))
  cat("移除变量:", best_var, "当前AIC:", AIC(current_model), "\n")
}

# 查看最终模型
summary(current_model)

方法二：利用glmulti包实现层次约束选择 #

glmulti包专门用于多模型推断与特征选择，支持直接开启层次约束，自动保证选中高阶项时必须保留对应低阶项：

1. 安装并加载包：install.packages("glmulti"); library(glmulti)
2. 定义模型公式，包含所有多项式项；
3. 运行特征选择时设置hierarchy = TRUE，指定评估准则（AIC或交叉验证）。

示例代码：

library(glmulti)
# 定义模型框架（用I()包裹多项式项）
model_formula <- y ~ x1 + I(x1^2) + x2 + I(x2^2)

# 带层次约束的向后选择，逻辑回归，AIC准则
glmulti_fit <- glmulti(model_formula, data = df, family = binomial, 
                       method = "h", 
                       hierarchy = TRUE, 
                       crit = "aic")

# 查看最优模型
summary(glmulti_fit)
best_model <- getModel(glmulti_fit)

方法三：遗传算法中加入层次约束惩罚 #

若使用遗传算法（如GA包），可通过自定义适应度函数实现约束：

1. 将特征选择编码为二进制向量（1表示保留变量，0表示移除）；
2. 计算适应度前先检查编码是否违反层次规则：若某高阶项被保留但对应低阶项被移除，给该个体极低的适应度（惩罚）；
3. 对符合约束的个体，用AIC或交叉验证性能（如ROC AUC）计算适应度。

示例代码片段：

library(GA)
# 变量列表
vars <- c("x1", "x1_sq", "x2", "x2_sq")
n_vars <- length(vars)

# 自定义适应度函数
fitness <- function(x) {
  selected_vars <- vars[x == 1]
  
  # 检查层次约束
  valid <- TRUE
  if("x1_sq" %in% selected_vars && !"x1" %in% selected_vars) valid <- FALSE
  if("x2_sq" %in% selected_vars && !"x2" %in% selected_vars) valid <- FALSE
  
  if(!valid) return(-Inf) # 违反约束，给予惩罚
  
  if(length(selected_vars) == 0) return(-Inf) # 至少保留一个变量
  
  # 构建模型并返回负AIC（GA默认最大化适应度）
  formula <- as.formula(paste0("y ~ ", paste(selected_vars, collapse = " + ")))
  model <- glm(formula, data = df, family = binomial)
  return(-AIC(model))
}

# 运行遗传算法
ga_fit <- ga(type = "binary", fitness = fitness, nBits = n_vars, 
             popSize = 50, maxiter = 100, seed = 123)

# 查看最优特征选择
best_solution <- vars[ga_fit@solution[1,] == 1]
best_solution

注意事项 #

• 变量命名需统一规则：比如用x1、x1_sq、x1_cu区分主效应、二次项、三次项，方便代码识别同特征的不同阶次；
• 若使用正交多项式（poly(x1, 2)），生成的变量为poly(x1,2)1、poly(x1,2)2，可通过正则表达式提取基变量来检查层次；
• 若需样本外评估，可将适应度函数替换为交叉验证的预测性能（如ROC AUC），替代AIC。

内容的提问来源于stack exchange，提问作者Aegis