嘿，这个需求我太熟了！手动写几十个交互组合模型简直是噩梦，给你一套R的自动化方案，轻松搞定所有模型的拟合和指标对比～

核心思路 #

我们先自动生成所有可能的双向交互项组合，再把这些组合和基础变量拼接成模型公式，循环拟合GLM后提取你需要的AIC、BIC、R²指标，最后整理成表格方便对比。

具体实现代码 #

假设你的数据集叫claim_data，响应变量是ClaimAmount（记得替换成你实际的响应变量名），基础变量是你提到的4个分类变量：

# 加载必要的包（如果没装先运行install.packages(c("dplyr","purrr"))安装）
library(dplyr)
library(purrr)

# 1. 定义变量
response_var <- "ClaimAmount"  # 替换为你的响应变量名
base_vars <- c("Gender", "Age", "Nationality", "RoomType")

# 2. 生成所有双向交互项
all_interactions <- combn(base_vars, 2, FUN = function(x) paste(x, collapse = "*"))

# 3. 生成所有非空的交互项子集（单个、两个...直到所有交互）
subset_indices <- unlist(lapply(1:length(all_interactions), function(k) {
  combn(1:length(all_interactions), k, simplify = FALSE)
}), recursive = FALSE)

# 4. 构建所有模型公式（包含基础模型）
# 先做基础模型公式
base_formula <- as.formula(paste(response_var, "~", paste(base_vars, collapse = " + ")))
# 生成带交互项的所有公式
interaction_formulas <- map(subset_indices, ~ as.formula(paste(
  response_var, "~", 
  paste(c(base_vars, all_interactions[.x]), collapse = " + ")
)))
# 合并所有公式
all_formulas <- c(base_formula, interaction_formulas)

# 5. 批量拟合模型并提取指标
model_results <- map_dfr(all_formulas, function(formula) {
  # 拟合GLM，注意family要根据你的数据类型调整（比如伽马、泊松）
  model <- glm(formula, data = claim_data, family = gaussian())
  
  # 提取交互项名称（方便查看）
  interaction_terms <- setdiff(all.vars(formula[[3]]), base_vars)
  interaction_str <- ifelse(length(interaction_terms) == 0, "None", paste(interaction_terms, collapse = ", "))
  
  # 返回结果
  tibble(
    Model = deparse(formula),  # 用公式作为模型名，直观易懂
    AIC = AIC(model),
    BIC = BIC(model),
    R_squared = summary(model)$r.squared,
    Included_Interactions = interaction_str
  )
})

# 6. 按AIC排序（AIC越小模型越优）
sorted_results <- model_results %>% arrange(AIC)

# 查看最终结果
print(sorted_results)

关键注意事项 #

• Family参数调整：如果你的索赔金额是正偏分布（通常都是），建议把family改成Gamma(link = "log")；如果是计数型索赔次数，用poisson()或者quasipoisson()，根据你的数据类型调整。
• 模型名称可读性：代码里用deparse(formula)直接把公式作为模型名称，能一眼看到模型包含的所有项，比编号更直观。
• 效率优化：如果数据集特别大，拟合所有模型可能需要点时间，但绝对比手动写模型高效N倍；也可以用parallel包做并行拟合，进一步提速。
• 扩展功能：如果需要查看交互项的显著性，可以在结果里加上summary(model)$coefficients里的p值，不过你要的AIC/BIC/R²已经足够做模型对比了。

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