滚动窗口多列解释变量LME模型构建问题：优化、存储与报错排查

1. 公式报错的解决方法 (invalid model formula) #

你的报错根源在于外层循环遍历的是列值而非列名：

• 原代码中 ind_variables <- data[c("value1","value2","value3")] 是数据框，循环时 col_name 会变成列的数值向量，导致 paste("rwi ~", scale(col_name)) 生成的公式包含一堆数字，完全不符合语法。
• 修正方式：将 ind_variables 改为列名字符向量，同时正确构造公式（把 scale() 作用于列名而非列值）：

ind_variables <- c("value1", "value2", "value3") # 改为列名向量
# 外层循环内的公式构造
model_formula <- formula(paste("rwi ~ scale(", col_name, ")"))

2. 将结果存储为数据框列表 #

初始化一个空列表，在每个自变量的循环结束后，将该变量对应的结果数据框存入列表，并用自变量名作为列表元素的名称：

# 初始化结果列表
results_list <- list()

# 外层循环内，替换原有的rbind逻辑
for (col_name in ind_variables) {
  # 预分配当前变量的结果数据框（比反复rbind高效）
  window_count <- length(seq(1950, 2020 - window_size + 1, by = step_size))
  var_results <- data.frame(
    start_year = integer(window_count),
    end_year = integer(window_count),
    estimate = numeric(window_count),
    se = numeric(window_count),
    p_value = numeric(window_count),
    variable = rep(col_name, window_count) # 可选：标记自变量
  )
  
  # 内层循环填充数据
  row_idx <- 1
  for (start_year in seq(1950, 2020 - window_size + 1, by = step_size)) {
    end_year <- start_year + window_size - 1
    subset_data <- data %>% filter(year >= start_year & year <= end_year)
    
    model <- lme(fixed = model_formula, random = ~1|treeID, 
                 data = subset_data, na.action = na.exclude, method = "ML")
    
    # 提取自变量的系数（排除截距项，取第2行）
    t_table <- summary(model)$tTable[2, ]
    var_results[row_idx, ] <- list(
      start_year = start_year,
      end_year = end_year,
      estimate = t_table["Value"],
      se = t_table["Std.Error"],
      p_value = t_table["p-value"],
      variable = col_name
    )
    row_idx <- row_idx + 1
  }
  
  # 将当前变量的结果存入列表
  results_list[[col_name]] <- var_results
}

之后你可以通过 results_list$value1、results_list$value2 分别访问每个自变量的结果。

3. 更高效的实现方式 #

核心优化点： #

1. 预分配数据框：避免循环中反复 rbind（每次 rbind 都会复制整个数据框，数据量大时极慢）。
2. 用 lapply 替代外层循环：更符合R的函数式编程风格，减少副作用：

# 定义单个自变量的处理函数
process_variable <- function(col_name) {
  window_seq <- seq(1950, 2020 - window_size + 1, by = step_size)
  window_count <- length(window_seq)
  
  var_results <- data.frame(
    start_year = integer(window_count),
    end_year = integer(window_count),
    estimate = numeric(window_count),
    se = numeric(window_count),
    p_value = numeric(window_count),
    variable = rep(col_name, window_count)
  )
  
  for (i in seq_along(window_seq)) {
    start_year <- window_seq[i]
    end_year <- start_year + window_size - 1
    subset_data <- data[data$year >= start_year & data$year <= end_year, ] # 用base R筛选比dplyr更快
    
    model <- lme(fixed = formula(paste("rwi ~ scale(", col_name, ")")), 
                 random = ~1|treeID, data = subset_data, 
                 na.action = na.exclude, method = "ML")
    
    t_table <- summary(model)$tTable[2, ]
    var_results[i, ] <- list(
      start_year = start_year,
      end_year = end_year,
      estimate = t_table["Value"],
      se = t_table["Std.Error"],
      p_value = t_table["p-value"],
      variable = col_name
    )
  }
  return(var_results)
}

# 批量处理所有自变量
results_list <- lapply(ind_variables, process_variable)
names(results_list) <- ind_variables # 给列表元素命名

3. 可选优化：如果数据量极大，可改用 data.table 进行窗口筛选，速度会比 dplyr 或base R更快；也可以用 parallel::mclapply 并行处理多个自变量（适合多核机器）。

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