如何在caret::train的glmStepAIC中维护回归层级并修复自定义模型报错

问题说明 #

使用caret::train调用glmStepAIC执行重复交叉验证时，默认MASS::stepAIC会在保留交互项的情况下移除对应主效应，违反层级回归原则。自定义caret模型试图修复该问题，但运行后出现RMSE、Rsquared、MAE全为NA/NaN的错误，原默认方法可正常运行。

原实现代码 #

基础caret逐步回归代码 #

splitRule <- trainControl(method="repeatedcv", number=10, repeats=5, verboseIter = FALSE)
model <- caret::train(
  OUTCOME ~ PredictorVariable1 + PredictorVariable2 + PredictorVariable3 + 
            PredictorVariable1 * PredictorVariable2,
  data = dataset$trainingSet,
  trControl = splitRule,
  na.action = na.omit,
  method = "glmStepAIC",
  direction = "both"
)
predictions_model <- predict(model, newdata = dataset$testingSet)

报错的自定义模型代码 #

hierarchicalStepAIC <- list(
  library = c("MASS", "stats"),
  type = "Regression",
  parameters = data.frame(parameter = "parameter", class = "character", label = "parameter"),
  grid = function(x, y, len = NULL, search = "grid") data.frame(parameter = "none"),
  
  fit = function(x, y, wts, param, lev, last, classProbs, ...) {
    # Ensure x is a data frame and combine with y
    if (!is.data.frame(x)) x <- as.data.frame(x)
    dat <- x
    dat$.outcome <- y
    
    # Use the formula passed via caret::train (accessed via dots)
    dots <- list(...)
    if (!"formula" %in% names(dots)) {
      stop("Formula not provided to custom model")
    }
    initial_formula <- dots$formula
    
    # Fit the initial GLM model
    initial_model <- glm(initial_formula, data = dat, family = "gaussian", na.action = na.omit)
    
    # Run stepAIC
    step_model <- tryCatch(
      stepAIC(initial_model, direction = "both", trace = FALSE),
      error = function(e) {
        warning("stepAIC failed: ", e$message)
        return(initial_model)  # Fallback to initial model
      }
    )
    
    # Extract terms from the final model
    final_terms <- terms(step_model$formula)
    final_labels <- attr(final_terms, "term.labels")
    
    # Identify interaction terms and required main effects
    interaction_terms <- grep(":", final_labels, value = TRUE)
    if (length(interaction_terms) > 0) {
      required_main_effects <- unique(unlist(strsplit(interaction_terms, ":")))
      missing_main_effects <- setdiff(required_main_effects, final_labels)
      
      # Add back missing main effects if any
      if (length(missing_main_effects) > 0) {
        updated_formula <- reformulate(c(final_labels, missing_main_effects), response = ".outcome")
        step_model <- glm(updated_formula, data = dat, family = "gaussian", na.action = na.omit)
      }
    }
    
    # Return the model
    step_model
  },
  
  predict = function(modelFit, newdata, submodels = NULL) {
    if (!is.data.frame(newdata)) newdata <- as.data.frame(newdata)
    predict(modelFit, newdata = newdata, type = "response")
  },
  
  prob = function(modelFit, newdata, submodels = NULL) NULL,
  levels = function(x) NULL
)

问题根源与修正方案 #

自定义模型报错的核心原因：

1. formula获取逻辑错误：caret向自定义fit函数传递的...中并不包含完整的原始公式，而是拆分后的x（预测变量）和y（响应变量），导致initial_formula无法正确构建。
2. 响应变量命名不匹配：用.outcome作为响应变量名，但caret在重采样时的数据集结构可能不兼容，导致后续预测或指标计算失败。
3. 未传递权重参数：忽略了caret可能传递的wts参数，部分重采样场景下会导致模型拟合异常。

修正后的自定义模型 #

