R语言扩展for循环实现多变量组合GAM模型报错排查

问题背景 #

计划开发循环批量运行多个GAM模型，生成包含AIC、平滑项p值、解释偏差等指标的排序结果表。已适配单个自变量的运行逻辑，现扩展至所有自变量组合（从1个到全部协变量），尝试双变量组合时触发错误：Error in all.vars1(gp$fake.formula[-2]) : can't handle [[ in formula

原始代码 #

library(mgcv)

## Example data
set.seed(0) 
dat <- gamSim(1,n=200,scale=2)
set.seed(1)
dat2 <- gamSim(1,n=200,scale=2)
names(dat2)[1:5] <- c("y1", paste0("x", 4:7))
d <- cbind(dat[, 1:5], dat2[, 1:5])

d_resp <- d[ c("y", "y1")]
d_pred <- d[, !(colnames(d) %in% c("y", "y1"))]

df_combinations2 <- lapply(1:(ncol(combn(1:ncol(d_pred), m = 2))), 
                           function(y) d_pred[, combn(1:ncol(d_pred), m = 2)[,y]]) 

## create a "matrix" list of dimensions i x j
results_m2 <-lapply(1:length(df_combinations2), matrix, data= NA, nrow=ncol(d_resp), ncol=2)

## for-loop
for(k in 1:length(df_combinations2)){
  for(i in 1:ncol(d_resp)){
  for(j in 1:ncol(df_combinations2[[k]])){
    results_m2[i, j][[1]] <- gam(d_resp[, i] ~ s(df_combinations2[[k]][,1])+s(df_combinations2[[k]][,2]))
  }
}}

错误原因 #

1. 公式解析失败：mgcv::gam()无法识别公式中的[[索引语法——直接在公式里使用df_combinations2[[k]][,1]这类动态对象索引，GAM的公式解析器需要明确的变量名，无法处理此类引用。
2. 循环逻辑冗余：内层j循环完全多余，每个双变量组合仅需运行一次包含两个平滑项的模型，无需遍历单个自变量。

修复方案 #

核心思路 #

• 放弃直接在公式中使用列表索引，改为动态构建字符串形式的公式，再转换为公式对象；
• 简化循环逻辑，每个因变量+双变量组合对应一个模型；
• 基于原始数据框d引用变量，避免拆分数据后的索引问题。

修复后的代码 #

library(mgcv)

## Example data
set.seed(0) 
dat <- gamSim(1,n=200,scale=2)
set.seed(1)
dat2 <- gamSim(1,n=200,scale=2)
names(dat2)[1:5] <- c("y1", paste0("x", 4:7))
d <- cbind(dat[, 1:5], dat2[, 1:5])

d_resp <- c("y", "y1")  # 直接存储变量名，而非数据框切片
d_pred <- colnames(d)[!colnames(d) %in% d_resp]

# 生成所有双变量组合的变量名对
combos2 <- combn(d_pred, m = 2, simplify = FALSE)

# 初始化结果列表：每个元素对应一个双变量组合，内部存储两个因变量的模型
results_m2 <- vector("list", length(combos2))
names(results_m2) <- sapply(combos2, paste, collapse = "+")

## 简化后的循环
for(k in seq_along(combos2)){
  vars <- combos2[[k]]
  # 动态构建公式字符串：因变量 ~ s(var1) + s(var2)
  for(resp in d_resp){
    formula_str <- paste(resp, "~", paste0("s(", vars, ")", collapse = "+"))
    model <- gam(as.formula(formula_str), data = d)
    results_m2[[k]][[resp]] <- model
  }
}

# 提取指标示例：提取每个模型的AIC和解释偏差
model_metrics <- lapply(results_m2, function(combo_models){
  lapply(combo_models, function(mod){
    data.frame(
      AIC = AIC(mod),
      dev_explained = 1 - mod$deviance/mod$null.deviance,
      smooth_pvals = sapply(mod$smooth, function(sm) sm$p.value)
    )
  })
})

扩展到全变量组合 #

若要覆盖1个到所有自变量的组合，可遍历组合长度实现：

# 生成所有长度1至最大变量数的组合
all_combos <- unlist(lapply(1:length(d_pred), function(m){
  combn(d_pred, m = m, simplify = FALSE)
}), recursive = FALSE)

# 初始化全组合结果列表
results_all <- vector("list", length(all_combos))
names(results_all) <- sapply(all_combos, paste, collapse = "+")

# 循环运行所有组合的模型
for(k in seq_along(all_combos)){
  vars <- all_combos[[k]]
  for(resp in d_resp){
    formula_str <- paste(resp, "~", paste0("s(", vars, ")", collapse = "+"))
    model <- gam(as.formula(formula_str), data = d)
    results_all[[k]][[resp]] <- model
  }
}

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