R语言循环遍历数据框计算分位数识别异常值的列引用问题

示例数据结构 #

现有如下结构的示例数据框：

structure(list(peso1_v00 = structure(c(97.8, 102.9, 78.4, 80.4, 
82.5, 92.6, 99, 66.8, 91.8, 70), label = "Peso: 1a determinación", format.spss = "F5.1"), 
    cintura1_v00 = structure(c(120, 123.2, 104, 106, 117.7, 115, 
    123, 93.5, 116.7, 109), label = "Cintura: 1a determinación", format.spss = "F5.1"), 
    tasis2_e_v00 = structure(c(139, 123, 138, 143, 160, 167, 
    139, 134, 161, 145), label = "TA: tensión arterial 2: sistólica", format.spss = "F4.0"), 
    tadias2_e_v00 = structure(c(74, 65, 75, 55, 75, 76, 65, 64, 
    62, 79), label = "TA: tensión arterial 2: diastólica", format.spss = "F4.0"), 
    p17_total_v00 = structure(c(7, 5, 9, 5, 8, 9, 4, 10, 10, 
    10), label = "Cuestionario de 17 puntos: Suma de puntuación de P17", format.spss = "F3.0"), 
    geaf_tot_v00 = structure(c(0, 1286.71, 5524.48, 419.58, 198.14, 
    839.16, 3188.81, 993.94, 2386.01, 3636.36), label = "AF: Gasto energético en actividad física total (MET•min/sem)", format.spss = "F8.2"), 
    glucosa_v00 = structure(c(116, 152, 113, 257, 108, 119, 108, 
    112, 141, 102), label = "Analítica: Glucosa en mg/dL", format.spss = "F4.0"), 
    albumi_v00 = structure(c(4.57, 4.71, 4.85, 4.43, 4.44, 4.59, 
    4.6, 4.56, 4.45, 4.59), label = "Analítica: Albúmina en g/dL", format.spss = "F6.2"), 
    coltot_v00 = structure(c(261, 197, 235, 168, 217, 248, 155, 
    254, 265, 326), label = "Analítica: Colesterol total en mg/dL", format.spss = "F4.0"), 
    hdl_v00 = structure(c(71, 43, 38, 36, 57, 49, 54, 60, 48, 
    60), label = "Analítica: Colesterol HDL en mg/dL", format.spss = "F4.0"), 
    ldl_calc_v00 = structure(c(147, 135, 161, NA, 145, 155, 73, 
    157, 172, 222), label = "Analítica: LDL calculado en mg/dL si trigli<=300", format.spss = "F4.0"), 
    trigli_v00 = structure(c(217, 93, 179, 315, 73, 219, 138, 
    186, 223, 220), label = "Analítica: Triglicéridos en mg/dL", format.spss = "F5.0"), 
    hba1c_v00 = structure(c(6.29, NA, 5.63, 12.17, NA, 6.11, 
    5.9, 5.68, NA, 6.16), label = "Analítica: Hemoglobina glicosilada (HbA1c %)", format.spss = "F5.2"), 
    i_hucpeptide_v00 = structure(c(854.34, NA, 1485.59, 4241.95, 
    NA, 847.89, 1524.39, 1265.48, NA, 290.22), label = "Hu C-peptide (72) IMIM S'han substituit en les següents var els codis de inf i sup a limit de detecció per el limit inf i sup de detecció", format.spss = "F9.2", display_width = 13L), 
    i_hughrelin_v00 = structure(c(681.94, NA, 480.11, 1587.73, 
    NA, 453, 263.93, 392.98, NA, 1327.91), label = "Hu Ghrelin (26) IMIM", format.spss = "F7.2", display_width = 10L), 
    i_hugip_v00 = structure(c(2.67, NA, 2.67, 2.67, NA, 2.67, 
    2.67, 2.67, NA, 2.67), label = "Hu GIP (14) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huglp1_v00 = structure(c(177.37, NA, 202.62, 519.02, NA, 
    200.13, 163.82, 20.29, NA, 14.14), label = "Hu GLP-1 (27) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huglucagon_v00 = structure(c(387.58, NA, 591.73, 855.29, 
    NA, 726.99, 430.35, 389.59, NA, 336.04), label = "Hu Glucagon (15) IMIM", format.spss = "F9.2", display_width = 11L), 
    i_huinsulin_v00 = structure(c(278.72, NA, 538.29, 1693.25, 
    NA, 299.75, 608.35, 397.7, NA, 129.17), label = "Hu Insulin (12) IMIM", format.spss = "F7.2", display_width = 10L), 
    i_huleptin_v00 = structure(c(2518.28, NA, 12175.88, 12369.5, 
    NA, 8409.76, 5998.71, 9298.52, NA, 5919.57), label = "Hu Leptin (78) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_hupai1_v00 = structure(c(3084.08, NA, 2650.85, 3202.18, 
    NA, 3085.25, 3410.73, 3109.79, NA, 1375.07), label = "Hu PAI-1 (61) IMIM", format.spss = "F7.2"), 
    i_huresistin_v00 = structure(c(4758.94, NA, 3594.11, 13564.63, 
    NA, 3221.72, 2864.01, 3630.63, NA, 2827.01), label = "Hu Resistin (65) IMIM", format.spss = "F8.2", display_width = 9L), 
    i_huvisfatin_v00 = structure(c(8.64, NA, 2081.59, 2363.58, 
    NA, 2989.72, 653.96, 1129.24, NA, 631.11), label = "Hu Visfatin (22) IMIM", format.spss = "F9.2", display_width = 6L), 
    col_rema_v00 = structure(c(43, 19, 36, NA, 15, 44, 28, 37, 
    45, 44), format.spss = "F8.2", display_width = 14L), homa_v00 = structure(c(1436.96, 
    NA, 2703.41, 19340.68, NA, 1585.34, 2920.08, 1979.66, NA, 
    585.57), format.spss = "F8.2", display_width = 10L), ln_trigli_v00 = structure(c(5.38, 
    4.53, 5.19, 5.75, 4.29, 5.39, 4.93, 5.23, 5.41, 5.39), format.spss = "F8.2", display_width = 15L), 
    ln_homa_v00 = structure(c(7.27, NA, 7.9, 9.87, NA, 7.37, 
    7.98, 7.59, NA, 6.37), format.spss = "F8.2", display_width = 13L), 
    ln_hba1c_v00 = structure(c(1.84, NA, 1.73, 2.5, NA, 1.81, 
    1.77, 1.74, NA, 1.82), format.spss = "F8.2", display_width = 14L), 
    ln_geaf_tot_v00 = structure(c(NA, 7.16, 8.62, 6.04, 5.29, 
    6.73, 8.07, 6.9, 7.78, 8.2), format.spss = "F8.2", display_width = 17L)), row.names = c(NA, 
-10L), class = "data.frame")

