library(foreign)
library(survey)
#setwd("~/PFOA")
#1999-2000
#baseurl00<-"ftp://ftp.cdc.gov/pub/Health_Statistics/NCHS/nhanes/1999-2000"
#files<-c("SSPFC_A.xpt","Lab13.xpt","SMQ.xpt","DEMO.xpt","BPX.xpt","LAB10.xpt","MCQ.xpt","LEXABPI.xpt")
#lapply(files, function(f) if (!file.exists(f)) download.file(paste(baseurl00,f,sep="/"),f))

pfc00<-read.xport("SSPFC_A.xpt")
chol00<-read.xport("Lab13.xpt")
smok00<-read.xport("SMQ.xpt")
demo00<-read.xport("DEMO.xpt")
bp00<-read.xport("BPX.xpt")
lab00<-read.xport("LAB10.xpt")
cond00<-read.xport("MCQ.xpt")
abi00<-read.xport("LEXABPI.xpt")
alldata00<-merge(pfc00,merge(chol00,merge(smok00,merge(demo00,merge(bp00,merge(lab00,merge(cond00,abi00)))))))

nhanes00<-svydesign(id=~SDMVPSU,strata=~SDMVSTRA,weights=~WTMEC2YR,data=alldata00,nest=TRUE)

nhanes00<-update(nhanes00, hadcvd=(MCQ160C==1) | (MCQ160E==1) | (MCQ160F==1), abi=(LEXRABPI+LEXLABPI)/2)
nhanes00<-update(nhanes00, haspad=ifelse(abi<0.9,1,ifelse(abi>1.5,NA,0)))
svyquantile(~SPFOA,nhanes00,quantiles=c(0.25,0.5,0.75))
nhanes00<-update(nhanes00, pfoa4=cut(SPFOA,c(0,3.7,5,6.8,Inf)))
nhanes00<-update(nhanes00, smoking=ifelse(SMQ020==2,0,ifelse(SMQ040 %in% c(1,2),1,2)))

mcvd00<-svyglm(hadcvd~pfoa4+DMDEDUC+LBXGH+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+BPXSAR+BPXDAR+factor(smoking)+LBXTC,family=binomial,design=nhanes00)
summary(mcvd00,df=degf(nhanes00))

mpad00<-svyglm(haspad~pfoa4+DMDEDUC+LBXGH+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+BPXSAR+BPXDAR+factor(smoking)+LBXTC,family=binomial,design=nhanes00)
summary(mpad00,df=degf(nhanes00))

mcvd00u<-svyglm(hadcvd~pfoa4+DMDEDUC+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+factor(smoking),family=binomial,design=nhanes00)

##2003-2004
##baseurl04<-"ftp://ftp.cdc.gov/pub/Health_Statistics/NCHS/nhanes/2003-2004"
##files<-c("L24PFC_C.xpt","l13_c.xpt","SMQ_C.xpt","demo_c.xpt","BPX_C.xpt","L10_C.xpt","MCQ_C.xpt","LEXAB_C.xpt")
##lapply(files, function(f) if (!file.exists(f)) download.file(paste(baseurl04,f,sep="/"),f))

pfc04<-read.xport("L24PFC_C.xpt")
chol04<-read.xport("l13_c.xpt")
smok04<-read.xport("SMQ_C.xpt")
demo04<-read.xport("demo_c.xpt")
bp04<-read.xport("BPX_C.xpt")
lab04<-read.xport("L10_C.xpt")
cond04<-read.xport("MCQ_C.xpt")
abi04<-read.xport("LEXAB_C.xpt")
alldata04<-merge(pfc04,merge(chol04,merge(smok04,merge(demo04,merge(bp04,merge(lab04,merge(cond04,abi04)))))))
nhanes04<-svydesign(id=~SDMVPSU,strata=~SDMVSTRA,weights=~WTSA2YR,data=alldata04,nest=TRUE)


