使用svyby和svymean计算单水平因子变量的CV与CI问题

使用survey包计算子区域的吸烟率及对应置信区间（CI）、变异系数（CV）时，常规代码如下：

svyby(~factor(SMOKING), ~COUNTY, na.rm=TRUE, dsgn, svymean, vartype=c("ci", "cv"))

但当某个子区域的SMOKING因子仅存在一个取值时（比如示例数据里的COUNTY=9，SMOKING只有取值1），会触发如下报错：

Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels

示例数据（行代表SMOKING取值，列代表COUNTY）：

1    3    5    7    9   11   13   15   17   19   21   23   25   27  777  999  Sum
  1     29   46   45   63    7   54   11   24   37   39   61   49   50   55   19    4  593
  2     18   36   28   69    0   28    8   18   21   24   29   50   27   45    5    5  411
  Sum   47   82   73  132    7   82   19   42   58   63   90   99   77  100   24    9 1004

解决方案：自定义包装函数处理单水平场景 #

无需提前过滤单水平子区域，可通过自定义包装函数自动适配两种场景：

library(survey)

# 自定义svymean包装函数，支持单水平因子
svymean_single_level <- function(formula, design, ...) {
  # 提取目标变量
  target_var <- all.vars(formula)[1]
  subset_data <- model.frame(design)[[target_var]]
  
  # 获取因子的所有原始水平（确保结果维度一致）
  all_levels <- levels(factor(subset_data))
  current_levels <- unique(subset_data)
  
  if (length(current_levels) == 1) {
    # 单水平场景：手动构造结果
    current_level <- as.character(current_levels[1])
    
    # 均值向量：对应水平为1（100%），其余为0（0%）
    mean_vals <- rep(0, length(all_levels))
    names(mean_vals) <- all_levels
    mean_vals[current_level] <- 1
    
    # 置信区间：方差为0，上下限与均值一致
    ci_vals <- matrix(rep(mean_vals, 2), nrow = length(all_levels), ncol = 2)
    colnames(ci_vals) <- c("ci_lower", "ci_upper")
    rownames(ci_vals) <- all_levels
    
    # 变异系数：方差为0时设为0（也可根据需求改为NA）
    cv_vals <- rep(0, length(all_levels))
    names(cv_vals) <- all_levels
    
    # 整理为svymean返回的格式
    result <- list(
      mean = mean_vals,
      var = diag(rep(0, length(all_levels))),
      ci = ci_vals,
      cv = cv_vals
    )
    class(result) <- "svystat"
  } else {
    # 多水平场景：正常调用svymean
    result <- svymean(formula, design, ...)
  }
  return(result)
}

# 使用自定义函数执行分组计算
svyby(~factor(SMOKING), ~COUNTY, na.rm=TRUE, dsgn, svymean_single_level, vartype=c("ci", "cv"))

函数逻辑说明 #

• 单水平场景：针对仅有一种取值的子区域，直接生成100%/0%的均值，置信区间与均值完全重合（因方差为0），变异系数设为0；
• 多水平场景：保持svymean的原有逻辑，正常计算均值、CI和CV。

这样就能让代码自动处理所有子区域，无需提前筛选数据。

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