原始错误代码与现象 #

用户定义的待处理列列表和原始循环代码如下：

colvars <- c("peso1_v00", "cintura1_v00", "tasis2_e_v00" , "tadias2_e_v00" ,   "p17_total_v00", "geaf_tot_v00"  ,   "glucosa_v00", "albumi_v00", "coltot_v00", "hdl_v00", "ldl_calc_v00", "trigli_v00", "hba1c_v00", "i_hucpeptide_v00", "i_hughrelin_v00", "i_hugip_v00", "i_huglp1_v00", "i_huglucagon_v00", "i_huinsulin_v00", "i_huleptin_v00", "i_hupai1_v00", "i_huresistin_v00", "i_huvisfatin_v00", "col_rema_v00", "homa_v00",  "ln_trigli_v00", "ln_homa_v00", "ln_hba1c_v00", "ln_geaf_tot_v00")

i=1
for (i in 1:length(colvars)) {
  
  q<-sapply(outliers_v00[ ,i], quantile, probs = c(0.25, 0.75), na.rm =T)
  iqr<-sapply(outliers_v00[ ,i], IQR, na.rm =T)
  up <-  q[ , i[2]]+1.5*iqr # Upper Range  
  low<- q[ ,i[1]]-1.5*iqr # Lower Range
  remained<- subset(outliers, outliers[ ,i] > (q[1] - 1.5*iqr) & outliers[ ,i] < (q[2]+1.5*iqr))
  eliminated<-subset(outliers_v00, outliers_v00[,i] < (q[1] - 1.5*iqr) & outliers_v00[,i] > (q[2]+1.5*iqr))

}

哪怕简化为仅计算分位数和IQR的版本，也无法得到正确结果：

i=1
for (i in 1:length(colvars)) {
  
  q<-sapply(outliers_v00[ ,i], quantile, probs = c(0.25, 0.75), na.rm =T)
  iqr<-sapply(outliers_v00[ ,i], IQR, na.rm =T)
}

问题表现为：运行后输出的Q25和Q75数值相同，和实际计算结果不符。单独计算对应列的分位数可以正常得到结果，但后续筛选留存样本（remained）和剔除异常样本（eliminated）的逻辑无法在循环中正常运行。

错误原因与修正方案 #

核心错误点 #

• 无需对单列调用sapply：outliers_v00[,i]本身就是向量，直接传入quantile和IQR函数即可，多余的sapply会导致返回值索引逻辑混乱
• 分位数索引逻辑错误：计算得到的q是长度为2的向量，直接用q[1]取25分位、q[2]取75分位即可，q[, i[2]]属于完全错误的索引写法
• 异常值筛选逻辑错误：值 < 下限 & 值 > 上限是不可能成立的逻辑，应该用|（或）连接两个判断条件，同时需要额外处理NA值
• 变量名混用：代码中交替使用outliers和outliers_v00两个变量名，会导致对象找不到或者计算逻辑错误
• 结果覆盖问题：循环内反复给remained和eliminated赋值，每次循环都会把上一次的结果覆盖，无法保留多列的处理结果

循环中用[,i]引用列的正确写法 #

如果要按列存储每列的异常/正常样本，参考如下修正代码：

# 提前初始化列表存储每列的处理结果
remained_list <- list()
eliminated_list <- list()

for (i in 1:length(colvars)) {
  # 直接取当前列的向量，加drop=TRUE避免tibble类数据框返回单列数据框
  current_col <- outliers_v00[, colvars[i], drop = TRUE]
  # 计算分位数和IQR
  q <- quantile(current_col, probs = c(0.25, 0.75), na.rm = TRUE)
  iqr <- IQR(current_col, na.rm = TRUE)
  up <-  q[2] + 1.5*iqr # 上限
  low <- q[1] - 1.5*iqr # 下限
  
  # 生成当前列的正常样本掩码：数值在范围内且不为NA
  col_mask <- current_col >= low & current_col <= up & !is.na(current_col)
  # 存储当前列的处理结果
  remained_list[[colvars[i]]] <- outliers_v00[col_mask, ]
  eliminated_list[[colvars[i]]] <- outliers_v00[!col_mask