nhanes04<-update(nhanes04, hadcvd=(MCQ160C==1) | (MCQ160E==1) | (MCQ160F==1), abi=(LEXRABPI+LEXLABPI)/2)
nhanes04<-update(nhanes04, haspad=ifelse(abi<0.85,1,ifelse(abi>1.5,NA,0)))
svyquantile(~LBXPFOA,nhanes04,quantiles=c(0.25,0.5,0.75),na.rm=TRUE)
nhanes04<-update(nhanes04, pfoa4=cut(LBXPFOA,c(0,2.8,4.2,6,Inf)))
nhanes04<-update(nhanes04, smoking=ifelse(SMQ020==2,0,ifelse(SMQ040 %in% c(1,2),1,2)))


mcvd04<-svyglm(hadcvd~pfoa4+DMDEDUC+LBXGH+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+BPXSAR+BPXDAR+factor(smoking)+LBXTC,family=binomial,design=nhanes04)
summary(mcvd04,df=degf(nhanes04))

mpad04<-svyglm(haspad~pfoa4+DMDEDUC+LBXGH+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+BPXSAR+BPXDAR+factor(smoking)+LBXTC,family=binomial,design=nhanes04)
summary(mpad04,df=degf(nhanes04))
mcvd04u<-svyglm(hadcvd~pfoa4+DMDEDUC+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+factor(smoking),family=binomial,design=nhanes04)


####2005-2006
##baseurl06<-"ftp://ftp.cdc.gov/pub/Health_Statistics/NCHS/nhanes/2005-2006"
##files<-c("PFC_D.xpt","TCHOL_d.xpt","SMQ_d.xpt","demo_d.xpt","BPX_D.xpt","GHB_D.xpt","MCQ_D.xpt")
##lapply(files, function(f) if (!file.exists(f)) download.file(paste(baseurl06,f,sep="/"),f))

pfc06<-read.xport("PFC_D.xpt")
chol06<-read.xport("TCHOL_d.xpt")
smok06<-read.xport("SMQ_d.xpt")
demo06<-read.xport("demo_d.xpt")
bp06<-read.xport("BPX_D.xpt")
lab06<-read.xport("GHB_D.xpt")
cond06<-read.xport("MCQ_d.xpt")
alldata06<-merge(pfc06,merge(chol06,merge(smok06,merge(demo06,merge(bp06,merge(lab06,cond06))))))
nhanes06<-svydesign(id=~SDMVPSU,strata=~SDMVSTRA,weights=~WTSA2YR,data=alldata06,nest=TRUE)


nhanes06<-update(nhanes06, hadcvd=(MCQ160C==1) | (MCQ160E==1) | (MCQ160F==1))
svyquantile(~LBXPFOA,nhanes06,quantiles=c(0.25,0.5,0.75),na.rm=TRUE)
nhanes06<-update(nhanes06, pfoa4=cut(LBXPFOA,c(0,2.6,4.2,6.2,Inf)))
nhanes06<-update(nhanes06, smoking=ifelse(SMQ020==2,0,ifelse(SMQ040 %in% c(1,2),1,2)))


mcvd06<-svyglm(hadcvd~pfoa4+DMDEDUC2+LBXGH+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+BPXSY1+BPXDI1+factor(smoking)+LBXTC,family=binomial,design=nhanes06)
summary(mcvd06,df=degf(nhanes06))

mcvd06u<-svyglm(hadcvd~pfoa4+DMDEDUC2+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+factor(smoking),family=binomial,design=nhanes06)


##2007-2008

##baseurl08<-"ftp://ftp.cdc.gov/pub/Health_Statistics/NCHS/nhanes/2007-2008"
##files<-c("PFC_E.xpt","TCHOL_e.xpt","SMQ_e.xpt","demo_e.xpt","BPX_E.xpt","GHB_E.xpt","MCQ_E.xpt")
##lapply(files, function(f) if (!file.exists(f)) download.file(paste(baseurl08,f,sep="/"),f))

pfc08<-read.xport("PFC_E.xpt")
chol08<-read.xport("TCHOL_e.xpt")
smok08<-read.xport("SMQ_e.xpt")
demo08<-read.xport("demo_e.xpt")
bp08<-read.xport("BPX_E.xpt")
lab08<-read.xport("GHB_E.xpt")
cond08<-read.xport("MCQ_e.xpt")
alldata08<-merge(pfc08,merge(chol08,merge(smok08,merge(demo08,merge(bp08,merge(lab08,cond08))))))
nhanes08<-svydesign(id=~SDMVPSU,strata=~SDMVSTRA,weights=~WTSC2YR,data=alldata08,nest=TRUE)


