library(survey)
##加权，定义survey数据集
NHANES_all, <- svydesign(id = ~psu,#id处传数据中的psu变量，代表告知指初级抽样单位，不需要处理直接传入
strata = ~strata,#strata处传strata变量，代表告知分层指标，不需要处理直接传入
weights = ~persWeight, #nhanes数据分析中weight需要提前计算
nest = TRUE,
data = nhanesAnalysis #数据集)
#取亚数据集
NHANES <- subset(NHANES_all, #survey数据集
inAnalysis #!duplicated(RXQ_RX_E$SEQN) 筛选条件
)
#计算加权均值及标准差，定义函数
getSummary <- function(varformula, byformula, design){
# Get mean, stderr, and unweighted sample size
c <- svyby(varformula, byformula, design, unwtd.count )
p <- svyby(varformula, byformula, design, svymean )
outSum <- left_join(select(c,-se), p)
outSum
}
#计算
getSummary(~Depression, ~one, NHANES)
#' By sex
getSummary(~Depression, ~Gender, NHANES)
#' By age
getSummary(~Depression, ~Age.Group, NHANES)
#' By sex and age
getSummary(~Depression, ~Gender + Age.Group, NHANES)
#data_2为原数据集，design_new为svydesign以后的新数据集
#计算数据中的年龄及种族加权及未加权的均值或比例，可以看出加权及未加权结果有很大差异，对数据进行基线信息描述时，应该使用加权结果。
#unweighted age and se
mean(data_2$age,na.rm=T)
#49.54916
# weighted age and se
svymean(~age, design_new, na.rm = TRUE)
#45.874
#' Proportion of unweighted interview sample
data_2 %>% count(race) %>%
mutate(prop= round(n / sum(n)*100, digits=1))
#' Proportion of weighted interview sample
data_2 %>% count(race, wt=WTMEC2YR) %>%
#svyglm分析
#使用常规的glm和weighted glm会对结果进行有偏估计，应该在构建survey数据库的基础上，进行svyglm分析，以下是三个方法的比较
#glm
Result2 <- glm(TT4~Bromoform+age+Gender+race+BMI+Education,
family = gaussian(), data=data_2)
summary(Result2)
#weighted glm
Result3 <- glm(TT4~Bromoform+age+Gender+race+BMI+Education,
family = gaussian(), data=data_2,weights =WTMEC2YR )
summary(Result3)
#survey-weighted glm
Result1 <- svyglm(TT4~Bromoform+age+Gender+race+BMI+Education,
family = gaussian(), data=data_2,design=design_new)
summary(Result1)
背景知识
创建样本权重是为了考虑复杂的调查设计（包括过度抽样）、调查无响应和后分层（以确保计算出的估计值代表美国平民非机构化人口NHANES复杂抽样导致需要）计算权重
The sample weights are created to account for the complex survey design (including oversampling), survey nonresponse, and post-stratification in order to ensure that calculated estimates are representative of the U.S. civilian noninstitutionalized population.
权重选择
①所有变量都是以in-home interview的方式收集，权重采用wtint2yr；
只有研究的所有变量都是以in-home interview的方式收集的，权重才用wtint2yr，如果还有其他方式收集的变量，接着往下看；
②部分变量是以MEC的方式收集的，权重采用wtmec2yr；
如果所有变量只有以in-home interview和MEC检查两种方式收集的，权重就用wtmec2yr，如果还有子样本变量（没有24-hour dietary recall变量），参考第③点选择权重；如果包含24-hour dietary recall变量则直接参考第④点选择权重；
③部分变量是调查子样本的一部分，则采用相应子样本权重；
如果研究的变量同时包括wtint2yr，wtmec2yr以及自样本权重这个三种权重，则选择相应子样本权重；
因为所有参与者都接受采访（in-home interview）（人群1），在接受采访的人中部分接受MEC检查（人群2），在接受MEC检查的人中，只有空腹8小时以上的人才检查了空腹甘油三酯（人群3：子样本人群），即人群3（子样本人群）<人群2<人群1，根据权重选择的核心原则，所以选择子样本变量对应的权重为最终权重。
④一些变量来自24小时饮食回忆（24-hour dietary recall）
24小时饮食回忆不属于子样本变量，但是完成这部分调查的参与者其权重比较特殊，一周中工作日和周末的饮食摄入量可能会存在差异，该权重可以调整这些差异。只要研究中有变量是以24-hour dietaryrecall方式收集的（不管是否包括wtint2yr，wtmec2yr或者子样本权重），最终权重都为wtdrd1（第一天）/wtdrd2（2天）。
结合周期计算权重
（合并多个周期时需合并权重，这里只针对1999年后的数据）
先根据以上原则选择相应的权重类型，然后根据合并的周期重新计算，以下以wtint2yr权重为例。
总原则：对于任何不包括1999-2000年的2001-2002年及以后的任何周期合并，权重都等于1/周期数*（相应权重）。
①如果只合并1999-2000和2001-2002四年（2个周期）的数据，则最终权重为wtint4yr（全部周期中，只有合并了1999-2002四年的数据才有wtint4yr这个权重类型，其他所有的都是wtint2yr）；
②如果合并1999-2004六年（3个周期）的数据，需分成1999-2002年（作为整体算2个周期）和2003-2004（1个周期）两部分考虑：
1999-2002年（2个周期）的权重为2/3*wtint4yr；
2003-2004年（1个周期）的权重为1/3*wtint2yr
③如果合并2001-2002和2003-2004四年（2个周期）的数据，则最终权重为1/2*wtint2yr；
④如果合并2001-2006六年（3个周期）的数据，因为没有包括1999-2000年的数据，所以最终权重为1/3*wtint2yr。
至此基本可满足大部分研究的权重选择和计算，计算好之后将最终的权重数据传入svydesign（）方法中weigths即完成了抽样方式的设置，这样就可以进行后续的分析了，包括根据纳入和排除标准选择最终纳入分析的人群以及统计分析。
文献
参考资料
1.NHANES数据库的介绍及使用（二） https://blog.csdn.net/weixin_40563866/article/details/120291633
- NHANES数据库权重计算 https://blog.csdn.net/qq_42458954/article/details/121296965
- NHANES官方文档 NHANES 教程 - 模块 3 - 加权 (cdc.gov)
