基于灵敏度、FPR、特异性统计变量选择算法的正确/过度设定回归模型数

我在k个合成数据集上分别运行了k次向后逐步回归（Backward Elimination, BE）和向前逐步回归（Forward Selection, FS），并编写代码计算了每次BE和FS回归所选变量的灵敏度（真阳性率，TPR）、假阳性率（FPR）和特异性（真阴性率，TNR）。针对BE逐步回归的实现代码如下：

### Benchmark 2: Run a Backward Elimination Stepwise Regression
### function on each of the  csvs.
### Assign the full models to their corresponding csvs and
### store these in the object "all_regressors_models"
library(parallel)
CL <- makeCluster(detectCores() - 1L)
clusterExport(CL, c('datasets'))
set.seed(11)      # for reproducibility
system.time(BE.fits <- parLapply(CL, datasets, \(X) {
    full_models <- lm(X$Y ~ ., X)
    back <- step(full_models, scope = formula(full_models), 
                    direction = 'back', trace = FALSE) }) )

BE_Coeffs <- lapply(seq_along(BE.fits), function(i) coef(BE.fits[[i]]))
stopCluster(CL)

IVs_Selected_by_BE <- lapply(seq_along(BE.fits), 
                             \(i) names(coef(BE.fits[[i]])[-1]))


### Count up how many Variables Selected match  the true 
### structural equation variables for that dataset in order
### to measure BE's performance.
# the True Positive Rate
Total_Positives <- lapply(True_Regressors, function(i) { length(i) })  
BE_TPs <- lapply(seq_along(datasets), \(i)
                                    sum(IVs_Selected_by_BE[[i]] %in% 
                                          True_Regressors[[i]]))
BM2_TPRs = lapply(seq_along(datasets), \(j)
                  j <- (BE_TPs[[j]]/Total_Positives[[j]]) )


# the False Positive Rate
BE_NNs <- lapply(True_Regressors, function(i) {30 - length(i)})
BE_FPs <- lapply(seq_along(datasets), \(i)
                          sum(!(IVs_Selected_by_BE[[i]] %in% 
                                True_Regressors[[i]]))) 
BE_FPRs = lapply(seq_along(datasets), \(j)
                  j <- (BM1_FPs[[j]])/BM1_NPs[[j]])


# the True Negative Rate
BE_TNRs <- lapply(BE_FPRs, \(i) 
                   i <- (1 - i))

注：代码中使用30 - length(i)的原因是每个数据集为500行31列，包含30个候选自变量列和1个因变量列。

需求 #

需要实现两个统计目标：

• 统计灵敏度=1且特异性=1的次数，对应BE选择的「正确设定」模型数量（无遗漏或冗余变量）；
• 统计灵敏度=1且FPR>0的次数，对应BE选择的「过度设定」模型数量（包含冗余变量但未遗漏核心变量）。

解决方案 #

第一步：修正原代码错误 #

原代码计算BE_FPRs时误用了BM1_FPs和BM1_NPs（应为BE相关变量），修正后的代码：

# 修正假阳性率计算
BE_FPRs <- lapply(seq_along(datasets), \(j) {
  BE_FPs[[j]] / BE_NNs[[j]]
})

第二步：将列表结果转为向量 #

为方便批量统计，先把列表格式的TPR、TNR、FPR转为向量：

# 列表转向量
BE_TPRs_vec <- unlist(BM2_TPRs)
BE_TNRs_vec <- unlist(BE_TNRs)
BE_FPRs_vec <- unlist(BE_FPRs)

第三步：执行统计 #

1. 统计正确设定模型数量

correct_spec_count <- sum(BE_TPRs_vec == 1 & BE_TNRs_vec == 1)

2. 统计过度设定模型数量

over_spec_count <- sum(BE_TPRs_vec == 1 & BE_FPRs_vec > 0)

上述代码通过向量化条件判断实现了类似ExcelSUMIF的统计逻辑，直接对所有数据集的结果进行批量计数。

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