循环构建线性回归模型：更换因变量与控制变量并纳入固定效应

第一步：修正并生成数据集 #

原代码中pivot_longer的匹配规则存在问题，starts_with("X")无法匹配年份列（2011/2012/2013），需改为starts_with("20")才能正确转换为长格式面板数据：

library(dplyr)
library(tidyr)

# 构造初始数据集
data <- data.frame(
  Country = c("USA","USA","USA","USA","India","India","India","India","China","China","China","China"),
  Indicator = rep(c("Population","GDP","Debt","Currency"), times=3),
  `2011` = rep(c(1,2,3,4), each=3),
  `2012` = rep(c(4,5,6,7), each=3),
  `2013` = rep(c(8,9,11,12), each=3)
)

# 转换为长格式并重塑宽格式
data <- data %>%
  pivot_longer(
    starts_with("20"),  # 修正匹配规则，选中年份列
    names_to = "Year",
    names_transform = list(Year = parse_number)
  ) %>%
  pivot_wider(names_from = Indicator, values_from = value) %>% 
  relocate(Year)

# 添加因变量y1-y5
data$y1 <- c(1,10,11,3,4,5,2,2,1)
data$y2 <- c(1,2,3,4,5,6,6,8,9)
data$y3 <- c(10,9,8,7,5,5,11,3,4)
data$y4 <- c(1,1,11,3,4,2,2,2,1)
data$y5 <- c(5,10,11,3,5,5,5,5,1)

第二步：批量构建模型并存储结果 #

使用plm包处理面板固定效应模型，同时规避共线性问题，循环遍历所有因变量、控制变量组合及固定效应选项：

library(plm)

# 定义因变量列表
y_vars <- paste0("y", 1:5)

# 定义控制变量候选集（核心自变量为GDP）
control_candidates <- c("Population", "Debt", "Currency", "Country")

# 生成所有控制变量组合（含无控制变量的情况）
control_combinations <- unlist(lapply(0:length(control_candidates), function(k) {
  combn(control_candidates, k, FUN = function(x) paste(x, collapse = " + "))
}), use.names = FALSE)
# 将空组合替换为"1"（表示仅保留截距项）
control_combinations[control_combinations == ""] <- "1"

# 定义固定效应选项
fe_options <- list(
  "无固定效应" = "",
  "仅国家固定效应" = "+ factor(Country)",
  "仅时间固定效应" = "+ factor(Year)",
  "双向固定效应" = "+ factor(Country) + factor(Year)"
)

# 创建嵌套列表存储所有模型摘要结果
model_results <- list()

# 循环拟合所有模型
for (y in y_vars) {
  model_results[[y]] <- list()
  for (controls in control_combinations) {
    # 构建基础公式（因变量~GDP+控制变量）
    base_formula <- if (controls == "1") {
      paste(y, "~ GDP")
    } else {
      paste(y, "~ GDP +", controls)
    }
    for (fe_label in names(fe_options)) {
      fe_term <- fe_options[[fe_label]]
      # 跳过控制变量含Country且同时加入国家固定效应的情况（避免完全共线性）
      if (grepl("Country", controls) && grepl("factor\\(Country\\)", fe_term)) {
        next
      }
      # 拼接完整公式
      full_formula <- as.formula(paste(base_formula, fe_term))
      
      # 拟合模型：有固定效应用plm的within模型，无固定效应用普通lm
      model <- if (fe_term != "") {
        plm(full_formula, data = data, model = "within", index = c("Country", "Year"))
      } else {
        lm(full_formula, data = data)
      }
      
      # 命名模型并存储摘要
      model_name <- paste0("控制变量: ", ifelse(controls == "1", "无", controls), " | ", fe_label)
      model_results[[y]][[model_name]] <- summary(model)
    }
  }
}

第三步：查看模型结果 #

通过嵌套列表索引查看指定模型的摘要，例如查看y1中仅GDP+双向固定效应的模型：

# 查看y1的目标模型结果
model_results[["y1"]][["控制变量: 无 | 双向固定效应"]]

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