NHANES数据(复杂调查数据)亚组交互函数2.2版（P for interaction）发布---用于一键生成交互效应表

写在前面的话，本函数只支持NHANES数据(复杂调查数据)的逻辑回归和线性回归，其他类型均不支持，请注意甄别，电子产品，买错不能退换。

在SCI文章中，交互效应表格（通常是表五）能为文章锦上添花，增加文章的信服力，增加结果的可信程度，还能进行数据挖掘。什么是亚组，通常就是特殊类型人群，比如男女，种族等，就是说你的数据放入特殊人群中结果还可靠吗？如果在各个特殊人群中，你的结果很稳定，说明你的结论很可靠。如果亚组的结论和你的数据数据结论相反，你可以拿来做个新论题。还可以比较不同亚组之间有无区别，比如做了心脏支架和没做支架的区别，可以发现很多新思路，易于数据挖掘。

在既往文章《代码+视频–NHANES数据(复杂调查数据)亚组交互函数1.6尝鲜版（P for interaction）发布—用于一键生成交互效应表》中，咱们发布了svy.scitb5函数，反响还不错，基本没啥大问题，有个小问题是下载到的NHANES数据需要整理一下，不能直接拿来分析，比如像这种带连接符的数据

不单是svy.scitb5识别不出来，survey包也是识别不出来，所以咱们一定不能偷懒，得到的nhanes数据一定要清洗一下，把它改成汉字或字符改成数字。

本次发布的2.2函数版本，构架变动挺大的，除了修正了前面版本的一些小问题，之前1.6版本的结果和手动验证还是有点差别，新版本改进了算法，并且增加了线性回归的森林图绘制，下面我来演示一下。先导入R包和数据

library(survey)
bc<-read.csv("E:/r/test/subtext.csv",sep=',',header=TRUE)
bc <- na.omit(bc)
str(bc)

我介绍一下数据，SEQN：序列号，RIAGENDR, # 性别， RIDAGEYR, # 年龄，RIDRETH1, # 种族，DMDMARTL, # 婚姻状况，WTINT2YR,WTMEC2YR, # 权重，SDMVPSU, # psu，SDMVSTRA，# strata，LBDGLUSI, #血糖mmol表示，LBDINSI, #胰岛素( pmmol/L)，PHAFSTHR #餐后血糖，LBXGH #糖化血红蛋白，SPXNFEV1, #FEV1：第一秒用力呼气量，SPXNFVC #FVC：用力肺活量，ml（估计肺容量），LBDGLTSI #餐后2小时血糖，factor.FVC是我把肺活量分为了2分类，方便用于测试。

把分类变量转成因子

bc$DMDMARTL<-ifelse(bc$DMDMARTL==1,1,0)
bc$RIAGENDR<-as.factor(bc$RIAGENDR)
bc$RIDRETH1<-as.factor(bc$RIDRETH1)
bc$DMDMARTL<-as.factor(bc$DMDMARTL)
bc$factor.FVC<-as.factor(bc$factor.FVC)

建立抽样调查函数

bcSvy2<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
                   nest=TRUE,data = bc)

设定协变量和交互变量，这里要注意一下cov1指协变量，Interaction为交互变量，交互变量，交互变量必须为分类变量

cov1<-c("RIAGENDR","LBDINSI","RIDRETH1")		
Interaction<-c("DMDMARTL","RIDRETH1")

导入函数svyscitb5函数

source("E:/r/test/2.2svyscitb5.R")

成功导入后应该显示13个函数

函数的主体为

svy.scitb5 (data,x,y,Interaction,cov=NULL,time=NULL,family=NULL,method=NULL,
                     ids=NULL,strata=NULL,weights=NULL,svydstr=NULL)

我来解释一下data是数据，必须数据框形式，x是你研究的目标变量，y是你的结局变量，Interaction是你的分层变量，这个必须是分类变量并转成因子，cov是你的协变量，在我的设定中cov是要包含Interaction的，这也符合我们的习惯，但是你也可以使用contain=F把它关掉。统一使用family=“svyglm”，svydstr这里填入调查函数，支持逻辑回归和线性回归，如果你使用过我写的scitb5函数或之前的版本，这个函数使用起来完全没有压力。

先来个X是连续变量Y是分类变量的

out<-svy.scitb5(data=bc,x="RIDAGEYR",y="factor.FVC",Interaction=Interaction,cov = cov1,family="svyglm",
              ids ="SDMVPSU",strata="SDMVSTRA",weights="WTMEC2YR",svydstr=bcSvy2,contain=F) #X连续Y分类

绘制森林图

这个森林图乍一看有点傻，其实修改也是有技巧的，主要是修改它的范围，因为我们看到它的OR其实变化很小，所以范围不要搞太大

plotforest(out,xticks=c(0.9,1,1.1))

上图是不是舒服一点了，所以修改范围也是很重要的。
下面我手动来验证一下对不对，就那婚姻这个分层变量来演示

dat1<-subset(bc,bc$DMDMARTL==0)
dat2<-subset(bc,bc$DMDMARTL==1)
bcSvy21<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
                    nest=TRUE,data = dat1)
bcSvy22<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
                    nest=TRUE,data = dat2)
fit1<-svyglm(factor.FVC ~ RIDAGEYR+RIAGENDR+RIDRETH1, family="quasibinomial",design = bcSvy21)
fit2<-svyglm(factor.FVC ~ RIDAGEYR+RIAGENDR+RIDRETH1,family="quasibinomial", design = bcSvy22)
exp(fit1$coefficients)
exp(fit2$coefficients)

咱们可以看到，手动计算和函数计算完全一样，结果的可靠性是没问题的。
再来个X是分类，Y是分类的

out<-svy.scitb5(data=bc,x="RIAGENDR",y="factor.FVC",Interaction=Interaction,cov = cov1,family="svyglm",
                 svydstr=bcSvy2,dec=5,pdec = 5,p.intervaue = 5,contain=F)  #x分类,Y分类

上图我要解释一下，做X分类变量的时候需要设定一个参考，在DMDMARTL这个分层比较的时候，RIAGENDR.1_DMDMARTL.0是被认定做参考的，什么意思呢？就是当RIAGENDR等于1和DMDMARTL等于0这个亚组的人群是被默认为做参考比较的，其他组都是和它进行比较，分类变量进行亚组交互的时候，分类最好不要太多，要不数据会很大，而且有些层分不到数据就会显示数据缺失NA。
绘制森林图

plotforest(out)

换个Y是连续变量的，这里生成的应该是β，我们可以看到，函数自己对数据类型进行了判断了

out<-svy.scitb5(data=bc,x="LBDINSI",y="SPXNFVC",Interaction=Interaction,cov = cov1,family="svyglm",
                svydstr=bcSvy2)  #x连续,Y连续

旧版本是做不了线性回归森林图的，我看了一下有的文章确实有做线性回归森林图，所以在新版本也开发了这个线性回归森林图，方法也是一样的

plotforest(out)

我们注意一下上图有红箭头，有小数点存在，这和我的算法有关，你可以自己调整一下

plotforest(out,xticks=c(-8,0,2))

如果你还不会也没有关系，可以参考一个旧视频操作，差不多的

代码+视频--NHANES数据(复杂调查数据)亚组交互函数1.4尝鲜版（P for interaction）发布---用于一键生成交互效应表

购买了本函数之前版本的无需重新购买，截图发给我就行

需要获取NHANES数据(复杂调查数据)亚组交互函数2.2版请看下面这篇文章

NHANES数据(复杂调查数据)亚组交互函数2.2版（P for interaction）发布—用于一键生成交互效应表