R语言survey包进行调查数据分析
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）以使其与已知的总体特征（如人口结构、地理分布等）相匹配，从而提高估计的准确性。
- 权重（Weights）
- 在抽样调查中，不同样本可能代表不同数量的总体个体（如某些群体被过采样或欠采样）。
- 权重用于调整样本数据，使其能更准确地反映总体特征。
- 校准（Calibration）
- 通过统计方法（如回归、迭代比例拟合等）调整权重，使样本的加权总和与已知的总体总量（如人口总数、性别比例等）一致。
- 预校准（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)。 - 这是校准权重的数学特性，确保推断的一致性。
注意事项:
- 仅用于已校准权重：若权重未校准，强制使用
calibrate.formula会导致标准误差错误。 - 子集数据需重新校准：
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）。
临床意义：年龄可能是肿瘤复发的风险因素，需进一步研究。
关键点
- 子总体估计不能简单地对数据子集进行分析，否则会低估标准误差
- 三种方法（子集法、回归法、比率法）在数学上是等价的
survey包中的subset()方法会自动处理权重调整- 这些方法适用于各种复杂调查设计，包括分层、聚类、多阶段抽样和校准估计