在R中对含缺失值的多列执行滚动回归的实现需求

需求说明 #

• 对数据框y的每一列与x执行窗口宽度为12的滚动回归
• 需计算每个窗口的回归系数（Beta）、截距（Intercept）和R²值
• 要求每个窗口至少包含12个有效观测值，不满足条件的位置填充NA
• 最终输出三个独立的结果矩阵/数据框：
  • 存储Beta值的矩阵，列名对应y的列（y1、y2、y3），默认尺寸为53×3（若需3×53可转置）
  • 存储截距值的矩阵，规格同上
  • 存储R²值的矩阵，规格同上
• 实际场景中y包含2000+列，需避免显式循环

示例数据 #

# 构造示例数据
y <- data.frame(y1 = c(NA,NA,NA, rnorm(50)),
                y2 = c(NA, NA, NA, NA, NA, rnorm(48)),
                y3 = rnorm(53))
x <- data.frame(x = c(NA, NA,NA, NA, NA, rnorm(48)))

尝试代码的问题 #

你提供的rollapply代码存在以下问题：

• FUN参数写法错误：rollapply的自定义函数仅接受一个参数，无法传入a,b两个参数
• 未提取目标统计量：直接返回lm对象无法得到可直接使用的数值结果
• 未校验窗口内有效观测数：无法保证每个窗口至少有12个有效样本

解决方案代码 #

library(zoo)

# 定义滚动回归统计量计算函数
roll_reg_stats <- function(y_col, x_col, window_width = 12) {
  # 单窗口处理逻辑
  single_window <- function(idx) {
    # 提取当前窗口的y和x值
    y_win <- y_col[idx]
    x_win <- x_col[idx]
    # 统计有效观测数（无NA的配对样本）
    valid_count <- sum(!is.na(y_win) & !is.na(x_win))
    
    # 有效数不足则返回NA
    if (valid_count < window_width) {
      return(c(Intercept = NA, Beta = NA, R2 = NA))
    }
    
    # 拟合回归模型（自动排除NA）
    model <- lm(y_win ~ x_win, na.action = na.exclude)
    # 提取截距、Beta和R²
    model_coefs <- coef(model)
    model_r2 <- summary(model)$r.squared
    
    return(c(Intercept = model_coefs[1], Beta = model_coefs[2], R2 = model_r2))
  }
  
  # 生成右对齐的滚动窗口索引
  window_indices <- rollapply(seq_along(y_col), 
                             width = window_width, 
                             FUN = identity, 
                             align = "right", 
                             fill = NA)
  
  # 对每个窗口应用处理逻辑
  results <- t(apply(window_indices, 1, function(idx) {
    if (any(is.na(idx))) return(c(NA, NA, NA))
    single_window(idx)
  }))
  
  return(results)
}

# 将x转换为向量便于处理
x_vector <- x$x

# 批量处理y的所有列（隐式循环，效率远高于显式for循环）
all_results <- lapply(y, function(col) roll_reg_stats(col, x_vector))

# 拆分并整理三个结果矩阵
beta_matrix <- do.call(cbind, lapply(all_results, function(res) res[, "Beta"]))
intercept_matrix <- do.call(cbind, lapply(all_results, function(res) res[, "Intercept"]))
r2_matrix <- do.call(cbind, lapply(all_results, function(res) res[, "R2"]))

# 设置列名与y保持一致
colnames(beta_matrix) <- colnames(y)
colnames(intercept_matrix) <- colnames(y)
colnames(r2_matrix) <- colnames(y)

# 若需3×53的矩阵（行对应统计量，列对应观测），可执行转置：
# beta_matrix <- t(beta_matrix)
# intercept_matrix <- t(intercept_matrix)
# r2_matrix <- t(r2_matrix)

方案说明 #

• 采用lapply批量处理y的列，避免显式for循环，适合2000+列的大规模数据
• 内置窗口有效性校验，确保只有有效观测数≥12的窗口才会拟合模型
• 直接提取所需统计量，输出结果为数值矩阵，便于后续分析
• 基于zoo包的rollapply实现滚动窗口逻辑，稳定性和效率有保障

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