在R的survey包中处理单元格计数为零的问题

问题背景 #

我正在学习R语言的survey包，需为58×4的列联表计算单元格比例（行百分比）和置信区间，同时考虑复杂抽样设计的方差估计。行对应加州58个县，列对应四类“性取向”变量。部分小县存在零单元格计数，使用svyby方法时触发错误：

data length is not a sub-multiple or multiple of the number of rows

可复现代码 #

# 设置试验次数和各结果概率
set.seed(16)
probs <- c(.94, .03, .02, .01)
n <- 10000

# 构造含零单元格的数据集
sogi <- data.frame(t(rmultinom(n, 1, prob = probs)))
names(sogi) <- c("straight","gay","bisexual","other")
df_sogi <- sogi %>% 
  mutate(sexual_orientation = 
           as.factor(case_when(
             straight == 1 ~ 0,
             gay == 1 ~ 1,
             bisexual == 1 ~ 2,
             other == 1 ~ 3))) %>%
  mutate(pw = runif(n, 0.2, 3.5)) %>%
  mutate(county = as.factor(rep(LETTERS[1:10], each = n/length(LETTERS[1:10])))) %>%
  mutate(ind = row_number()) %>%
  select(ind, county, sexual_orientation, pw) 

# 查看原始计数
table(df_sogi$county, df_sogi$sexual_orientation)

# 创建调查设计对象
df_svy_z1 <- svydesign(id=~ind, data = df_sogi, weights = ~pw)

# 估计单元格均值
svyby(formula = ~county, by = ~sexual_orientation, design = df_svy_z1, FUN = svymean)


# 修改数据集，构造零单元格
df_sogi$sexual_orientation <- ifelse(df_sogi$county %in% c("G","I") & df_sogi$sexual_orientation==3, 
                                     0, df_sogi$sexual_orientation)

# 查看修改后的计数
table(df_sogi$county, df_sogi$sexual_orientation)

df_svy_z2 <- svydesign(id=~ind, data = df_sogi, weights = ~pw)
# 运行此代码会触发报错
svyby(~factor(county), ~factor(sexual_orientation), df_svy_z2, svymean)

问题原因 #

当部分县在某类性取向中无观测值（零单元格）时，factor(sexual_orientation)在不同县的子集中会丢失对应水平，导致svyby合并结果时出现维度不匹配，触发报错。

解决方案 #

方案1：强制保留因子的完整水平 #

修改数据后，重新指定sexual_orientation的所有可能水平，确保每个子集的因子维度一致：

# 修改数据后重置因子水平，保留全部4类性取向
df_sogi$sexual_orientation <- factor(df_sogi$sexual_orientation, levels = c(0,1,2,3))

# 重新创建调查设计并运行svyby
df_svy_z2 <- svydesign(id=~ind, data = df_sogi, weights = ~pw)
svyby(~county, ~sexual_orientation, df_svy_z2, svymean)

方案2：直接计算行百分比与置信区间 #

若需直接得到行百分比，可结合svytable和svyciprop实现：

# 生成加权列联表
svytab <- svytable(~county + sexual_orientation, df_svy_z2)
# 计算行百分比（按县分组）
row_percent <- prop.table(svytab, margin = 1) * 100

# 批量计算每个单元格的置信区间
library(purrr)
library(tidyr)

# 生成所有县和性取向的组合
all_combinations <- expand.grid(
  county = levels(df_sogi$county),
  sexual_orientation = levels(df_sogi$sexual_orientation),
  stringsAsFactors = FALSE
)

# 计算每个组合的置信区间
ci_results <- all_combinations %>%
  mutate(
    ci = pmap(list(county, sexual_orientation), function(c, so) {
      # 构造子集设计
      subset_design <- subset(df_svy_z2, county == c)
      # 若子集无数据，返回NA
      if (nrow(subset_design$variables) == 0) {
        return(tibble(lower = NA, upper = NA))
      }
      # 计算比例的置信区间（logit方法更适合极端比例）
      prop_est <- svyciprop(~I(sexual_orientation == so), subset_design, method = "logit")
      tibble(
        lower = confint(prop_est)[1] * 100,
        upper = confint(prop_est)[2] * 100
      )
    })
  ) %>%
  unnest(ci)

# 合并百分比与置信区间结果
final_results <- cbind(as.data.frame(row_percent), ci_results[,3:4])
colnames(final_results) <- c("row_percent", "ci_lower", "ci_upper")

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