在R语言中计算加权均值置信区间的方法正确性及工具咨询

问题背景 #

我正在尝试用R语言计算加权均值的置信区间，示例数据如下：

t <- tibble(value = c(3, 6, 2, 8, 4), weights = c(0.9, 0.8, 1, 1.2, 1.1))
t
# A tibble: 5 × 2
  value weights
  <dbl>   <dbl>
1     3     0.9
2     6     0.8
3     2     1  
4     8     1.2
5     4     1.1

我用stats::weighted.mean(value)估算加权均值，用Hmisc::wtd.var(t$value)估算方差，再通过sqrt(Hmisc::wtd.var(t$value))得到标准差，标准误则是sqrt(Hmisc::wtd.var(t$value))/sqrt(nrow(t))。

请问我下面这种估算置信区间的方式是否正确？有没有现成的函数可以直接完成这个操作？

mean = stats::weighted.mean(value)
sd = sqrt(Hmisc::wtd.var(t$value))
n = nrow(t),
se = sd / sqrt(n)

ci_lower = mean - qt(0.975, df = n - 1) * se,
ci_upper = mean + qt(0.975, df = n - 1) * se

问题分析与正确做法 #

现有代码的问题 #

1. 加权方差参数缺失：Hmisc::wtd.var必须传入权重参数，仅传t$value会报错，正确写法是Hmisc::wtd.var(t$value, t$weights)。
2. 标准误计算错误：加权均值的标准误不能直接用加权标准差除以sqrt(n)——权重会改变有效样本量，正确公式为sqrt(wtd_var / n_eff)，其中n_eff是有效样本量。
3. 自由度选择不当：不能直接用原始样本量n-1作为自由度，应基于有效样本量计算，有效样本量公式为sum(wts)^2 / sum(wts^2)。

手动正确计算代码 #

library(tibble)
library(Hmisc)

t <- tibble(value = c(3, 6, 2, 8, 4), weights = c(0.9, 0.8, 1, 1.2, 1.1))

# 加权均值
wtd_mean <- weighted.mean(t$value, t$weights)
# 加权方差（传入权重）
wtd_var <- wtd.var(t$value, t$weights)
# 有效样本量
n_eff <- sum(t$weights)^2 / sum(t$weights^2)
# 标准误
se <- sqrt(wtd_var / n_eff)
# 95%置信区间（用有效样本量减1做自由度）
ci_lower <- wtd_mean - qt(0.975, df = n_eff - 1) * se
ci_upper <- wtd_mean + qt(0.975, df = n_eff - 1) * se

# 输出结果
data.frame(
  加权均值 = wtd_mean,
  标准误 = se,
  置信区间下限 = ci_lower,
  置信区间上限 = ci_upper,
  有效样本量 = n_eff
)

现成函数推荐 #

1. survey包（专业加权统计） #

survey包专为加权数据的统计推断设计，自动处理有效样本量和标准误：

library(survey)

# 创建加权调查设计对象
design <- svydesign(ids = ~1, weights = ~weights, data = t)
# 计算加权均值及95%置信区间
svymean(~value, design, conf.level = 0.95)

2. weights包的wtd.t.test函数 #

直接计算加权均值的置信区间，输出包含完整统计量：

library(weights)

wtd.t.test(t$value, weight = t$weights, conf.int = 0.95)

内容的提问来源于stack exchange，提问作者T. C. Nobel