R语言survey包进行调查数据分析

5

survey包是R语言中用于复杂调查数据分析的重要工具,它专门设计用于处理具有分层、聚类或多阶段抽样设计的调查数据,能够正确计算标准误和进行统计推断。

survey包主要用于:

  • 处理复杂调查设计(分层、聚类、权重)
  • 计算考虑设计效应的统计量
  • 进行回归分析和假设检验
  • 估计总体参数和标准误

核心功能

  • 调查设计对象创建:svydesign(), svrepdesign()
  • 描述性统计:svymean(), svytotal(), svyquantile()
  • 回归分析:svyglm(), svyolr()
  • 假设检验:svyttest(), svychisq()
  • 可视化:svyboxplot(), svyplot()

示例数据集(API,he Academic Performance Index)

library(survey)
data(api)
dplyr::glimpse(apipop)
## Rows: 6,194
## Columns: 37
## $ cds      <chr> "01611190130229", "01611190132878", "01611196000004", "016111…
## $ stype    <fct> H, H, M, E, E, E, E, E, M, E, E, M, E, E, E, E, H, E, M, E, H…
## $ name     <chr> "Alameda High", "Encinal High", "Chipman Middle", "Lum (Donal…
## $ sname    <chr> "Alameda High", "Encinal High", "Chipman Middle", "Lum (Donal…
## $ snum     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ dname    <chr> "Alameda City Unified", "Alameda City Unified", "Alameda City…
## $ dnum     <int> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 6…
## $ cname    <chr> "Alameda", "Alameda", "Alameda", "Alameda", "Alameda", "Alame…
## $ cnum     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ flag     <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ pcttest  <int> 96, 99, 99, 99, 99, 100, 100, 99, 99, 99, 98, 99, 99, 98, 100…
## $ api00    <int> 731, 622, 622, 774, 811, 780, 808, 739, 795, 650, 668, 701, 7…
## $ api99    <int> 693, 589, 572, 732, 784, 725, 765, 667, 792, 580, 621, 651, 7…
## $ target   <int> 5, 11, 11, 3, 1, 4, 2, 7, 1, 11, 9, 7, 4, 11, NA, NA, 1, NA, …
## $ growth   <int> 38, 33, 50, 42, 27, 55, 43, 72, 3, 70, 47, 50, 70, 79, 47, 3,…
## $ sch.wide <fct> Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Y…
## $ comp.imp <fct> Yes, No, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Ye…
## $ both     <fct> Yes, No, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Ye…
## $ awards   <fct> Yes, No, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Ye…
## $ meals    <int> 14, 20, 55, 35, 15, 25, 22, 50, 10, 71, 48, 32, 23, 76, 8, 3,…
## $ ell      <int> 16, 18, 25, 26, 9, 18, 9, 35, 10, 41, 21, 22, 14, 26, 13, 14,…
## $ yr.rnd   <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ mobility <int> 9, 13, 20, 21, 11, 12, 8, 18, 11, 17, 15, 14, 5, 21, 11, 9, 8…
## $ acs.k3   <int> NA, NA, NA, 20, 20, 20, 20, 20, NA, 20, 20, NA, 21, 19, 19, 1…
## $ acs.46   <int> NA, NA, 26, 30, 29, 29, 26, 31, 29, 25, 29, 30, 28, 26, 28, 2…
## $ acs.core <int> 25, 27, 27, NA, NA, NA, NA, NA, 29, NA, NA, 28, NA, NA, NA, N…
## $ pct.resp <int> 91, 84, 86, 96, 96, 87, 90, 82, 92, 91, 91, 89, 91, 92, 93, 9…
## $ not.hsg  <int> 6, 11, 11, 3, 3, 6, 4, 11, 2, 16, 5, 2, 4, 13, 2, 1, 4, 2, 1,…
## $ hsg      <int> 16, 20, 31, 22, 9, 11, 27, 35, 12, 30, 32, 22, 14, 31, 8, 4, …
## $ some.col <int> 22, 29, 30, 29, 29, 28, 25, 27, 24, 31, 34, 30, 22, 38, 14, 9…
## $ col.grad <int> 38, 31, 20, 31, 26, 41, 32, 21, 34, 20, 22, 32, 37, 13, 43, 4…
## $ grad.sch <int> 18, 9, 8, 15, 33, 13, 12, 6, 28, 3, 7, 14, 24, 4, 33, 40, 36,…
## $ avg.ed   <dbl> 3.45, 3.06, 2.82, 3.32, 3.76, 3.44, 3.23, 2.75, 3.74, 2.66, 2…
## $ full     <int> 85, 90, 80, 96, 95, 90, 93, 93, 89, 100, 96, 97, 91, 100, 100…
## $ emer     <int> 16, 10, 12, 4, 5, 5, 0, 11, 16, 0, 0, 3, 9, 0, 0, 0, 9, 0, 6,…
## $ enroll   <int> 1278, 1113, 546, 330, 233, 276, 180, 363, 850, 263, 369, 792,…
## $ api.stu  <int> 1090, 840, 472, 272, 216, 247, 167, 292, 782, 219, 330, 646, …

