使用svyby()处理调查加权数据时确保比例置信下限非负

你当前使用survey包计算加权分组比例的置信区间时，出现了下限为负的情况，原代码及输出如下：

原代码 #

svy_design <- svydesign(ids = ~block, strata = ~strata_study, weights = ~wt, data = data, nest = TRUE, fpc = NULL)

svyby(~binary_outcome,by=~binary_exposure, design=svy_design, FUN=svymean, data=data,keep.var = TRUE, vartype = "ci")

运行结果 #

binary_exposure binary_outcome0 binary_outcome1 ci_l.binary_outcome0 ci_l.binary_outcome1 ci_u.binary_outcome0 ci_u.binary_outcome1
1               1        0.9858986       0.01410138             0.9782720           0.006474771             0.9935252            0.02172798
2               2        0.9264685       0.07353151             0.8489539          -0.003983028             1.0039830            0.15104605

出现负下限的原因是svymean默认使用正态近似计算置信区间，当比例接近0或1时，正态分布的假设会导致区间超出0-1的合理范围。以下是三种解决方法：

方法1：手动截断置信区间 #

直接将超出0-1范围的置信区间边界替换为0或1，适合快速处理结果：

# 假设结果存储在result变量中
result <- svyby(~binary_outcome,by=~binary_exposure, design=svy_design, FUN=svymean, data=data,keep.var = TRUE, vartype = "ci")

# 把所有下限小于0的替换为0，上限大于1的替换为1
result[, grep("ci_l", colnames(result))] <- pmax(result[, grep("ci_l", colnames(result))], 0)
result[, grep("ci_u", colnames(result))] <- pmin(result[, grep("ci_u", colnames(result))], 1)

方法2：使用变换后的置信区间 #

通过对数或logit变换将比例转换为无界尺度，计算CI后再反变换回原尺度，避免超出0-1范围：

针对接近0的比例：对数变换 #

# 自定义计算对数变换置信区间的函数
log_ci_fun <- function(x, design) {
  mean_est <- svymean(x, design, keep.var = TRUE)
  p <- coef(mean_est)
  se <- sqrt(vcov(mean_est))
  
  # 对数变换计算CI
  log_p <- log(p)
  log_se <- se / p
  log_ci_l <- log_p - qnorm(0.975) * log_se
  log_ci_u <- log_p + qnorm(0.975) * log_se
  
  # 反变换回概率尺度
  ci_l <- exp(log_ci_l)
  ci_u <- exp(log_ci_u)
  
  return(c(mean = p, ci_l = ci_l, ci_u = ci_u))
}

# 分组计算
svyby(~binary_outcome1, by=~binary_exposure, design=svy_design, FUN=log_ci_fun)

通用0-1比例：logit变换 #

# 自定义logit变换置信区间函数
logit_ci_fun <- function(x, design) {
  mean_est <- svymean(x, design, keep.var = TRUE)
  p <- coef(mean_est)
  se <- sqrt(vcov(mean_est))
  
  # logit变换：log(p/(1-p))
  logit_p <- log(p / (1 - p))
  # logit尺度的标准误
  logit_se <- se / (p * (1 - p))
  
  logit_ci_l <- logit_p - qnorm(0.975) * logit_se
  logit_ci_u <- logit_p + qnorm(0.975) * logit_se
  
  # 反变换回概率
  ci_l <- exp(logit_ci_l) / (1 + exp(logit_ci_l))
  ci_u <- exp(logit_ci_u) / (1 + exp(logit_ci_u))
  
  return(c(mean = p, ci_l = ci_l, ci_u = ci_u))
}

# 分组计算
svyby(~binary_outcome1, by=~binary_exposure, design=svy_design, FUN=logit_ci_fun)

方法3：使用svyciprop专门计算比例置信区间 #

survey包提供的svyciprop函数针对比例优化了置信区间计算，默认使用logit方法，自动确保CI在0-1范围内：

# 自定义函数调用svyciprop
prop_ci_fun <- function(x, design) {
  ci_result <- svyciprop(x, design, method = "logit")
  return(c(mean = coef(ci_result), ci_l = confint(ci_result)[1], ci_u = confint(ci_result)[2]))
}

# 分组计算
svyby(~binary_outcome1, by=~binary_exposure, design=svy_design, FUN=prop_ci_fun)

如果你需要同时计算binary_outcome0和binary_outcome1的CI，可以分别调用函数后合并结果。

内容的提问来源于stack exchange，提问作者s.stats