R语言基于group_by(.add=TRUE)的分位数嵌套分组计算正确性验证

我正在用R做数据分析，构建了模拟患者数据集：

set.seed(123)
library(dplyr)

Patient_ID = 1:5000
gender <- c("Male","Female")
gender <- sample(gender, 5000, replace=TRUE, prob=c(0.45, 0.55))
Gender <- as.factor(gender)


status <- c("Immigrant","Citizen")
status <- sample(status, 5000, replace=TRUE, prob=c(0.3, 0.7))
Status  <- as.factor(status )

Height = rnorm(5000, 150, 10)
Weight = rnorm(5000, 90, 10)
Hospital_Visits = sample.int(20,  5000, replace = TRUE)

################

disease <- c("Yes","No")
disease <- sample(disease, 5000, replace=TRUE, prob=c(0.4, 0.6))
Disease <- as.factor(disease)

###################
my_data = data.frame(Patient_ID, Gender, Status, Height, Weight, Hospital_Visits, Disease)

之前我用ntile实现了嵌套分组的疾病比例计算，代码如下：

# e.g. using 3 ntiles

my_data %>% 
  group_by(Gender, Status) %>%
  mutate(Height_ntile = ntile(Height, 3),
         Height_range = paste(min(Height), max(Height), sep = "-")) %>%
  group_by(Height_ntile, Height_range, .add = TRUE) %>%
  mutate(Weight_ntile = ntile(Weight, 3),
         Weight_range = paste(min(Weight), max(Weight), sep = "-")) %>%
  group_by(Weight_ntile, Weight_range, .add = TRUE) %>%
  mutate(Hospital_Visits_ntile = ntile(Hospital_Visits, 3),
         Hospital_range = paste(min(Hospital_Visits), max(Hospital_Visits), sep = "-")) %>%
  group_by(Hospital_Visits_ntile, Hospital_range, .add = TRUE) %>%
  summarize(percent_disease = mean(Disease == "Yes"), 
            count = n(),
            .groups = "drop")

现在我想改成用分位数实现，写了下面的代码，能运行但不确定是否正确——比如结果里出现了-Inf值，麻烦帮忙验证：

my_data %>% 
  group_by(Gender, Status) %>%
  mutate(Height_group = cut(Height, breaks = c(-Inf, 
                                               quantile(Height, c(0.33, 0.67)), 
                                               Inf)),
         Height_range = paste(min(Height), max(Height), sep = "-")) %>%
  group_by(Height_group, Height_range, .add = TRUE) %>%
  mutate(Weight_group = cut(Weight, breaks = c(-Inf, 
                                               quantile(Weight, c(0.33, 0.67)), 
                                               Inf)),
         Weight_range = paste(min(Weight), max(Weight), sep = "-")) %>%
  group_by(Weight_group, Weight_range, .add = TRUE) %>%
  mutate(Hospital_Visits_group = cut(Hospital_Visits, breaks = c(-Inf, 
                                                                quantile(Hospital_Visits, c(0.33, 0.67)), 
                                                                Inf)),
         Hospital_range = paste(min(Hospital_Visits), max(Hospital_Visits), sep = "-")) %>%
  group_by(Hospital_Visits_group, Hospital_range, .add = TRUE) %>%
  summarize(percent_disease = mean(Disease == "Yes"), 
            count = n(),
            .groups = "drop")

代码问题分析与修正 #

1. -Inf出现的原因 #

你用c(-Inf, quantile(...), Inf)作为cut的breaks时，cut默认会将左闭右开区间的左端点显示为-Inf，这是格式问题而非逻辑错误，但可以通过自定义标签优化输出。另外，若分组内的变量值完全重复，quantile会返回重复分位点，可能导致cut报错，需提前处理。

2. 分组逻辑的偏差与优化 #

原ntile代码的嵌套分组逻辑是正确的，但分位数实现有两处可优化：

• 用probs = c(1/3, 2/3)替代近似值0.33、0.67，确保和ntile的三分位逻辑完全对齐
• 增加分位点去重处理，避免cut因重复breaks报错
• 自定义分组标签，让结果更易读

修正后的代码 #

my_data %>% 
  group_by(Gender, Status) %>%
  mutate(
    # 计算准确三分位数，处理重复分位点，自定义分组标签
    Height_quantiles = quantile(Height, probs = c(1/3, 2/3)),
    Height_breaks = unique(c(-Inf, Height_quantiles, Inf)),
    Height_group = cut(Height, breaks = Height_breaks,
                       labels = c("Low", "Medium", "High")),
    Height_range = paste(round(min(Height), 2), round(max(Height), 2), sep = "-")
  ) %>%
  group_by(Height_group, Height_range, .add = TRUE) %>%
  mutate(
    Weight_quantiles = quantile(Weight, probs = c(1/3, 2/3)),
    Weight_breaks = unique(c(-Inf, Weight_quantiles, Inf)),
    Weight_group = cut(Weight, breaks = Weight_breaks,
                       labels = c("Low", "Medium", "High")),
    Weight_range = paste(round(min(Weight), 2), round(max(Weight), 2), sep = "-")
  ) %>%
  group_by(Weight_group, Weight_range, .add = TRUE) %>%
  mutate(
    Hospital_quantiles = quantile(Hospital_Visits, probs = c(1/3, 2/3)),
    Hospital_breaks = unique(c(-Inf, Hospital_quantiles, Inf)),
    Hospital_Visits_group = cut(Hospital_Visits, breaks = Hospital_breaks,
                                labels = c("Low", "Medium", "High")),
    Hospital_range = paste(min(Hospital_Visits), max(Hospital_Visits), sep = "-")
  ) %>%
  group_by(Hospital_Visits_group, Hospital_range, .add = TRUE) %>%
  summarize(
    percent_disease = mean(Disease == "Yes"), 
    count = n(),
    .groups = "drop"
  )

验证正确性的方法 #

对比原ntile代码和修正后的分位数代码结果：

1. 嵌套分组的样本量（count列）应大致接近——ntile是强制均分样本，分位数是按实际分布划分，样本量可能略有差异，但整体趋势一致
2. 同一嵌套分组下的percent_disease数值趋势应相同，说明分组逻辑的一致性

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