加州学校学生表现数据集(API,he Academic Performance Index):
API基于加州学校学生的标准化考试成绩计算。数据涵盖所有学生数不少于100人的学校及其概率抽样样本。

数据结构:
完整人口数据(apipop)包含6194条记录,37个变量,主要包括:

  • cds:唯一标识码
  • stype:学校类型(小学/中学/高中)
  • name/sname:学校名称
  • snum:学校编号
  • dname/dnum:学区名称及编号
  • cname/cnum:县名称及编号
  • pcttest:测试学生比例
  • api00/api99:2000年和1999年API分数
  • target:API变化目标
  • growth:API变化值
  • sch.wide、comp.imp、both:是否达成各类增长目标
  • awards:是否符合奖项资格
  • meals:符合补助餐学生比例
  • ell:英语学习者比例
  • yr.rnd:全年制学校
  • mobility:首次入学学生比例
  • acs.k3、acs.46:K-3年级与4-6年级平均班级人数
  • acs.core:核心课程数
  • pct.resp:已知家长教育水平比例
  • not.hsg、hsg、some.col、col.grad、grad.sch:家长不同教育水平比例
  • avg.ed:家长平均教育水平
  • full、emer:合格教师及紧急资格教师比例
  • enroll:注册学生数
  • api.stu:测试学生数

其他数据集包含抽样权重(pw)及有限总体修正系数(fpc)。

补充说明:

  • apipop:完整数据集
  • apisrs:简单随机样本
  • apiclus1:学区聚类样本(权重有误,应为757/15)
  • apistrat:分层抽样(按学校类型)
  • apiclus2:两阶段聚类样本(学区内学校)

创建调查设计对象

survey包的核心是定义调查设计(如分层、聚类、权重等)。常用函数为svydesign()

  • svydesign()函数说明
svydesign(ids, probs = NULL, strata = NULL, variables = NULL, fpc = NULL,
          data = NULL, nest = FALSE, check.strata = !nest, weights = NULL,
          pps = FALSE, ...)

重要参数

  • ids 指定聚类 ID 的公式或数据框,按从大到小的层级排列。~0~1表示无聚类。

  • probs 指定聚类抽样概率的公式或数据框。

  • strata 指定分层的公式或向量,NULL 表示无分层。

  • variables 指定调查变量的公式或数据框。若为 NULL,则使用 data 参数。

  • fpc 有限总体校正(Finite Population Correction),详见下文。

  • weights 作为 probs 替代的抽样权重(公式或向量)。

  • data 用于查找公式变量的数据框、数据库表名或 imputationList 对象。

  • nest 若为 TRUE,重新标记聚类 ID 以确保其在层内嵌套。

  • pps

    • "brewer":使用 Brewer 近似法进行无放回 PPS 抽样。
    • "overton":使用 Overton 近似法。
    • HR 类对象:使用 Hartley-Rao 近似法。
    • ppsmat 类对象:使用 Horvitz-Thompson 估计量。
  • calibrate.formula 指定权重已校准的模型公式。

  • variance 对于无放回 PPS 抽样,使用 variance = "YG" 表示 Yates-Grundy 估计量,而非默认的 Horvitz-Thompson 估计量。

示例:使用apiclus1数据集(加州学校数据)创建聚类抽样设计:

data(api)  # 加载内置数据集
dclus1 <- svydesign(
  id = ~dnum,          # 聚类变量(学区编号)
  weights = ~pw,       # 抽样权重
  data = apiclus1,     # 数据集
  fpc = ~fpc           # 有限总体校正(FPC)
)

常用分析

1. 计算均值与标准误