hierarchicalStepAIC <- list(
  library = c("MASS", "stats"),
  type = "Regression",
  parameters = data.frame(parameter = "direction", class = "character", label = "Step Direction"),
  grid = function(x, y, len = NULL, search = "grid") {
    # 支持指定stepAIC的direction参数，可选"both"、"backward"、"forward"
    data.frame(direction = c("both"))
  },
  
  fit = function(x, y, wts, param, lev, last, classProbs, ...) {
    # 确保x是数据框，构建完整数据集
    if (!is.data.frame(x)) x <- as.data.frame(x)
    dat <- cbind(x, .outcome = y)
    
    # 从x的列名构建初始公式（包含所有预测变量，包括交互项）
    pred_vars <- colnames(x)
    initial_formula <- reformulate(pred_vars, response = ".outcome")
    
    # 处理权重参数
    glm_args <- list(formula = initial_formula, data = dat, family = "gaussian", na.action = na.omit)
    if (!is.null(wts)) glm_args$weights <- wts
    
    # 拟合初始GLM
    initial_model <- do.call(glm, glm_args)
    
    # 运行带约束的逐步回归
    step_model <- tryCatch(
      stepAIC(initial_model, direction = param$direction, trace = FALSE),
      error = function(e) {
        warning("stepAIC failed: ", e$message)
        return(initial_model)
      }
    )
    
    # 检查并修复层级约束
    final_terms <- terms(step_model)
    final_labels <- attr(final_terms, "term.labels")
    interaction_terms <- grep(":", final_labels, value = TRUE)
    
    if (length(interaction_terms) > 0) {
      # 提取交互项对应的主效应
      required_mains <- unique(unlist(strsplit(interaction_terms, ":")))
      # 找出被移除的主效应
      missing_mains <- setdiff(required_mains, final_labels)
      
      if (length(missing_mains) > 0) {
        # 更新公式，添加缺失的主效应
        updated_terms <- c(final_labels, missing_mains)
        updated_formula <- reformulate(updated_terms, response = ".outcome")
        
        # 重新拟合模型，保留权重参数
        refit_args <- glm_args
        refit_args$formula <- updated_formula
        step_model <- do.call(glm, refit_args)
      }
    }
    
    return(step_model)
  },
  
  predict = function(modelFit, newdata, submodels = NULL) {
    if (!is.data.frame(newdata)) newdata <- as.data.frame(newdata)
    # 确保新数据包含模型所需的所有变量（包括可能补加的主效应）
    missing_vars <- setdiff(all.vars(modelFit$formula)[-1], colnames(newdata))
    if (length(missing_vars) > 0) {
      warning("New data missing required variables: ", paste(missing_vars, collapse = ", "))
    }
    predict(modelFit, newdata = newdata, type = "response")
  },
  
  prob = function(modelFit, newdata, submodels = NULL) NULL,
  levels = function(x) levels(x$.outcome)
)

关键修改点 #

• 公式构建逻辑：基于x的列名动态生成初始公式，确保兼容caret的重采样数据拆分逻辑。
• 权重参数传递：通过do.call传递wts参数，适配caret的加权重采样场景。
• 层级约束修复：在补加主效应后，使用相同的glm参数重新拟合，保证模型一致性。
• 参数化支持：将direction设为可调整的参数，方便后续扩展。

使用示例 #

# 定义重采样规则
splitRule <- trainControl(method="repeatedcv", number=10, repeats=5, verboseIter = FALSE)

# 训练模型
model <- caret::train(
  OUTCOME ~ PredictorVariable1 + PredictorVariable2 + PredictorVariable3 + 
            PredictorVariable1 * PredictorVariable2,
  data = dataset$trainingSet,
  trControl = splitRule,
  na.action = na.omit,
  method = hierarchicalStepAIC,
  direction = "both"  # 传递给stepAIC的direction参数
)

# 预测
predictions_model <- predict(model, newdata = dataset$testingSet)

# 查看模型结果
print(model)

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