R循环调用线性模型残差方差检验函数列值不匹配问题

我开发了一个可对线性模型残差生成方差检验结果或诊断摘要的自定义函数：

• 方差检验模式下，需要传入数据集、待检验列名、分组阈值三类输入参数
• 遍历线性模型列表循环调用该函数时，输出不符合预期，判断为传入的列名与对应阈值没有正确匹配

原始问题代码 #

perTest_nl <- function(model, ...) {
  pr <- list(...)
    if (is(pr[[1]],'list') == TRUE) {
      if (length(pr[[1]]) != 3) {
        stop('You need three parameters, i.e. list("data", "column", "value")')
      } else {
        dataset <- pr[[1]][[1]] 
        column <- pr[[1]][[2]]
        value <- pr[[1]][[3]]
        dataset <- dataset %>% data.frame()
        variance_test <-
            var.test(residuals(model)[dataset[,column] > value], residuals(model)[dataset[,column] < value])
        return(variance_test)
      }
    } else {
      sumry <-
        parse(text = "summary(lm(sqrt(abs(residuals(model))) ~ fitted(model)))")
      
      if (pr[[1]] == 'summary') {
        model_diagnostic <-
            eval(sumry)$coefficients %>%
              data.frame() %>%
              .[2, ] %>%
              add_column(
                F_statistic = eval(sumry)$fstatistic[1],
                p_value = 1 - pf(
                  eval(sumry)$fstatistic[[1]],
                  eval(sumry)$fstatistic[[2]],
                  eval(sumry)$fstatistic[[3]]
                ),
                df = eval(sumry)$fstatistic[[3]],
                RSE = eval(sumry)$sigma
              ) %>% `colnames<-`(
                c(
                  'Estimate',
                  'Std.Error',
                  't.value',
                  'Pr(>|t|)',
                  'F.statistic',
                  'P.value',
                  'df',
                  'RSE'
                )
              )
        return(model_diagnostic)
      } 
    }
  }  

library(faraway)
data(savings)
library(tidyverse)

predictors <- c("pop15", "pop75", "dpi", "ddpi")
scan_values <- list(pop15 = 35, pop75 = 2.5, dpi = 2000, ddpi = 7)

lm1 <- lm(sr ~ pop15 +pop75+dpi, savings)
lm2 <- lm(sr ~ pop15 +pop75+dpi+ddpi, savings)
lmod1<- list(lm1, lm2)

val_list <- c()
for(i in predictors){
    val_list[[i]]<-map(lmod1, function(x)perTest_nl(x, list(savings, i, scan_values)))
}

异常表现 #

循环运行得到的ddpi分组下第二个模型检验结果：

F = 0.96324, num df = 9, denom df = 39, p-value = 0.9685

手动调用perTest_nl(lmod1[[2]], list(savings, 'ddpi', 7))得到的正确结果：

F = 0.7781, num df = 5, denom df = 43, p-value = 0.8581

显式传入scan_values最后一组键值（ddpi=7）遍历模型时，结果和手动调用一致，说明循环迭代中参数传递存在匹配错误。

核心问题有两个：

1. 传入阈值参数时错误传入了整个scan_values列表，而非当前循环变量i对应的单个阈值。手动调用时第三个参数是数值7，但循环中传入的是长度为4的完整阈值列表，导致子集筛选的比较逻辑完全失效。
2. 结果列表初始化用c()不够严谨，虽然不会直接导致计算错误，但在列表嵌套赋值时容易出现类型转换异常。

只需要修改循环部分的参数传入逻辑，将第三个参数从完整的scan_values改为当前预测变量对应的阈值scan_values[[i]]，同时将结果列表初始化为空列表list()即可：

# 初始化空列表存储结果
val_list <- list()
for(i in predictors){
    # 第三个参数传入当前i对应的单个阈值，而非整个列表
    val_list[[i]] <- map(lmod1, function(x) perTest_nl(x, list(savings, i, scan_values[[i]])))
}

如果遇到旧版本purrr的惰性求值捕获问题，可以把匿名函数改成强制传参的写法，彻底避免环境绑定错误：

val_list <- list()
for(i in predictors){
    cur_col <- i
    cur_val <- scan_values[[i]]
    val_list[[i]] <- map(lmod1, function(x) perTest_nl(x, list(savings, cur_col, cur_val)))
}

修改后运行val_list[['ddpi']][[2]]，即可得到和手动调用完全一致的检验结果。

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