library(mediation)
library(gmodels)
library(foreign)
library(dplyr)
library(RNHANES)
library(survey)
library(reshape)
library(vctrs)
library(tidyverse)
library(RNHANES)
library(survey)
library(testthat)
library(plyr)
library(epiDisplay)
library(tableone)
library(vctrs)
library(tidyverse)
library(DMwR2)
library(plyr)
library(pROC)
library(car)
library(VIM)
library(mice)
library(Rcpp)
library(pROC)
library(leaps)
library(MASS)
library(tidyverse)
library(caret)
library(stringi)
options(survey.lonely.psu="adjust")
setwd("E:/BaiduSyncdisk")

#################################数据下载
pfas_files <- tribble(
  ~cycle,      ~file_name,
  "2007-2008", "DR1TOT_E",
  "2009-2010", "DR1TOT_F",
  "2011-2012", "DR1TOT_G",
  "2013-2014", "DR1TOT_H",
)
EPH_files <- tribble(
  ~cycle,      ~file_name,

  "2013-2014", "PFAS_H",
  "2013-2014", "SSPFAS_H",
  "2011-2012", "PFC_G",
  "2009-2010", "PFC_F",
  "2007-2008", "PFC_E",
)
pfas2_files <- tribble(             ###我这里没用DR2的数据，也没加载
  ~cycle,      ~file_name,
  "2007-2008", "DR2TOT_E",
  "2005-2006", "DR2TOT_D",
  "2009-2010", "DR2TOT_F",
  "2011-2012", "DR2TOT_G",
  "2013-2014", "DR2TOT_H",
)


bmi_files <- tribble(
  ~cycle,      ~file_name,
  "2005-2006", "BMX_D",
  "2007-2008", "BMX_E",
  "2009-2010", "BMX_F",
  "2011-2012", "BMX_G",
  "2013-2014", "BMX_H",
)

DIQ_files <- tribble(
  ~cycle,      ~file_name,
  "2005-2006", "DIQ_D",
  "2007-2008", "DIQ_E",
  "2009-2010", "DIQ_F",
  "2011-2012", "DIQ_G",
  "2013-2014", "DIQ_H",
  )

ALQ_files <- tribble(
  ~cycle,      ~file_name,
  "2005-2006", "ALQ_D",
  "2007-2008", "ALQ_E",
  "2009-2010", "ALQ_F",
  "2011-2012", "ALQ_G",
  "2013-2014", "ALQ_H",
)

SMQ_files <- tribble(
  ~cycle,      ~file_name,
  "2005-2006", "SMQ_D",
  "2007-2008", "SMQ_E",
  "2009-2010", "SMQ_F",
  "2011-2012", "SMQ_G",
  "2013-2014", "SMQ_H",
)

OCQ_files <- tribble(
  ~cycle,      ~file_name,
  "2005-2006", "OCQ_D",
  "2007-2008", "OCQ_E",
  "2009-2010", "OCQ_F",
  "2011-2012", "OCQ_G",
  "2013-2014", "OCQ_H",
)

MCQ_files <- tribble(
  ~cycle,      ~file_name,
  "2005-2006", "MCQ_D",
  "2007-2008", "MCQ_E",
  "2009-2010", "MCQ_F",
  "2011-2012", "MCQ_G",
  "2013-2014", "MCQ_H",
)
BPX_files <- tribble(
  ~cycle,      ~file_name,
  "2007-2008", "BPX_E",
  "2009-2010", "BPX_F",
  "2011-2012", "BPX_G",
  "2013-2014", "BPX_H",
)
BPQ_files <- tribble(
  ~cycle,      ~file_name,
  "2007-2008", "BPQ_E",
  "2009-2010", "BPQ_F",
  "2011-2012", "BPQ_G",
  "2013-2014", "BPQ_H",
)


BIOPRO_files <- tribble(
  ~cycle,      ~file_name,
  "2007-2008", "BIOPRO_E",
  "2009-2010", "BIOPRO_F",
  "2011-2012", "BIOPRO_G",
  "2013-2014", "BIOPRO_H",
)

COT_files <- tribble(
  ~cycle,      ~file_name,
  "2013-2014", "COT_H",
  "2011-2012", "COTNAL_G",
  "2009-2010", "COTNAL_F",
  "2007-2008", "COTNAL_E",
  "2005-2006", "COT_D",
)
CRP_files <- tribble(
  ~cycle,      ~file_name,
  "2013-2014", "CBC_H",
  "2011-2012", "CBC_G",
  "2009-2010", "CBC_F",
  "2007-2008", "CBC_E",
  "2005-2006", "CBC_D",
)
####################################################数据下载，代码来源于EHP论文，需注意该代码仅可应用于2015年以前的数据批量下载。且需要进行检查。可能会缺损。建议自己下载用它读取。
summary(d$SDDSRVYR)

nhanes_pfas_data <-          nhanes_load_data(file_name = pfas_files$file_name,
                                              year = pfas_files$cycle,
                                              destination = "E:/BaiduSyncdisk/NHANES202110/data",  ###这里是储存数据原始文件的地方
                                              demographics = T)
nhanes_pfas2_data <-          nhanes_load_data(file_name = pfas2_files$file_name,
                                               year = pfas2_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",
                                               demographics = F)
nhanes_EPH_data <-          nhanes_load_data(file_name = EPH_files$file_name,
                                               year = EPH_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",
                                               demographics = F)
nhanes_diq_data<-          nhanes_load_data(file_name = DIQ_files$file_name,
                                            year = DIQ_files$cycle,
                                            destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)