svymean(~api00, dclus1)  # 计算变量api00的均值及标准误
##         mean     SE
## api00 644.17 23.542

mean:变量api00的加权均值(644.17),反映总体平均水平。 SE:均值的标准误(23.542),衡量估计的精确性。值越小,估计越可靠。

2. 分位数计算

svyquantile(~api00, dclus1, quantile = c(0.25, 0.5, 0.75), ci = TRUE)
## $api00
##      quantile ci.2.5 ci.97.5       se
## 0.25      552    492     627 31.47166
## 0.5       652    561     714 35.66788
## 0.75      719    696     777 18.88300
## 
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"

quantile:25%、50%(中位数)、75%分位数的估计值。 ci.2.5和ci.97.5:分位数的95%置信区间下限和上限。 se:分位数的标准误,用于评估估计的变异性。

3. 总数估计

svytotal(~stype, dclus1)  # 按学校类型(stype)计算总数
##          total      SE
## stypeE 4873.97 1333.32
## stypeH  473.86  158.70
## stypeM  846.17  167.55

total:各类别(如学校类型stypeE、stypeH)的加权总数,反映总体规模。 SE:总数的标准误,用于评估估计的可靠性。

4. 比率估计

svyratio(~api.stu, ~enroll, dclus1)  # 计算api.stu与enroll的比率
## Ratio estimator: svyratio.survey.design2(~api.stu, ~enroll, dclus1)
## Ratios=
##            enroll
## api.stu 0.8497087
## SEs=
##              enroll
## api.stu 0.008386297

Ratios:分子(api.stu)与分母(enroll)的加权比率(0.85),表示平均每个学生参与测试的比例。 SEs:比率的标准误,衡量比率估计的精度

5. 子集分析

使用subset()对特定子群进行分析:

svyratio(~api.stu, ~enroll, design = subset(dclus1, stype == "H"))  # 仅高中(stype=H)
## Ratio estimator: svyratio.survey.design2(~api.stu, ~enroll, design = subset(dclus1, 
##     stype == "H"))
## Ratios=
##            enroll
## api.stu 0.8300683
## SEs=
##             enroll
## api.stu 0.01472607

结果仅针对高中(stype==“H”)子群,输出形式同比率估计,但标准误已根据子群设计调整,避免低估。

预校准权重(pre-calibrated weights)

预校准权重(pre-calibrated weights) 是指在抽样调查或数据分析中,预先调整样本权重(sampling weights)以使其与已知的总体特征(如人口结构、地理分布等)相匹配,从而提高估计的准确性。

  1. 权重(Weights)
    • 在抽样调查中,不同样本可能代表不同数量的总体个体(如某些群体被过采样或欠采样)。
    • 权重用于调整样本数据,使其能更准确地反映总体特征。
  2. 校准(Calibration)
    • 通过统计方法(如回归、迭代比例拟合等)调整权重,使样本的加权总和与已知的总体总量(如人口总数、性别比例等)一致。
  3. 预校准(Pre-calibration)
    • 权重在数据发布前已由调查机构(如官方统计部门)完成校准,可直接使用,无需再调整。

为什么使用预校准权重?

  • 减少偏差:确保样本估计(如均值、比例)与总体真实值更接近。
  • 提高效率:避免用户自行校准的复杂计算。
  • 标准化:不同研究使用相同的校准基准,结果可比性更强。

1. 基础设计对象(未校准)

library(survey)
data(api)

# 创建基础调查设计对象(未使用校准信息)
dclus1 <- svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc)

# 检查权重总和(已知总体总数应为6194)
sum(weights(dclus1))  # 输出:6194(但标准误差非零)
## [1] 6194
# 添加辅助变量"one"(全1列,用于总量估计)
dclus1 <- update(dclus1, one = rep(1, nrow(apiclus1)))

# 估计总体总数(标准误差错误,应为0)
svytotal(~one, dclus1)  # 输出:total=6194, SE=1442.9)
##     total     SE
## one  6194 1442.9
  • id = ~dnum 指定聚类变量(学区编号)。
  • weights = ~pw 使用数据集中的预校准权重 pw
  • fpc = ~fpc 包含有限总体校正(finite population correction)。
  • 由于未告知 R 权重已校准,svytotal 计算的标准误差不正确(应为0)。

2. 显式校准(手动校准)

# 手动校准到已知总体总量(6194)
cal_dclus1 <- calibrate(dclus1, formula = ~1, population = sum(weights(dclus1)))

