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()
方法会自动处理权重调整- 这些方法适用于各种复杂调查设计,包括分层、聚类、多阶段抽样和校准估计