nhanes_BPX_data <-            nhanes_load_data(file_name = BPX_files$file_name,
                                               year = BPX_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data")

nhanes_bmi_data <-            nhanes_load_data(file_name = bmi_files$file_name,
                                               year = bmi_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data")

nhanes_ALQ_data <-            nhanes_load_data(file_name = ALQ_files$file_name,
                                               year = ALQ_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)

nhanes_SMQ_data <-            nhanes_load_data(file_name = SMQ_files$file_name,
                                               year = SMQ_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)

nhanes_OCQ_data <-            nhanes_load_data(file_name = OCQ_files$file_name,
                                               year = OCQ_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)

nhanes_MCQ_data <-            nhanes_load_data(file_name = MCQ_files$file_name,
                                               year = MCQ_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)
nhanes_BPQ_data <-            nhanes_load_data(file_name = BPQ_files$file_name,
                                               year = BPQ_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)
nhanes_BIOPRO_data <-            nhanes_load_data(file_name = BIOPRO_files$file_name,
                                               year = BIOPRO_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)
nhanes_COT_data <-            nhanes_load_data(file_name = COT_files$file_name,
                                                  year = COT_files$cycle,
                                                  destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)
nhanes_CRP_data <-            nhanes_load_data(file_name = CRP_files$file_name,
                                               year = CRP_files$cycle,
                                               destination = "E:/BaiduSyncdisk/NHANES202110/data",demographics = F)
sspfas_reduced2 <- nhanes_EPH_data$SSPFAS_H

sspfas_reduced<-sspfas_reduced2 [,c("SEQN", "SSNPFOA", "SSBPFOA", "SSNPFOS", "SSMPFOS", "SDNPFOAL", "SDBPFOAL", "SDNPFOSL", "SDMPFOSL")]
nhanes_EPH_data$PFAS_H <- left_join(nhanes_EPH_data$PFAS_H, sspfas_reduced, by = "SEQN")
nhanes_EPH_data$PFAS_H$LBXPFOA <- rowSums(nhanes_EPH_data$PFAS_H[,c('SSNPFOA', 'SSBPFOA')])
nhanes_EPH_data$PFAS_H$LBDPFOAL <- ifelse(nhanes_EPH_data$PFAS_H$SDNPFOAL == 0 | nhanes_EPH_data$PFAS_H$SDBPFOAL == 0, 0, 1)
nhanes_EPH_data$PFAS_H$LBXPFOS <- rowSums(nhanes_EPH_data$PFAS_H[,c('SSNPFOS', 'SSMPFOS')])
nhanes_EPH_data$PFAS_H$LBDPFOSL <- ifelse(nhanes_EPH_data$PFAS_H$SDNPFOSL == 0 | nhanes_EPH_data$PFAS_H$SDMPFOSL == 0, 0, 1)
nhanes_EPH_data$SSPFAS_H <- NULL
######################################################################整理数据step1 ，合并各个cycle数据
nhanes_combined_pfas_data <- bind_rows(nhanes_pfas_data) 
nhanes_combined_pfas2_data <- bind_rows(nhanes_pfas2_data) 
nhanes_combined_bmi_data <- bind_rows(nhanes_bmi_data) 
nhanes_combined_SMQ_data <- bind_rows(nhanes_SMQ_data) 
nhanes_combined_ALQ_data <- bind_rows(nhanes_ALQ_data)
nhanes_combined_OCQ_data <- bind_rows(nhanes_OCQ_data) 
nhanes_combined_MCQ_data <- bind_rows(nhanes_MCQ_data)
nhanes_combined_BPQ_data <- bind_rows(nhanes_BPQ_data)
nhanes_combined_BIOPRO_data <- bind_rows(nhanes_BIOPRO_data)
nhanes_combined_EPH_data <- bind_rows(nhanes_EPH_data)
nhanes_combined_COT_data <- bind_rows(nhanes_COT_data)
nhanes_combined_CRP_data <- bind_rows(nhanes_CRP_data)
######################################################STEP2，根据SEQN合并分类数据
nhanes_data <- reduce(list(nhanes_combined_pfas_data,
                           nhanes_combined_pfas2_data,
                           nhanes_combined_ALQ_data,
                           nhanes_combined_SMQ_data,
                           nhanes_combined_OCQ_data,
                           nhanes_combined_MCQ_data,
                           nhanes_combined_BIOPRO_data,
                           nhanes_combined_bmi_data,
                           nhanes_combined_EPH_data,
                           nhanes_combined_BPQ_data,
                           nhanes_combined_COT_data,
                           nhanes_combined_CRP_data
                           ),
                      left_join,
                      by = "SEQN")

##################################################################计算PHEN age
DR1TOT_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DR1TOT_J.xpt")
DR1TOT_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DR1TOT_I.xpt")
DR1TOT<-bind_rows(DR1TOT_J,DR1TOT_I)

DEMO_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DEMO_J.xpt")
DEMO_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DEMO_I.xpt")
DEMO<-bind_rows(DEMO_J,DEMO_I)

DR2TOT_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DR2TOT_J.xpt")
DR2TOT_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DR2TOT_I.xpt")
DR2TOT<-bind_rows(DR2TOT_J,DR2TOT_I)

PFAS_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/PFAS_J.xpt")
PFAS_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/PFAS_I.xpt")
PFAS<-bind_rows(PFAS_J,PFAS_I)

PFAS$LBXPFOA <- rowSums(PFAS[,c('LBXNFOA', 'LBXBFOA')])
PFAS$LBDPFOAL <- ifelse(PFAS$LBDNFOAL == 0 | PFAS$LBDBFOAL == 0, 0, 1)
PFAS$LBXPFOS <- rowSums(PFAS[,c('LBXNFOS', 'LBXMFOS')])
PFAS$LBDPFOSL <- ifelse(PFAS$LBDNFOSL == 0 | PFAS$LBDMFOSL == 0, 0, 1)


BMX_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BMX_J.xpt")
BMX_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BMX_I.xpt")
BMX<-bind_rows(BMX_J,BMX_I)

DIQ_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DIQ_J.xpt")
DIQ_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/DIQ_I.xpt")
DIQ<-bind_rows(DIQ_J,DIQ_I)

ALQ_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/ALQ_J.xpt")
ALQ_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/ALQ_I.xpt")
ALQ<-bind_rows(ALQ_J,ALQ_I)

SMQ_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/SMQ_J.xpt")
SMQ_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/SMQ_I.xpt")
SMQ<-bind_rows(SMQ_J,SMQ_I)

OCQ_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/OCQ_J.xpt")
OCQ_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/OCQ_I.xpt")
OCQ<-bind_rows(OCQ_J,OCQ_I)

MCQ_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/MCQ_J.xpt")
MCQ_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/MCQ_I.xpt")
MCQ<-bind_rows(MCQ_J,MCQ_I)

BPX_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BPX_J.xpt")
BPX_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BPX_I.xpt")
BPX<-bind_rows(BPX_J,BPX_I)

BPQ_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BPQ_J.xpt")
BPQ_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BPQ_I.xpt")
BPQ<-bind_rows(BPQ_J,BPQ_I)

BIOPRO_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BIOPRO_J.xpt")
BIOPRO_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/BIOPRO_I.xpt")
BIOPRO<-bind_rows(BIOPRO_J,BIOPRO_I)

COT_J <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/COT_J.xpt")
COT_I <- read.xport("E:/BaiduSyncdisk/NHANES202110/data/COT_I.xpt")
COT<-bind_rows(COT_J,COT_I)

EPDMO3 <- reduce(list(DR1TOT,
                      DR2TOT,
                      DEMO,
                      PFAS,
                      BMX,
                      DIQ,
                      BIOPRO,
                      SMQ,
                      ALQ,
                      OCQ,
                      MCQ,
                      BPX,
                      COT),
                 left_join,
                 by = "SEQN")


EPDMO2 <- bind_rows(nhanes_data,EPDMO3)







########################################  计算WTS

EPDMO2 <- EPDMO2 %>%
  mutate(WTS = case_when(
    !is.na(.$WTDR2D.x) ~ WTDR2D.x/6
  )) %>%
  filter(!is.na(WTS)) 
EPDMO2<-EPDMO2[EPDMO2$WTS!=0,]
EPDMO2$DIQ[EPDMO2$DIQ010==1]<-2
EPDMO2$DIQ[EPDMO2$DIQ010==2]<-0
EPDMO2$DIQ[EPDMO2$DIQ010==3]<-2
EPDMO2$DIQ1[EPDMO2$DIQ==2]<-1
EPDMO2$DIQ1[EPDMO2$DIQ==0]<-0

####################################################  计算DII
EPDMO2$DR1TN3<-EPDMO2$DR1TP205+EPDMO2$DR1TP225+EPDMO2$DR1TP226
EPDMO2$DR1TN6<-EPDMO2$DR1TP182+EPDMO2$DR1TP183+EPDMO2$DR1TP204
EPDMO2$DR2TN3<-EPDMO2$DR2TP205+EPDMO2$DR2TP225+EPDMO2$DR2TP226
EPDMO2$DR2TN6<-EPDMO2$DR2TP182+EPDMO2$DR2TP183+EPDMO2$DR2TP204
EPDMO2$DR1TCAFF<-EPDMO2$DR1TCAFF/1000
EPDMO2$DR2TCAFF<-EPDMO2$DR2TCAFF/1000
NUT1<-c("DR1TVB12","DR1TVB6","DR1TBCAR","DR1TCARB","DR1TCHOL","DR1TKCAL",
        "DR1TTFAT","DR1TFIBE","DR1TFOLA","DR1TIRON","DR1TMAGN","DR1TMFAT","DR1TNIAC","DR1TPROT","DR1TPFAT",
        "DR1TVB2","DR1TSFAT","DR1TSELE","DR1TVB1","DR1TVARA","DR1TVC","DR1TVD","DR1TATOC","DR1TZINC","DR1TN3","DR1TN6") 
NUT2<-c("DR2TVB12","DR2TVB6","DR2TBCAR","DR2TCARB","DR2TCHOL","DR2TKCAL",
        "DR2TTFAT","DR2TFIBE","DR2TFOLA","DR2TIRON","DR2TMAGN","DR2TMFAT","DR2TNIAC","DR2TPROT","DR2TPFAT",
        "DR2TVB2","DR2TSFAT","DR2TSELE","DR2TVB1","DR2TVARA","DR2TVC","DR2TVD","DR2TATOC","DR2TZINC","DR2TN3","DR2TN6") 
NUT3<-c("TALCO","TVB12","TVB6","TBCAR","TCARB","TCHOL","TKCAL",
        "TTFAT","TFIBE","TFOLA","TIRON","TMAGN","TMFAT","TNIAC","TPROT","TPFAT",
        "TVB2","TSFAT","TSELE","TVB1","TVARA","TVC","TVD","TATOC","TZINC","TN3","TN6") 

DII1 <- function(VAR,DATA,digit = 3){
  varc<-paste0(VAR,"C")
  DIIUSE<-DII[,VAR]
  ZSCOREXXX<-(DATA[,VAR]-DIIUSE[1])/DIIUSE[2]
  P_XXX<-as.data.frame((pnorm(ZSCOREXXX, mean = 0, sd = 1, lower.tail = T)*2-1)*DIIUSE[3])
  names(P_XXX)=varc
  return(P_XXX)
}
DRADD<- function(dat, var,digit = 3){    #线性回归函数
  var1<-paste0("DR1",var)
  var2<-paste0("DR2",var)
  data1<-dat[,var1]
  data2<-dat[,var2]
  data22<-cbind(data1,data2)
  data3<-as.data.frame(rowMeans(data22,na.rm=T))
  names(data3)<-var
  return(data3)
}


DATADII<-lapply(NUT3, function(x) DRADD(var =x,dat=EPDMO2))
newdata <- do.call(cbind, DATADII)
EPDMO2<-cbind(EPDMO2,newdata)
DII<-read.table("E:/BaiduSyncdisk/XUQI/DIIR2.csv",header=T,sep=",")  
des4<-nhanes_survey_design(EPDMO2, "WTS")


A1<-lapply(NUT1, function(x) DII1( VAR = x, DATA = EPDMO2))
A2 <- do.call(cbind, A1)
A3<-as.data.frame(rowSums(A2))
names(A3)<-"DII1"
EPDMO2<-cbind(EPDMO2,A3,by=0)
A1<-lapply(NUT2, function(x) DII1( VAR = x, DATA = EPDMO2))
A2 <- do.call(cbind, A1)
A3<-as.data.frame(rowSums(A2))
names(A3)<-"DII2"

EPDMO2<-cbind(EPDMO2,A3,by=0)
EPDMO2$DII<-rowMeans(EPDMO2[,c("DII1","DII2")],na.rm=T)
summary(EPDMO2$DII1)

des4<-nhanes_survey_design(EPDMO2, "WTS")




######################################################可以初步清理下数据 数据清理模式与PMID: 36152103完全相同，可参考 ###只入选大于等于20岁，一般成人就这样

#######################################种族和教育程度清理，根据常规使用方法分类
EPDMO2$DMDEDUC2[EPDMO2$DMDEDUC2==1]<-2
EPDMO2$DMDEDUC2[EPDMO2$DMDEDUC2==7]<-NA
EPDMO2$DMDEDUC2[EPDMO2$DMDEDUC2==9]<-NA
EPDMO2$RIDRETH1[EPDMO2$RIDRETH1==2]<-1






##############################################################整理饮酒协变量


EPDMO2$ALQ120UU[EPDMO2$ALQ120U==1]<-1
EPDMO2$ALQ120UU[EPDMO2$ALQ120U==2]<-1/30*7
EPDMO2$ALQ120UU[EPDMO2$ALQ120U==3]<-1/365*7
EPDMO2$ALQ120UU[EPDMO2$ALQ120U==7]<-NA
EPDMO2$ALQ120UU[EPDMO2$ALQ120U==9]<-NA
EPDMO2$drink2<-(EPDMO2$ALQ120Q)*(EPDMO2$ALQ120UU)
EPDMO2$drink2[EPDMO2$ALQ120Q==0] <-0
EPDMO2$DRINK4[EPDMO2$ALQ120Q==0] <-0
EPDMO2$DRINK4[EPDMO2$drink2>0&EPDMO2$drink2<=4] <-1
EPDMO2$DRINK4[EPDMO2$drink2>4] <-2
EPDMO2$DRINK4[EPDMO2$ALQ111==2] <-0
EPDMO2$DRINK4[EPDMO2$ALQ121>3&EPDMO2$ALQ121<11] <-1
EPDMO2$DRINK4[EPDMO2$ALQ121==0] <-1
EPDMO2$DRINK4[EPDMO2$ALQ121==1] <-2
EPDMO2$DRINK4[EPDMO2$ALQ121==2] <-2



EPDMO2$DRINK4<-as.numeric(EPDMO2$DRINK4)

######################################################################### 整理吸烟协变量
EPDMO2$smoke[EPDMO2$SMQ020==1&EPDMO2$SMQ040==3]<-3
EPDMO2$smoke[EPDMO2$SMQ020==1&EPDMO2$SMQ040==1]<-1
EPDMO2$smoke[EPDMO2$SMQ020==1&EPDMO2$SMQ040==2]<-1
EPDMO2$smoke[EPDMO2$SMQ020==2]<-2
EPDMO2$smoke<-as.factor(EPDMO2$smoke)
###################################################

EPDMO21<-EPDMO2[complete.cases(EPDMO2$DII), ]
EPDMO22<-EPDMO21[complete.cases(EPDMO21$LBXPFOS), ]
EPDMO22<-EPDMO22[complete.cases(EPDMO22$LBXPFDE), ]
EPDMO23<-EPDMO22[EPDMO22$RIDAGEYR>=20,]


del <- which(EPDMO23$RIDEXPRG==1)
EPDMO24<-EPDMO23[-del, ]
EPDMO250<-EPDMO24[EPDMO24$TKCAL>800, ]
EPDMO250<-EPDMO250[EPDMO250$TKCAL<5000, ]
EPDMO251<-EPDMO250[complete.cases(EPDMO250$DMDEDUC2), ]
EPDMO25<-EPDMO251[complete.cases(EPDMO251$RIDRETH1), ]
EPDMO26<-EPDMO25[complete.cases(EPDMO25$BMXBMI), ]
EPDMO27<-EPDMO26[complete.cases(EPDMO26$smoke), ]
EPDMO28<-EPDMO27[complete.cases(EPDMO27$LBXCOT), ]
EPDMO29<-EPDMO28[complete.cases(EPDMO28$TKCAL), ]
EPDMO210<-EPDMO29[complete.cases(EPDMO29$INDFMPIR), ]
EPDMO211<-EPDMO210[complete.cases(EPDMO210$DRINK4), ]
table(EPDMO29$SDDSRVYR)
summary(is.na(EPDMO28$LBXPFHS))
summary(is.na(EPDMO28$LBXPFOS))
summary(is.na(EPDMO28$LBXPFDE))
summary(is.na(EPDMO28$LBXPFOA))
summary(is.na(EPDMO28$LBXPFNA))
nrow(EPDMO23)-nrow(EPDMO22)
nrow(EPDMO24)-nrow(EPDMO23)
nrow(EPDMO250)-nrow(EPDMO24)
nrow(EPDMO29)
nrow(EPDMO250)-nrow(EPDMO251)
nrow(EPDMO251)-nrow(EPDMO25)
nrow(EPDMO25)-nrow(EPDMO26)
nrow(EPDMO26)-nrow(EPDMO27)
nrow(EPDMO27)-nrow(EPDMO28)
nrow(EPDMO28)-nrow(EPDMO29)
nrow(EPDMO29)-nrow(EPDMO210)
nrow(EPDMO210)-nrow(EPDMO211)

##################################################################
d<-EPDMO29  ####change model1(29) or model 2(211)
d$LBXPFAS<-d$LBXPFOA+d$LBXPFOS+d$LBXPFNA+d$LBXPFHS+d$LBXPFDE

d$AGE[d$RIDAGEYR>=60]<-2
d$AGE[d$RIDAGEYR<60]<-1
d$AGE<-as.factor(d$AGE)
d$PIR[d$INDFMPIR<=1.3]<-1
d$PIR[d$INDFMPIR>1.3&d$INDFMPIR<=3.5]<-2
d$PIR[d$INDFMPIR>3.5]<-3
d$BMI[d$BMXBMI<25]<-1
d$BMI[d$BMXBMI>=25&d$BMXBMI<30]<-2
d$BMI[d$BMXBMI>=30]<-3
d$RIDRETH1<-as.factor(d$RIDRETH1)
PCT2<-paste0("LBX", c("PFOA", "PFNA","PFDE","PFOS","PFHS","PFAS")) 
des4<-nhanes_survey_design(d, "WTS")
dive4<-function(varss,design,data) {                   ###定义分类函数，根据4分位分类为1,2,3,4,，并命名为PCT2中的名字+4，即DII4
  varc=paste0("~",varss)
  vart=as.formula(varc)
  varx=paste0(varss,"4")
  testfor3<-oldsvyquantile(vart, design, c(.333,.6666),na.rm =T)
  testfor5<-cbind("-1",testfor3,"1000000000")
  varn=as.numeric(unlist(data[varss]))
  x1<-cut(varn,breaks = testfor5,labels = c(1,2,3),include_lowest=T)
  x2<-as.data.frame(x1)
  names(x2)=varx
  return(x2)
}
dive2<-function(varss,design,data) {                   ###定义分类函数，根据4分位分类为1,2,3,4,，并命名为PCT2中的名字+4，即DII4
  varc=paste0("~",varss)
  vart=as.formula(varc)
  varx=paste0(varss,"4")
  testfor3<-oldsvyquantile(vart, design, .5,na.rm =T)
  testfor5<-cbind("-1",testfor3,"1000000000")
  varn=as.numeric(unlist(data[varss]))
  x1<-cut(varn,breaks = testfor5,labels = c(1,2),include_lowest=T)
  x2<-as.data.frame(x1)
  names(x2)=varx
  return(x2)
}
LOGadjust<-function(varss,data) {
  varc<-paste0(varss,"L")
  x1<-data[,varss]
  x2<-log(x1)
  x2=as.data.frame(as.numeric(x2))
  names(x2)=varc
  return(x2)
}

#此处修改
NUM<-function(varss,data) {
  varc<-paste0(varss,"N")
  x1<-data[,varss]
  x2=as.data.frame(as.numeric(x1))
  names(x2)=varc
  return(x2)
  
}
MARKER<-c("LBDLYMNO","LBDNENO","LBXSTB","LBXSAL","LBXSAPSI","LBXSIR")

d$LBXSTB[d$LBXSTB==0]<-0.00001
DI<-c("DII","TKCAL","DII1","DII2")
PCT42<-paste0(PCT2,"4")
########## 定义PCT2后面+4为PCT42，并运行分类函数
DVC41<-lapply(PCT2, function(x) dive4(design=des4,data=d, varss = x))
newdata1 <- do.call(cbind, DVC41)
d<-cbind(d,newdata1) 
DVC41<-lapply(DI, function(x) dive4(design=des4,data=d, varss = x))
newdata1 <- do.call(cbind, DVC41)
d<-cbind(d,newdata1)    
DVC41<-lapply(PCT42, function(x) NUM(data=d, varss = x))
newdata1 <- do.call(cbind, DVC41)
d<-cbind(d,newdata1)  
DVC41<-lapply(PCT2, function(x) LOGadjust(data=d, varss = x))
newdata1 <- do.call(cbind, DVC41)
d<-cbind(d,newdata1)  ###合并数据
DVC41<-lapply(MARKER, function(x) LOGadjust(data=d, varss = x))
newdata1 <- do.call(cbind, DVC41)
d<-cbind(d,newdata1)  ###合并数据
PCT42N<-paste0(PCT42,"N")
PCT2L<-paste0(PCT2,"L")  
MARKERL<-paste0(MARKER,"L")  
#d$LBXSASSIL<-d$LBXSASSI*2
d$DII3[d$DII<=1.238328] <-0
d$DII3[d$DII>1.238328] <-1
e<-d[d$DII4==1,]
f<-d[d$DII4==2,]
g<-d[d$DII4==3,]

des4<-nhanes_survey_design(d, "WTS") 
des5<-nhanes_survey_design(e, "WTS") 
des6<-nhanes_survey_design(f, "WTS") 
des7<-nhanes_survey_design(g, "WTS") 
des8<-nhanes_survey_design(h, "WTS") 


write.csv(d,file = "E:/BaiduSyncdisk/nhanes2023/DATA MODEL1.csv",
          quote = F,row.names = T)
############################################################# 
logitMultiVar <- function(dess, group, var, adjvar,dat,digit = 3){            
  if(length(adjvar) == 1){
    formu <- as.formula(paste0(group, " ~ ", var, "+", adjvar))
  }else{
    formu <- as.formula(paste0(group, " ~ ", var, "+", paste(adjvar, collapse = "+")))
  }
  dat[[group]] <- as.factor(dat[[group]])
  subgroup <- levels(as.factor(dat[[group]]))
  subgroup1 <- paste0(subgroup[2], " vs ", subgroup[1])
  fit <- svyglm(formu, design =dess, family = quasibinomial(),control=list(maxit=100))
  unisum <- summary(fit)
  OR <- exp(coef(fit))[2:4]
  OR <- round(OR, digit=3)
  ci <- exp(confint(fit))[2:4,1:2]
  ci <- round(ci, digit=3)
  ci<-as.data.frame(ci)
  p <- unisum$coefficients[2:4,"Pr(>|t|)"]
  p <- ifelse(p < 0.001, "< 0.001", round(p, 3))
  vartcbbb <- names(exp(coef(fit)))[2:4]
  OR<-as.data.frame(OR)
  p<-as.data.frame(p)
  vartcbbb2<-as.data.frame(vartcbbb)
  row.names(vartcbbb2)<-vartcbbb
  result <- cbind(vartcbbb2,  OR, ci, p)
  names(result) <- c("var", "OR", "OR1","OR2", "p.val")
  return(result)
}
logitMultiVar2<- function(dess, group, var, adjvar,dat,digit = 3){    #线性回归函数
  if(length(adjvar) == 1){
    formu <- as.formula(paste0(group, " ~ ", var, "+", adjvar))
  }else{
    formu <- as.formula(paste0(group, " ~ ", var, "+", paste(adjvar, collapse = "+")))
  }
  dat[[group]] <- as.factor(dat[[group]])
  subgroup <- levels(as.factor(dat[[group]]))
  subgroup1 <- paste0(subgroup[2], " vs ", subgroup[1])
  fit <- svyglm(formu, design =dess,control=list(maxit=100),na.action = na.omit)
  unisum <- summary(fit)
  OR <- coef(fit)[2:4]
  OR <- round(OR, digit=3)
  ci <- confint(fit)[2:4,1:2]
  ci <- round(ci, digit=3)
  ci<-as.data.frame(ci)
  p <- unisum$coefficients[2:4,"Pr(>|t|)"]
  p <- ifelse(p < 0.001, "< 0.001", round(p, 3))
  vartcbbb <- names(coef(fit))[2:4]
  OR<-as.data.frame(OR)
  p<-as.data.frame(p)
  vartcbbb2<-as.data.frame(vartcbbb)
  row.names(vartcbbb2)<-vartcbbb
  result <- cbind(vartcbbb2,  OR, ci, p)
  names(result) <- c("var", "B", "B1","B2", "p.val")
  return(result)
}
logitMultiVar3<- function(dess, group, var, adjvar,dat,digit = 3){    #线性回归函数
    formu <- as.formula(paste0(group, " ~ ", var, "+", paste(adjvar, collapse = "+")))
  fit <- svyglm(formu, design =dess,control=list(maxit=100),na.action = na.omit)
  unisum <- summary(fit)
  OR <- coef(fit)[2]
  OR <- round(OR, digit=3)
  ci <- confint(fit)[2,1:2]
  ci <- round(ci, digit=3)
  ci<-t(as.data.frame(ci))
  p <- unisum$coefficients[2,"Pr(>|t|)"]
  p <- ifelse(p < 0.001, "< 0.001", round(p, 3))
  vartcbbb <- names(coef(fit))[2]
  OR<-as.data.frame(OR)
  p<-as.data.frame(p)
  vartcbbb2<-as.data.frame(vartcbbb)
  row.names(vartcbbb2)<-vartcbbb
  result <- cbind(vartcbbb2,  OR, ci, p)
  names(result) <- c("var", "B", "B1","B2", "p.val")
  return(result)
} 
logitMultiVar4<- function(dess, group, var, adjvar,dat,digit = 3){    #线性回归函数
  if(length(adjvar) == 1){
    formu <- as.formula(paste0(group, " ~ ", var, "+", adjvar))
  }else{
    formu <- as.formula(paste0(group, " ~ ", var, "+", paste(adjvar, collapse = "+")))
  }
  dat[[group]] <- as.factor(dat[[group]])
  subgroup <- levels(as.factor(dat[[group]]))
  subgroup1 <- paste0(subgroup[2], " vs ", subgroup[1])
  fit <- svyglm(formu, design =dess,control=list(maxit=100),na.action = na.omit)
  unisum <- summary(fit)
  OR <- coef(fit)[2]
  OR <- round(OR, digit=3)
  ci <- confint(fit)[2,1:2]
  ci <- round(ci, digit=3)
  ci<-t(as.data.frame(ci))
  p <- unisum$coefficients[2,"Pr(>|t|)"]
  p <- round(p, digit=3)
  vartcbbb <- names(coef(fit))[2]
  OR<-as.data.frame(OR)
  p<-as.data.frame(p)
  vartcbbb2<-as.data.frame(vartcbbb)
  row.names(vartcbbb2)<-vartcbbb
  result <- cbind(vartcbbb2,  OR, ci, p)
  names(result) <- c("var", "B", "B1","B2", "p.val")
  return(result)
} 


###########STEP2 运行




stat<-function(var,dat,digit = 3){    #线性回归函数
  varc<-paste0(var,"S")
  A1<-dat[,var]
  A1 <- ifelse(A1 == 0, 0.00000001, A1)
  A2<-as.data.frame(log(A1))
  names(A2)<-varc
  return(A2)
  }


des4<-nhanes_survey_design(d, "WTS")


DVC41<-lapply(NUT3, function(x) dive2(design=des4,data=d, varss = x))
newdata1 <- do.call(cbind, DVC41)
d<-cbind(d,newdata1)
NUT34<-paste0(NUT3,"4")

DVC41<-lapply(NUT3, function(x) stat(dat=d, var=x))
newdata1 <- do.call(cbind, DVC41)
d<-cbind(d,newdata1)  ###合并数据
NUT3S<-paste0(NUT3,"S")

des4<-nhanes_survey_design(d, "WTS") 
d$LBDPFAS2<-d$LBDPFDEL+d$LBDPFNAL+d$LBDPFOAL+d$LBDPFNAL+d$LBDPFOSL
table(d$LBDPFAS2)
d$LBDPFAS2<-as.factor(d$LBDPFAS2)
des4<-nhanes_survey_design(d, "WTS") 
prop<-(1-svymean(~LBDPFAS2,des4,na.rm=T))*100
prop
PFAS<-oldsvyquantile(~LBXPFAS, des4, c(.25,.50,.75),na.rm =T)
PFAS
summary(d$DII)



##########################################################End of data clean
########################################################################PART1 linear regression DII vs PFAS
adjvar32<-c("RIDAGEYR","RIAGENDR","RIDRETH1","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
logitResM1<-lapply(PCT2L, function(x) logitMultiVar3(des4, group = x, var ="DII",dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
row.names(logitResMDat1)<-PCT2L

logitResM1<-lapply(PCT2L, function(x) logitMultiVar2(des4, group = x, var ="DII4",dat=d,adjvar=adjvar32))
logitResMDat2 <- do.call(rbind, logitResM1)
logitResM1<-lapply(PCT2L, function(x) logitMultiVar3(des4, group = x, var ="as.numeric(DII4)",dat=d,adjvar=adjvar32))
logitResMDat3 <- do.call(rbind, logitResM1)
logitResMDat3
logitResMDat<-rbind(logitResMDat1,logitResMDat2,logitResMDat3)
logitResMDat
summary(d$DII)

write.csv(logitResMDat,file = "E:/BaiduSyncdisk/nhanes2023/Table2 M1.csv",
          quote = F,row.names = T)
####################################################################################################### PART2 linear regression DII food vs PFAS
adjvar32<-c("RIDAGEYR","RIAGENDR","RIDRETH1","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
logitResM1<-lapply(NUT34, function(x) logitMultiVar4(des4, group = "LBXPFOAL", var =x,dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
logitResMDat0<-logitResMDat1[logitResMDat1$p.val<0.05,]
logitResM1<-lapply(NUT34, function(x) logitMultiVar4(des4, group = "LBXPFNAL", var =x,dat=d,adjvar=adjvar32))
logitResMDat12 <- do.call(rbind, logitResM1)
logitResMDat2<-logitResMDat12[logitResMDat12$p.val<0.05,]
logitResM1<-lapply(NUT34, function(x) logitMultiVar4(des4, group = "LBXPFDEL", var =x,dat=d,adjvar=adjvar32))
logitResMDat13 <- do.call(rbind, logitResM1)
logitResMDat3<-logitResMDat13[logitResMDat13$p.val<0.05,]
logitResM1<-lapply(NUT34, function(x) logitMultiVar4(des4, group = "LBXPFOSL", var =x,dat=d,adjvar=adjvar32))
logitResMDat14 <- do.call(rbind, logitResM1)
logitResMDat4<-logitResMDat14[logitResMDat14$p.val<0.05,]
logitResM1<-lapply(NUT34, function(x) logitMultiVar4(des4, group = "LBXPFHSL", var =x,dat=d,adjvar=adjvar32))
logitResMDat15 <- do.call(rbind, logitResM1)
logitResMDat5<-logitResMDat15[logitResMDat15$p.val<0.05,]
logitResM1<-lapply(NUT34, function(x) logitMultiVar4(des4, group = "LBXPFASL", var =x,dat=d,adjvar=adjvar32))
logitResMDat16 <- do.call(rbind, logitResM1)


logitResMDat1$OUTCOME<-"PFOA"
logitResMDat12$OUTCOME<-"PFNA"
logitResMDat13$OUTCOME<-"PFDA"
logitResMDat14$OUTCOME<-"PFOS"
logitResMDat15$OUTCOME<-"PFHxS"
logitResMDat16$OUTCOME<-"Sum-PFAS"
adjvar3322<-c("RIDAGEYR","RIAGENDR","RIDRETH1","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
adjvar3333<-c("RIDAGEYR","RIAGENDR","RIDRETH1","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR")

logitResMDatAAA<-rbind(logitResMDat1,logitResMDat12,logitResMDat13,
                       logitResMDat14,logitResMDat15,logitResMDat16)
logitResMDatAAA2<-logitResMDatAAA[logitResMDatAAA$var!="TALCO42",]
logitResMDatAAA2<-logitResMDatAAA2[logitResMDatAAA2$var!="TKCAL42",]
logitResM222<-lapply(PCT2L, function(x) logitMultiVar4(des4, group = x, var ="TALCO4",dat=d,adjvar=adjvar3322))
logitResMDat1222 <- do.call(rbind, logitResM222)
logitResMDat1222$OUTCOME<-c("PFOA","PFNA","PFDA","PFOS","PFHxS","Sum-PFAS")
logitResM1333<-lapply(PCT2L, function(x) logitMultiVar4(des4, group = x, var ="TKCAL4",dat=d,adjvar=adjvar3333))
logitResMDat1333 <- do.call(rbind, logitResM1333)
logitResMDat1333$OUTCOME<-c("PFOA","PFNA","PFDA","PFOS","PFHxS","Sum-PFAS")
logitResMDatALL<-rbind(logitResMDatAAA2,logitResMDat1222,logitResMDat1333)
code<-read.table("E:/BaiduSyncdisk/nhanes2023/T3CODE.csv",header=T,sep=",")  
logitResMDatALL$var <- stri_replace_all_regex(logitResMDatALL$var,code$A,code$B,vectorize=F)
logitResMDatALL2<-logitResMDatALL[,c(1:2,5:6)]

write.csv(logitResMDatALL2,file = "E:/BaiduSyncdisk/nhanes2023/DATA FIGURE2 M1.csv",
          quote = F,row.names = T)
logitResMDatALL$BETA<-paste0(logitResMDatALL$B," ","(",logitResMDatALL$B1,", ",logitResMDatALL$B2,")")
logitResMDatALL3<-logitResMDatALL[,c("OUTCOME","var","BETA","p.val")]
library(openxlsx)
library(xlsx)
write.xlsx(logitResMDatALL3,file = "E:/BaiduSyncdisk/nhanes2023/Table S1M1.xlsx",
          quote = F,row.names = T)





#############################################################################################PART3 Figure 2
logitMultiVar5<- function(dess, group, var, adjvar,dat,digit = 3){    #线性回归函数
  formu <- as.formula(paste0(group, " ~ ", var, "*","DII","+", paste(adjvar, collapse = "+")))
  fit <- svyglm(formu, design =dess,control=list(maxit=100),na.action = na.omit)
  unisum <- summary(fit)
  p <- unisum$coefficients[,"Pr(>|t|)"]
  p <- ifelse(p < 0.001, "< 0.001", round(p, 3))
  vartcbbb <- names(coef(fit))
  p<-as.data.frame(p)
  vartcbbb2<-as.data.frame(vartcbbb)
  row.names(vartcbbb2)<-vartcbbb
  names(p) <- c( "Interaction P")
  p = tail(p, n = 1)
  return(p)
} 
  MARKERL<-c("LBDLYMNOL","LBDNENOL","LBXSTBL","LBXSAL","LBXSAPSIL","LBXSIRL")
Markers<-c("ALC","","ANC","","TBIL","","ALB","","ALP","","IRON")
logitResM0<-lapply(MARKERL, function(x) logitMultiVar4(des4, group = x, var ="LBXPFOSL",dat=d,adjvar=adjvar32))
logitResMDat0 <- do.call(rbind, logitResM0)
logitResM1<-lapply(MARKERL, function(x) logitMultiVar4(des5, group = x, var ="LBXPFOSL",dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
logitResM2<-lapply(MARKERL, function(x) logitMultiVar4(des6, group = x, var ="LBXPFOSL",dat=d,adjvar=adjvar32))
logitResMDat2 <- do.call(rbind, logitResM2)
logitResM3<-lapply(MARKERL, function(x) logitMultiVar4(des7, group = x, var ="LBXPFOSL",dat=d,adjvar=adjvar32))
logitResMDat3 <- do.call(rbind, logitResM3)
logitResM4<-lapply(MARKERL, function(x) logitMultiVar5(des4, group = x, var ="LBXPFOSL",dat=d,adjvar=adjvar32))
logitResMDat4 <- do.call(rbind, logitResM4)
logitResMDat4
logitResMDat<-cbind(logitResMDat0[,2:4],logitResMDat1[,2:4],logitResMDat2[,2:4],logitResMDat3[,2:4],logitResMDat4)

logitResMDat2<-logitResMDat


logitResMDat2<-rbind(logitResMDat2[1,],"",logitResMDat2[2,],"",logitResMDat2[3,],"",logitResMDat2[4,],"",
                     logitResMDat2[5,],"",logitResMDat2[6,])
logitResMDat2<-cbind(Markers,logitResMDat2)
write.csv(logitResMDat2,file = "E:/BaiduSyncdisk/nhanes2023/MAPFOSM1.csv",
          quote = F,row.names = F)


logitResM0<-lapply(MARKERL, function(x) logitMultiVar4(des4, group = x, var ="LBXPFOAL",dat=d,adjvar=adjvar32))
logitResMDat0 <- do.call(rbind, logitResM0)
logitResM1<-lapply(MARKERL, function(x) logitMultiVar4(des5, group = x, var ="LBXPFOAL",dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
logitResM2<-lapply(MARKERL, function(x) logitMultiVar4(des6, group = x, var ="LBXPFOAL",dat=d,adjvar=adjvar32))
logitResMDat2 <- do.call(rbind, logitResM2)
logitResM3<-lapply(MARKERL, function(x) logitMultiVar4(des7, group = x, var ="LBXPFOAL",dat=d,adjvar=adjvar32))
logitResMDat3 <- do.call(rbind, logitResM3)
logitResM4<-lapply(MARKERL, function(x) logitMultiVar5(des4, group = x, var ="LBXPFOAL",dat=d,adjvar=adjvar32))
logitResMDat4 <- do.call(rbind, logitResM4)
logitResMDat<-cbind(logitResMDat0[,2:4],logitResMDat1[,2:4],logitResMDat2[,2:4],logitResMDat3[,2:4],logitResMDat4)

logitResMDat2<-logitResMDat



logitResMDat2<-rbind(logitResMDat2[1,],"",logitResMDat2[2,],"",logitResMDat2[3,],"",logitResMDat2[4,],"",
                     logitResMDat2[5,],"",logitResMDat2[6,])
logitResMDat2<-cbind(Markers,logitResMDat2)

write.csv(logitResMDat2,file = "E:/BaiduSyncdisk/nhanes2023/MAPFOAM1.csv",
          quote = F,row.names = F)


logitResM0<-lapply(MARKERL, function(x) logitMultiVar4(des4, group = x, var ="LBXPFNAL",dat=d,adjvar=adjvar32))
logitResMDat0 <- do.call(rbind, logitResM0)
logitResM1<-lapply(MARKERL, function(x) logitMultiVar4(des5, group = x, var ="LBXPFNAL",dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
logitResM2<-lapply(MARKERL, function(x) logitMultiVar4(des6, group = x, var ="LBXPFNAL",dat=d,adjvar=adjvar32))
logitResMDat2 <- do.call(rbind, logitResM2)
logitResM3<-lapply(MARKERL, function(x) logitMultiVar4(des7, group = x, var ="LBXPFNAL",dat=d,adjvar=adjvar32))
logitResMDat3 <- do.call(rbind, logitResM3)
logitResM4<-lapply(MARKERL, function(x) logitMultiVar5(des4, group = x, var ="LBXPFNAL",dat=d,adjvar=adjvar32))
logitResMDat4 <- do.call(rbind, logitResM4)
logitResMDat<-cbind(logitResMDat0[,2:4],logitResMDat1[,2:4],logitResMDat2[,2:4],logitResMDat3[,2:4],logitResMDat4)

logitResMDat2<-logitResMDat


logitResMDat2<-rbind(logitResMDat2[1,],"",logitResMDat2[2,],"",logitResMDat2[3,],"",logitResMDat2[4,],"",
                     logitResMDat2[5,],"",logitResMDat2[6,])
logitResMDat2<-cbind(Markers,logitResMDat2)

write.csv(logitResMDat2,file = "E:/BaiduSyncdisk/nhanes2023/MAPFNAM1.csv",
          quote = F,row.names = F)
logitResM0<-lapply(MARKERL, function(x) logitMultiVar4(des4, group = x, var ="LBXPFASL",dat=d,adjvar=adjvar32))
logitResMDat0<- do.call(rbind, logitResM0)
logitResM1<-lapply(MARKERL, function(x) logitMultiVar4(des5, group = x, var ="LBXPFASL",dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
logitResM2<-lapply(MARKERL, function(x) logitMultiVar4(des6, group = x, var ="LBXPFASL",dat=d,adjvar=adjvar32))
logitResMDat2 <- do.call(rbind, logitResM2)
logitResM3<-lapply(MARKERL, function(x) logitMultiVar4(des7, group = x, var ="LBXPFASL",dat=d,adjvar=adjvar32))
logitResMDat3 <- do.call(rbind, logitResM3)
logitResM4<-lapply(MARKERL, function(x) logitMultiVar5(des4, group = x, var ="LBXPFASL",dat=d,adjvar=adjvar32))
logitResMDat4 <- do.call(rbind, logitResM4)
logitResMDat<-cbind(logitResMDat0[,2:4],logitResMDat1[,2:4],logitResMDat2[,2:4],logitResMDat3[,2:4],logitResMDat4)

logitResMDat2<-logitResMDat



logitResMDat2<-rbind(logitResMDat2[1,],"",logitResMDat2[2,],"",logitResMDat2[3,],"",logitResMDat2[4,],"",
                     logitResMDat2[5,],"",logitResMDat2[6,])
logitResMDat2<-cbind(Markers,logitResMDat2)

write.csv(logitResMDat2,file = "E:/BaiduSyncdisk/nhanes2023/MAPFASM1.csv",
          quote = F,row.names = F)
logitResM0<-lapply(MARKERL, function(x) logitMultiVar4(des4, group = x, var ="LBXPFDEL",dat=d,adjvar=adjvar32))
logitResMDat0<- do.call(rbind, logitResM0)
logitResM1<-lapply(MARKERL, function(x) logitMultiVar4(des5, group = x, var ="LBXPFDEL",dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
logitResM2<-lapply(MARKERL, function(x) logitMultiVar4(des6, group = x, var ="LBXPFDEL",dat=d,adjvar=adjvar32))
logitResMDat2 <- do.call(rbind, logitResM2)
logitResM3<-lapply(MARKERL, function(x) logitMultiVar4(des7, group = x, var ="LBXPFDEL",dat=d,adjvar=adjvar32))
logitResMDat3 <- do.call(rbind, logitResM3)
logitResM4<-lapply(MARKERL, function(x) logitMultiVar5(des4, group = x, var ="LBXPFDEL",dat=d,adjvar=adjvar32))
logitResMDat4 <- do.call(rbind, logitResM4)
logitResMDat<-cbind(logitResMDat0[,2:4],logitResMDat1[,2:4],logitResMDat2[,2:4],logitResMDat3[,2:4],logitResMDat4)

logitResMDat2<-logitResMDat



logitResMDat2<-rbind(logitResMDat2[1,],"",logitResMDat2[2,],"",logitResMDat2[3,],"",logitResMDat2[4,],"",
                     logitResMDat2[5,],"",logitResMDat2[6,])
logitResMDat2<-cbind(Markers,logitResMDat2)

write.csv(logitResMDat2,file = "E:/BaiduSyncdisk/nhanes2023/MAPFDAM1.csv",
          quote = F,row.names = F)


logitResM0<-lapply(MARKERL, function(x) logitMultiVar4(des4, group = x, var ="LBXPFHSL",dat=d,adjvar=adjvar32))
logitResMDat0<- do.call(rbind, logitResM0)
logitResM1<-lapply(MARKERL, function(x) logitMultiVar4(des5, group = x, var ="LBXPFHSL",dat=d,adjvar=adjvar32))
logitResMDat1 <- do.call(rbind, logitResM1)
logitResM2<-lapply(MARKERL, function(x) logitMultiVar4(des6, group = x, var ="LBXPFHSL",dat=d,adjvar=adjvar32))
logitResMDat2 <- do.call(rbind, logitResM2)
logitResM3<-lapply(MARKERL, function(x) logitMultiVar4(des7, group = x, var ="LBXPFHSL",dat=d,adjvar=adjvar32))
logitResMDat3 <- do.call(rbind, logitResM3)
logitResM4<-lapply(MARKERL, function(x) logitMultiVar5(des4, group = x, var ="LBXPFHSL",dat=d,adjvar=adjvar32))
logitResMDat4 <- do.call(rbind, logitResM4)
logitResMDat<-cbind(logitResMDat0[,2:4],logitResMDat1[,2:4],logitResMDat2[,2:4],logitResMDat3[,2:4],logitResMDat4)

logitResMDat2<-logitResMDat




logitResMDat2<-rbind(logitResMDat2[1,],"",logitResMDat2[2,],"",logitResMDat2[3,],"",logitResMDat2[4,],"",
                     logitResMDat2[5,],"",logitResMDat2[6,])
logitResMDat2<-cbind(Markers,logitResMDat2)

write.csv(logitResMDat2,file = "E:/BaiduSyncdisk/nhanes2023/MAPFHSM1.csv",
          quote = F,row.names = F)
########################################################################################################################
logitMultiVar6<- function(dess, group, var,varsub, adjvar,dat,digit = 3){    #线性回归函数
  formu <- as.formula(paste0(group, " ~ ", var, "*",varsub,"+", paste(adjvar, collapse = "+")))
  fit <- svyglm(formu, design =dess,control=list(maxit=100),na.action = na.omit)
  unisum <- summary(fit)
  p <- unisum$coefficients[,"Pr(>|t|)"]
  p<-tail(p, n = 1)
  p <- ifelse(p < 0.001, "< 0.001", round(p, 3))
  vartcbbb <- names(coef(fit))
  vartcbbb<-tail(vartcbbb, n = 1)
  p<-as.data.frame(p)
  vartcbbb2<-as.data.frame(vartcbbb)
  row.names(vartcbbb2)<-vartcbbb
  names(p) <- c( "Interaction P")
  return(p)}



adjvar32<-c("RIDAGEYR","RIAGENDR","RIDRETH1","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
adjvar33<-c("RIDAGEYR","RIDRETH1","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
adjvar34<-c("RIDAGEYR","RIAGENDR","RIDRETH1","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
adjvar35<-c("RIDAGEYR","RIAGENDR","RIDRETH1","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
adjvar36<-c("RIDAGEYR","RIAGENDR","smoke","DMDEDUC2","BMXBMI","LBXCOT","SDDSRVYR","TKCAL")
PCT2LL<-paste0("LBX", c("PFOAL", "PFNAL","PFOSL","PFASL")) 
dage1<-d[d$AGE==1,]
dage2<-d[d$AGE==2,]
dsex1<-d[d$RIAGENDR==1,]
dsex2<-d[d$RIAGENDR==2,]
dPIR1<-d[d$PIR==1,]
dPIR2<-d[d$PIR==2,]
dPIR3<-d[d$PIR==3,]
dsmoke1<-d[d$smoke==1,]
dsmoke2<-d[d$smoke==2,]
dsmoke3<-d[d$smoke==3,]
dDRINK1<-d[d$DRINK4==0,]
dDRINK2<-d[d$DRINK4==1,]
dDRINK3<-d[d$DRINK4==2,]
dRACE1<-d[d$RIDRETH1==1,]
dRACE2<-d[d$RIDRETH1==3,]
dRACE3<-d[d$RIDRETH1==4,]
dRACE4<-d[d$RIDRETH1==5,]
dBMI1<-d[d$BMI==1,]
dBMI2<-d[d$BMI==2,]
dBMI3<-d[d$BMI==3,]
age1<-nhanes_survey_design(dage1, "WTS")
age2<-nhanes_survey_design(dage2, "WTS")
sex1<-nhanes_survey_design(dsex1, "WTS")
sex2<-nhanes_survey_design(dsex2, "WTS")
PIR1<-nhanes_survey_design(dPIR1, "WTS")
PIR2<-nhanes_survey_design(dPIR2, "WTS")
PIR3<-nhanes_survey_design(dPIR3, "WTS")
smoke1<-nhanes_survey_design(dsmoke1, "WTS")
smoke2<-nhanes_survey_design(dsmoke2, "WTS")
smoke3<-nhanes_survey_design(dsmoke3, "WTS")
DRINK1<-nhanes_survey_design(dDRINK1, "WTS")
DRINK2<-nhanes_survey_design(dDRINK2, "WTS")
DRINK3<-nhanes_survey_design(dDRINK3, "WTS")
RACE1<-nhanes_survey_design(dRACE1, "WTS")
RACE2<-nhanes_survey_design(dRACE2, "WTS")
RACE3<-nhanes_survey_design(dRACE3, "WTS")
RACE4<-nhanes_survey_design(dRACE4, "WTS")
BMI1<-nhanes_survey_design(dBMI1, "WTS")
BMI2<-nhanes_survey_design(dBMI2, "WTS")
BMI3<-nhanes_survey_design(dBMI3, "WTS")

logitResMage1<-lapply(PCT2LL, function(x) logitMultiVar3(age1, group=x, var ="DII",dat=d,adjvar=adjvar32))
age1 <- do.call(rbind, logitResMage1)


logitResMage2<-lapply(PCT2LL, function(x) logitMultiVar3(age2, group=x, var ="DII",dat=d,adjvar=adjvar32))
age2 <- do.call(rbind, logitResMage2)

logitResMsex1<-lapply(PCT2LL, function(x) logitMultiVar3(sex1, group=x, var ="DII",dat=d,adjvar=adjvar33))
sex1 <- do.call(rbind, logitResMsex1)


logitResMsex2<-lapply(PCT2LL, function(x) logitMultiVar3(sex2, group=x, var ="DII",dat=d,adjvar=adjvar33))
sex2 <- do.call(rbind, logitResMsex2)



logitResMPIR1<-lapply(PCT2LL, function(x) logitMultiVar3(PIR1, group=x, var= "DII",dat=d,adjvar=adjvar32))
PIR1 <- do.call(rbind, logitResMPIR1)
logitResMPIR2<-lapply(PCT2LL, function(x) logitMultiVar3(PIR2, group=x, var= "DII",dat=d,adjvar=adjvar32))
PIR2 <- do.call(rbind, logitResMPIR2)
logitResMPIR3<-lapply(PCT2LL, function(x) logitMultiVar3(PIR3, group=x, var= "DII",dat=d,adjvar=adjvar32))
PIR3 <- do.call(rbind, logitResMPIR3)

logitResMsmoke1<-lapply(PCT2LL, function(x) logitMultiVar3(smoke1, group=x, var= "DII",dat=d,adjvar=adjvar34))
smoke1 <- do.call(rbind, logitResMsmoke1)
logitResMsmoke2<-lapply(PCT2LL, function(x) logitMultiVar3(smoke2, group=x, var= "DII",dat=d,adjvar=adjvar34))
smoke2 <- do.call(rbind, logitResMsmoke2)
logitResMsmoke3<-lapply(PCT2LL, function(x) logitMultiVar3(smoke3, group=x, var= "DII",dat=d,adjvar=adjvar34))
smoke3 <- do.call(rbind, logitResMsmoke3)

logitResMDRINK1<-lapply(PCT2LL, function(x) logitMultiVar3(DRINK1, group=x, var= "DII",dat=d,adjvar=adjvar35))
DRINK1 <- do.call(rbind, logitResMDRINK1)
logitResMDRINK2<-lapply(PCT2LL, function(x) logitMultiVar3(DRINK2, group=x, var= "DII",dat=d,adjvar=adjvar35))
DRINK2 <- do.call(rbind, logitResMDRINK2)
logitResMDRINK3<-lapply(PCT2LL, function(x) logitMultiVar3(DRINK3, group=x, var= "DII",dat=d,adjvar=adjvar35))
DRINK3 <- do.call(rbind, logitResMDRINK3)

logitResMRACE1<-lapply(PCT2LL, function(x) logitMultiVar3(RACE1, group=x, var= "DII",dat=d,adjvar=adjvar36))
RACE1 <- do.call(rbind, logitResMRACE1)
logitResMRACE2<-lapply(PCT2LL, function(x) logitMultiVar3(RACE2, group=x, var= "DII",dat=d,adjvar=adjvar36))
RACE2 <- do.call(rbind, logitResMRACE2)
logitResMRACE3<-lapply(PCT2LL, function(x) logitMultiVar3(RACE3, group=x, var= "DII",dat=d,adjvar=adjvar36))
RACE3 <- do.call(rbind, logitResMRACE3)
logitResMRACE4<-lapply(PCT2LL, function(x) logitMultiVar3(RACE4, group=x, var= "DII",dat=d,adjvar=adjvar36))
RACE4 <- do.call(rbind, logitResMRACE4)

logitResMBMI1<-lapply(PCT2LL, function(x) logitMultiVar3(BMI1, group=x, var= "DII",dat=d,adjvar=adjvar32))
BMI1 <- do.call(rbind, logitResMBMI1)
logitResMBMI2<-lapply(PCT2LL, function(x) logitMultiVar3(BMI2, group=x, var= "DII",dat=d,adjvar=adjvar32))
BMI2 <- do.call(rbind, logitResMBMI2)
logitResMBMI3<-lapply(PCT2LL, function(x) logitMultiVar3(BMI3, group=x, var= "DII",dat=d,adjvar=adjvar32))
BMI3 <- do.call(rbind, logitResMBMI3)

PFOASUB<-rbind("",age2[1,2:4],age1[1,2:4],"",sex1[1,2:4],sex2[1,2:4],"",RACE1[1,2:4],RACE2[1,2:4],RACE3[1,2:4],RACE4[1,2:4],"",
               BMI1[1,2:4],BMI2[1,2:4],BMI3[1,2:4],"",PIR1[1,2:4],PIR2[1,2:4],PIR3[1,2:4]
               ,"",smoke1[1,2:4],smoke2[1,2:4],smoke3[1,2:4],"",DRINK1[1,2:4],DRINK2[1,2:4],DRINK3[1,2:4])

SUBVAL<-c("AGE","RIAGENDR","as.numeric(RIDRETH1)","BMI","PIR","as.numeric(smoke)","DRINK4")
SUBGROUPS<-c("Age","<60 years","≥60 years","Sex","Male","Female","Race","Hispanic","Non-Hispanic White",
             "Non-Hispanic Black","Other Race","BMI","<25.0","25.0-29.9","≥30.0","PIR","<1.3","1.3-3.5",">3.5",
             "Smoke Status","Current" ,"Never","Former","Drinking Status","Never","≤4 drinks per week",">4 drinks per week")

PFOASUB<-cbind(SUBGROUPS, PFOASUB)
logitResMINTERPFOA<-lapply(SUBVAL, function(x) logitMultiVar6(des4, varsub = x, group = "LBXPFOAL",var="DII",dat=d,adjvar=adjvar32))
ERPFOA <- do.call(rbind, logitResMINTERPFOA)
ERPFOA2<-as.data.frame(rbind(ERPFOA[1,],"","",ERPFOA[2,],"","",ERPFOA[3,],"","","","",ERPFOA[4,],"","","",ERPFOA[5,],"","",""
               ,ERPFOA[6,],"","","",ERPFOA[7,],"","",""))
names(ERPFOA2)<-"Interaction P"
PFOASUB2<-cbind(PFOASUB,ERPFOA2)
write.csv(PFOASUB2,file = "E:/BaiduSyncdisk/nhanes2023/PFOASUBM1.csv",
          quote = F,row.names = F)

PFNASUB<-rbind("",age2[2,2:4],age1[2,2:4],"",sex1[2,2:4],sex2[2,2:4],"",RACE1[2,2:4],RACE2[2,2:4],RACE3[2,2:4],RACE4[2,2:4],"",
               BMI1[2,2:4],BMI2[2,2:4],BMI3[2,2:4],"",PIR1[2,2:4],PIR2[2,2:4],PIR3[2,2:4]
               ,"",smoke1[2,2:4],smoke2[2,2:4],smoke3[2,2:4],"",DRINK1[2,2:4],DRINK2[2,2:4],DRINK3[2,2:4])


PFNASUB<-cbind(SUBGROUPS, PFNASUB)
logitResMINTERPFNA<-lapply(SUBVAL, function(x) logitMultiVar6(des4, varsub = x, group = "LBXPFNAL",var="DII",dat=d,adjvar=adjvar32))
ERPFNA <- do.call(rbind, logitResMINTERPFNA)
ERPFNA2<-as.data.frame(rbind(ERPFNA[1,],"","",ERPFNA[2,],"","",ERPFNA[3,],"","","","",ERPFNA[4,],"","","",ERPFNA[5,],"","",""
                             ,ERPFNA[6,],"","","",ERPFNA[7,],"","",""))
names(ERPFNA2)<-"Interaction P"
PFNASUB2<-cbind(PFNASUB,ERPFNA2)
write.csv(PFNASUB2,file = "E:/BaiduSyncdisk/nhanes2023/PFNASUBM1.csv",
          quote = F,row.names = F)


PFOSSUB<-rbind("",age2[3,2:4],age1[3,2:4],"",sex1[3,2:4],sex2[3,2:4],"",RACE1[3,2:4],RACE2[3,2:4],RACE3[3,2:4],RACE4[3,2:4],"",
               BMI1[3,2:4],BMI2[3,2:4],BMI3[3,2:4],"",PIR1[3,2:4],PIR2[3,2:4],PIR3[3,2:4]
               ,"",smoke1[3,2:4],smoke2[3,2:4],smoke3[3,2:4],"",DRINK1[3,2:4],DRINK2[3,2:4],DRINK3[3,2:4])
PFOSSUB<-cbind(SUBGROUPS, PFOSSUB)
logitResMINTERPFOS<-lapply(SUBVAL, function(x) logitMultiVar6(des4, varsub = x, group = "LBXPFOSL",var="DII",dat=d,adjvar=adjvar32))
ERPFOS <- do.call(rbind, logitResMINTERPFOS)
ERPFOS2<-as.data.frame(rbind(ERPFOS[1,],"","",ERPFOS[2,],"","",ERPFOS[3,],"","","","",ERPFOS[4,],"","","",ERPFOS[5,],"","",""
                             ,ERPFOS[6,],"","","",ERPFOS[7,],"","",""))
names(ERPFOS2)<-"Interaction P"
PFOSSUB2<-cbind(PFOSSUB,ERPFOS2)

write.csv(PFOSSUB2,file = "E:/BaiduSyncdisk/nhanes2023/PFOSSUBM1.csv",
          quote = F,row.names = F)



PFASSUB<-rbind("",age2[4,2:4],age1[4,2:4],"",sex1[4,2:4],sex2[4,2:4],"",RACE1[4,2:4],RACE2[4,2:4],RACE3[4,2:4],RACE4[4,2:4],"",
               BMI1[4,2:4],BMI2[4,2:4],BMI3[4,2:4],"",PIR1[4,2:4],PIR2[4,2:4],PIR3[4,2:4]
               ,"",smoke1[4,2:4],smoke2[4,2:4],smoke3[4,2:4],"",DRINK1[4,2:4],DRINK2[4,2:4],DRINK3[4,2:4])
PFASSUB<-cbind(SUBGROUPS, PFASSUB)
logitResMINTERPFAS<-lapply(SUBVAL, function(x) logitMultiVar6(des4, varsub = x, group = "LBXPFASL",var="DII",dat=d,adjvar=adjvar32))
ERPFAS <- do.call(rbind, logitResMINTERPFAS)
ERPFAS2<-as.data.frame(rbind(ERPFAS[1,],"","",ERPFAS[2,],"","",ERPFAS[3,],"","","","",ERPFAS[4,],"","","",ERPFAS[5,],"","",""
                             ,ERPFAS[6,],"","","",ERPFAS[7,],"","",""))
names(ERPFAS2)<-"Interaction P"
PFASSUB2<-cbind(PFASSUB,ERPFAS2)

write.csv(PFASSUB2,file = "E:/BaiduSyncdisk/nhanes2023/PFASSUBM1.csv",
          quote = F,row.names = F)

##############################################################################################################
library(epiDisplay)
library(corrplot)
library(Hmisc)
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat) 
  data.frame( row = rownames(cormat)[row(cormat)[ut]], 
              column = rownames(cormat)[col(cormat)[ut]], cor =(cormat)[ut], p = pmat[ut] )
}
d2<-d[,PCT2L]
names(d2)<- c("ln(PFOA)", "ln(PFNA)","ln(PFDA)","ln(PFOS)","ln(PFHxS)","ln(Sum-PFAS)")
b2<-asNumericMatrix(d2)
flattenCorrMatrix(res3$r, res3$P)
coraa<-rcorr(b2, type ="pearson")
cortext<-flattenCorrMatrix(coraa$r, coraa$P)
cortext$P<-round(cortext$P,digits=3)


colors2 <- colorRampPalette(c("#32B897", "#FFA500", "#F27970"))(30)
corraar<-coraa$r
corraar[corraar < 0] <- 0
pdf("E:/BaiduSyncdisk/nhanes2023/FIGURES1M1.pdf",width =16, height =10)
corrplot(corr=corraar,method = "ellipse",type="upper",tl.pos = "d",p.mat = coraa$P, sig.level = .01,insig = "blank",is.corr=T,cl.lim=c(0, 1),
         col = colors2,tl.col="black",tl.cex=1.2)
corrplot(corr=corraar,p.mat = coraa$P, sig.level = .01,add=TRUE, type="lower", method="number",order="AOE",
         diag=FALSE,tl.pos="n", cl.pos="n",insig = "blank",col = colors2,number.cex = 1.5)
dev.off()

