在R中运行多组Poisson回归时的对比错误排查

你的推测方向是对的，问题确实出在ethgp上，但不是全局的因子水平问题——而是在group_by(time, gender, age)拆分出的部分子组里，ethgp只有单一水平。虽然全局来看ethgp有M和NM两个水平，但当你按年份、性别、年龄分组后，某些细分的子组可能只包含其中一个族群的观测（比如2010年40-44岁男性组里全是M族群，没有NM的数据），这时候glm就没法完成族群间的对比计算，直接抛出了对比错误。

下面是具体的验证方法和解决思路：

1. 先确认哪些分组出了问题 #

你可以用这段代码快速定位所有ethgp水平数不足2的子组：

df %>% 
  group_by(time, gender, age) %>% 
  summarise(ethgp_unique = n_distinct(ethgp), .groups = "drop") %>% 
  filter(ethgp_unique < 2)

运行后会列出所有符合条件的子组，比如某年份+性别+年龄组合下只有M族群的数据，这些就是导致报错的根源。

2. 解决办法：根据需求选择处理方式 #

方式一：过滤掉无效子组（推荐） #

如果这些只有单一族群的子组对你的分析没有价值，可以直接在建模前过滤掉，只保留包含两个族群的子组：

Poisson_test <- df %>% 
  group_by(time, gender, age) %>% 
  filter(n_distinct(ethgp) >= 2) %>%  # 仅保留有完整族群数据的子组
  do({
    model = glm(asr_round ~ relevel(ethgp, ref = 2), family = "poisson", data = .)
    data.frame(nlRR_MNM = coef(model)[[2]], SE_MNM = summary(model)$coefficients[,2][2])
  })

方式二：保留所有子组并标记缺失值 #

如果你需要保留所有分组（比如后续要做可视化或补全），可以在do函数里加入判断逻辑，给单水平子组返回NA或自定义标记：

Poisson_test <- df %>% 
  group_by(time, gender, age) %>% 
  do({
    # 检查当前子组的ethgp水平数
    if(n_distinct(.$ethgp) < 2){
      data.frame(nlRR_MNM = NA, SE_MNM = NA)
    } else {
      model = glm(asr_round ~ relevel(ethgp, ref = 2), family = "poisson", data = .)
      data.frame(nlRR_MNM = coef(model)[[2]], SE_MNM = summary(model)$coefficients[,2][2])
    }
  })

补充说明 #

为什么全局检查ethgp是2水平还会报错？因为group_by之后，每个子组都是独立的小数据集，glm是针对每个小数据集单独建模的——只有当某个子组里同时存在M和NM两个族群时，才能计算两者的率比和置信区间，单一水平的子组没有对比的基础，自然会触发对比错误。

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