如何在R中使用rollapply()实现滚动窗口极大似然估计(MLE)

你遇到的报错核心是mle2()拟合失败返回了空系数，触发了类有效性校验，具体有4个可修复的错误点：

• 负似然函数硬编码调用全局数据集CFbasin2，没有使用当前滚动窗口的切片数据，相当于每次拟合都用全量3768行数据而非20行窗口数据，同时存在列名不匹配问题：你的数据集列名为Qdeep，但代码中提取的是Qdp
• 负似然函数中对误差项均值做无约束估计，和模型自带的截距项b0完全共线性，导致参数不可识别，拟合无法收敛
• 滚动窗口内没有定义初始值dw_coef、up_coef，不合理的初始值会直接导致拟合失败返回空系数
• rollapply默认要求返回可拼接的向量/矩阵结果，直接返回mle2的S4对象无法完成后续的结果合并

第一步：重写负似然函数 #

把窗口数据作为入参，删除冗余的误差均值估计，将误差标准差作为待估参数：

regress.ll <- function(b0, b1, r0, r1, r2, b2, sd_err, window_data) {
  x1 <- window_data$Qhaw
  x2 <- window_data$Qdeep
  x3 <- window_data$Qlit
  y <- window_data$Qkl
  y.hat <- b0 + b1*(r0 + r1*x1 + r2*x2) + b2*x3
  errors <- y - y.hat
  neg.loklik <- -sum(dnorm(errors, mean = 0, sd = sd_err, log = TRUE))
  return(neg.loklik)
}

第二步：改写滚动窗口逻辑 #

每个窗口先做普通线性回归得到合理初始值，拟合成功后提取系数返回：

library(zoo)
library(bbmle)

roll_mle_result <- rollapply(CFbasin2,
          width = 20,
          FUN = function(Z) {
            # 把当前窗口转成数据框，保留列名
            window_df <- as.data.frame(Z)
            # 先做线性回归获取初始值：适配你的模型形式拆解得到线性回归系数
            lm_fit <- lm(Qkl ~ Qhaw + Qdeep + Qlit, data = window_df)
            lm_coef <- coef(lm_fit)
            # 给待估参数设置合理初始值，避免收敛失败
            start_vals <- list(
              b0 = lm_coef[1]/2, 
              b1 = 1, 
              r0 = lm_coef[1]/2, 
              r1 = lm_coef[2], 
              r2 = lm_coef[3], 
              b2 = lm_coef[4],
              sd_err = sigma(lm_fit)
            )
            # 拟合MLE，加异常捕获避免个别窗口收敛失败中断整个流程
            mle.results <- tryCatch({
              mle2(minuslogl = regress.ll, 
                   data = list(window_data = window_df),
                   start = start_vals,
                   control = list(maxit = 1000, trace = FALSE))
            }, error = function(e) {
              return(rep(NA, 7))
            })
            # 提取拟合系数返回，适配rollapply的输出格式
            if(inherits(mle.results, "mle2")) {
              return(coef(mle.results))
            } else {
              return(mle.results)
            }
          }, by.column=FALSE, align="right")

结果说明 #

输出的roll_mle_result是一个矩阵，每一行对应一个窗口的7个待估参数（b0、b1、r0、r1、r2、b2、sd_err），收敛失败的窗口会填充NA。如果需要提取每个窗口的完整mle2对象，可以用lapply遍历窗口索引实现。

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