nhanes08<-update(nhanes08, hadcvd=(MCQ160C==1) | (MCQ160E==1) | (MCQ160F==1))
svyquantile(~LBXPFOA,nhanes08,quantiles=c(0.25,0.5,0.75),na.rm=TRUE)
nhanes08<-update(nhanes08, pfoa4=cut(LBXPFOA,c(0,2.6,4.2,6.2,Inf)))
nhanes08<-update(nhanes08, smoking=ifelse(SMQ020==2,0,ifelse(SMQ040 %in% c(1,2),1,2)))


mcvd08<-svyglm(hadcvd~pfoa4+DMDEDUC2+LBXGH+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+BPXSY1+BPXDI1+factor(smoking)+LBXTC,family=binomial,design=nhanes08)
summary(mcvd08,df=degf(nhanes08))

mcvd08u<-svyglm(hadcvd~pfoa4+DMDEDUC2+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+factor(smoking),family=binomial,design=nhanes08)


##2009-2010

##baseurl10<-"ftp://ftp.cdc.gov/pub/Health_Statistics/NCHS/nhanes/2009-2010"
##files<-c("PFC_F.xpt","TCHOL_f.xpt","SMQ_f.xpt","demo_f.xpt","BPX_F.xpt","GHB_F.xpt","MCQ_F.xpt")
##lapply(files, function(f) if (!file.exists(f)) download.file(paste(baseurl10,f,sep="/"),f))

pfc10<-read.xport("PFC_F.xpt")
chol10<-read.xport("TCHOL_f.xpt")
smok10<-read.xport("SMQ_f.xpt")
demo10<-read.xport("demo_f.xpt")
bp10<-read.xport("BPX_F.xpt")
lab10<-read.xport("GHB_F.xpt")
cond10<-read.xport("MCQ_f.xpt")
alldata10<-merge(pfc10,merge(chol10,merge(smok10,merge(demo10,merge(bp10,merge(lab10,cond10))))))
nhanes10<-svydesign(id=~SDMVPSU,strata=~SDMVSTRA,weights=~WTSC2YR,data=alldata10,nest=TRUE)


nhanes10<-update(nhanes10, hadcvd=(MCQ160C==1) | (MCQ160E==1) | (MCQ160F==1))
svyquantile(~LBXPFOA,nhanes10,quantiles=c(0.25,0.5,0.75),na.rm=TRUE)
nhanes10<-update(nhanes10, pfoa4=cut(LBXPFOA,c(0,2.6,4.2,6.2,Inf)))
nhanes10<-update(nhanes10, smoking=ifelse(SMQ020==2,0,ifelse(SMQ040 %in% c(1,2),1,2)))


mcvd10<-svyglm(hadcvd~pfoa4+DMDEDUC2+LBXGH+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+BPXSY1+BPXDI1+factor(smoking)+LBXTC,family=binomial,design=nhanes10)
summary(mcvd10,df=degf(nhanes10))

mcvd10u<-svyglm(hadcvd~pfoa4+DMDEDUC2+RIDAGEEX+factor(RIDRETH1)+RIAGENDR+factor(smoking),family=binomial,design=nhanes10)

models<-list(mcvd00,mcvd04,mcvd06,mcvd08,mcvd10)
coefs<-sapply(models,function(m) coef(m)[2:4])
ses<-sapply(models, function(m) SE(m)[2:4])

library(rmeta)
q4<-meta.summaries(coefs[3,],ses[3,],method="random",logscale=TRUE)
summary(q4)
metaplot(coefs[3,],ses[3,],col=meta.colors(box=c(2,2,1,1,1)),logeffect=TRUE)

modelsu<-list(mcvd00u,mcvd04u,mcvd06u,mcvd08u,mcvd10u)
coefsu<-sapply(modelsu,function(m) coef(m)[2:4])
sesu<-sapply(modelsu, function(m) SE(m)[2:4])
metaplot(coefsu[3,],sesu[3,],col=meta.colors(box=c(2,2,1,1,1)),logeffect=TRUE)
q4u<-meta.summaries(coefsu[3,],sesu[3,],method="random",logscale=TRUE)
summary(q4u)