# 重新估计总数(标准误差正确)
svytotal(~one, cal_dclus1)  # 输出:total=6194, SE=0
##     total SE
## one  6194  0
# 检查权重是否变化(应不变)
summary(weights(cal_dclus1) / weights(dclus1))  # 输出:全为1
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1
  • calibrate() 显式将设计对象校准到总体总量(~1 表示仅校准到总数)。
  • population 参数指定已知总量(此处用样本权重总和替代,因已知 sum(pw)=6194)。
  • 权重未改变,但标准误差被修正为0(反映校准特性)。

3. 直接使用 calibrate.formula(推荐)

# 在创建设计对象时指定校准变量(~1表示校准到总数)
precal_dclus1 <- svydesign(
  id = ~dnum, 
  weights = ~pw, 
  data = apiclus1, 
  fpc = ~fpc, 
  calibrate.formula = ~1  # 关键:声明权重已校准到总数
)

# 验证标准误差
precal_dclus1 <- update(precal_dclus1, one = rep(1, nrow(apiclus1)))
svytotal(~one, precal_dclus1)  # 输出:total=6194, SE=0
##     total SE
## one  6194  0
  • calibrate.formula = ~1 直接在设计中声明权重已校准到总体总数。
  • 效果与手动校准相同,但更简洁,避免后续重复调用 calibrate()

4. 校准对均值与总量的影响

# 未校准设计(错误的标准误差比例)
enroll_t <- svytotal(~enroll, dclus1)  # SE=932235
enroll_m <- svymean(~enroll, dclus1)   # SE=45.191
SE(enroll_t) / 6194  # 输出:150.506 ≠ SE(enroll_m)
##          enroll
## enroll 150.5061
# 校准后设计(标准误差比例正确)
cenroll_t <- svytotal(~enroll, precal_dclus1)  # SE=279915
cenroll_m <- svymean(~enroll, precal_dclus1)   # SE=45.191
SE(cenroll_t) / 6194  # 输出:45.191 == SE(cenroll_m)
##          enroll
## enroll 45.19137
  • 校准后,总量的标准误差与均值的标准误差满足 SE(total)/N = SE(mean)
  • 这是校准权重的数学特性,确保推断的一致性。

注意事项:

  1. 仅用于已校准权重:若权重未校准,强制使用 calibrate.formula会导致标准误差错误。
  2. 子集数据需重新校准
subset_data <- subset(apiclus1, sch.wide == "Yes")
# 错误:子集权重未重新校准!
subset_design <- svydesign(id = ~dnum, weights = ~pw, data = subset_data, 
                            calibrate.formula = ~1)
# 正确做法:对子集重新校准或使用原始设计对象。

子总体估计

在调查数据分析中估计子总体(domain)的均值或总和,并计算其标准误差。

1. 直接使用survey包的子集功能

这是最简单的方法,使用subset()函数创建子总体的调查设计对象,然后使用svymean()计算均值。

library(survey)
data(fpc)
dfpc <- svydesign(id=~psuid, strat=~stratid, weight=~weight, data=fpc, nest=TRUE)
dsub <- subset(dfpc, x>4)
svymean(~x, design=dsub)
##    mean     SE
## x 6.195 0.7555

均值估计:子总体(x>4)的x均值为 6.195。
标准误:均值的标准误差为 0.7555,反映估计的精确度。
统计意义:若假设真实均值为0,t = 6.195/0.7555 ≈8.2,远大于临界值(通常|t|>2为显著),说明估计值显著不为零。

2. 回归方法

通过回归模型(无截距项)估计子总体均值,二元协变量的系数即为各子总体的均值估计。

summary(svyglm(x ~ I(x>4) + 0, design=dfpc))
## 
## Call:
## svyglm(formula = x ~ I(x > 4) + 0, design = dfpc)
## 
## Survey design:
## svydesign(id = ~psuid, strat = ~stratid, weight = ~weight, data = fpc, 
##     nest = TRUE)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## I(x > 4)FALSE   3.3143     0.3117   10.63 0.000127 ***
## I(x > 4)TRUE    6.1950     0.7555    8.20 0.000439 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.557379)
## 
## Number of Fisher Scoring iterations: 2

回归系数:

I(x>4)FALSE:x≤4的子总体均值为 3.3143(SE=0.3117)。 I(x>4)TRUE:x>4的子总体均值为 6.195(SE=0.7555),与直接子集法结果一致。 显著性:两组的p值均远小于0.001(***),说明均值显著差异。

3. 比率估计方法

将子总体均值视为一个比率估计问题,分子为子总体内的变量值,分母为指示是否属于子总体的指标。

svyratio(~I(x*(x>4)), ~as.numeric(x>4), dfpc)
## Ratio estimator: svyratio.survey.design2(~I(x * (x > 4)), ~as.numeric(x > 4), 
##     dfpc)
## Ratios=
##                as.numeric(x > 4)
## I(x * (x > 4))             6.195
## SEs=
##                as.numeric(x > 4)
## I(x * (x > 4))         0.7555129

比率估计: 将x>4的均值视为分子(x * I(x>4))与分母(I(x>4))的比率,结果同样为 6.195(SE=0.7555)。

三种方法(子集、回归、比率)的估计值和标准误完全一致,说明方法的正确。

复杂调查设计分析示例

1. 校准(GREG)估计

data(api)
dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
pop.totals <- c('(Intercept)'=6194, stypeH=755, stypeM=1018)
gclus1 <- calibrate(dclus1, ~stype+api99, c(pop.totals, api99=3914069))
svymean(~api00, subset(gclus1, comp.imp=="Yes"))
##         mean     SE
## api00 672.05 6.5182

均值估计:完成计算机项目(comp.imp=“Yes”)的学校API均值为 672.05。
标准误:6.5182 较小,说明估计精度较高(可能因样本量大或设计效应低)。
实际意义:该均值显著高于未完成项目的学校(comp.impNo的均值=649.706,SE=12.563),可能暗示计算机项目对学术表现有正向影响。

2. 两阶段抽样

data(mu284)
dmu284 <- svydesign(id=~id1+id2, fpc=~n1+n2, data=mu284)
svymean(~y1, subset(dmu284, y1>40))
##     mean     SE
## y1 48.69 1.1088

均值估计:y1>40的子总体均值为 48.69(SE=1.1088)。
设计影响:两阶段抽样通过fpc参数校正了有限总体效应,标准误更准确。

3. 分层两阶段抽样

library("survival")
data(nwtco)
set.seed(100)
nwtco$incc2<-as.logical(with(nwtco, ifelse(rel | instit==2,1,rbinom(nrow(nwtco),1,.1))))
dccs8<-twophase(id=list(~seqno,~seqno), strata=list(NULL,~interaction(rel,stage,instit)), data=nwtco, subset=~incc2)
svymean(~rel, subset(dccs8,age>36))
##       mean     SE
## rel 0.1845 0.0115
svyratio(~I(rel*as.numeric(age>36)), ~as.numeric(age>36), dccs8)
## Ratio estimator: svyratio.twophase2(~I(rel * as.numeric(age > 36)), ~as.numeric(age > 
##     36), dccs8)
## Ratios=
##                               as.numeric(age > 36)
## I(rel * as.numeric(age > 36))            0.1845032
## SEs=
##                               as.numeric(age > 36)
## I(rel * as.numeric(age > 36))           0.01151891
summary(svyglm(rel~I(age>36)+0, dccs8))
## 
## Call:
## svyglm(formula = rel ~ I(age > 36) + 0, design = dccs8)
## 
## Survey design:
## twophase2(id = id, strata = strata, probs = probs, fpc = fpc, 
##     subset = subset, data = data, pps = pps)
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## I(age > 36)FALSE 0.102666   0.007609   13.49   <2e-16 ***
## I(age > 36)TRUE  0.184503   0.011519   16.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1200944)
## 
## Number of Fisher Scoring iterations: 2

概率估计:诊断年龄>36个月的儿童复发概率为 18.45%(SE=1.15%)。
显著性:与低年龄组(I(age>36)FALSE的均值=10.27%,SE=0.76%)相比,高龄组复发风险显著更高(t=13.49, p<2e-16)。
临床意义:年龄可能是肿瘤复发的风险因素,需进一步研究。

关键点

  1. 子总体估计不能简单地对数据子集进行分析,否则会低估标准误差
  2. 三种方法(子集法、回归法、比率法)在数学上是等价的
  3. survey包中的subset()方法会自动处理权重调整
  4. 这些方法适用于各种复杂调查设计,包括分层、聚类、多阶段抽样和校准估计