#####################Spearman’s correlation analysis######################
library(corrplot)
library(psych)
library(haven)
long <- read_dta("D:/Stata17/ado/personal/long.dta")
long <- as.data.frame(long)
data1 <- subset(long,select=c(deng,pm25,pm10,qiwen,qiya,shidu))
colnames(data1) = c("ALAN","PM2.5","PM10","MeanTemp","Ap","RH")
M<- cor(data1, use="pairwise.complete.obs", method="pearson")
P<- round(as.matrix(corr.test(data1, use="pairwise.complete.obs",method = "pearson")$p),3)

tiff(file = "correlationplot.tiff", width = 3000, height = 3000, res = 300)
corrplot(M, p.mat = P, method = "circle", type = "lower",
         sig.level = c(.001, .01, .05),  tl.pos="lt", tl.col="black", tl.cex=2,  tl.offset=0.6,cl.pos="r",
         insig = "label_sig", pch.cex = 2, pch.col="red",cl.cex = 2)
corrplot(M,  type="upper", method="number",
         col="coral4",  tl.pos="n", cl.pos="n", number.cex = 1.5, add=T,diag=F)
dev.off()


##############################################################################
##############################################################################
###########Lag association between outdoor ALAN and metabolicdisease###################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)

wide_men <- read_dta("D:/Stata17/ado/personal/Data/wide_men.dta")
wide_men <- as.data.frame(wide_men)
wide_men$time = ymd(wide_men$time)
wide_men$doy <- weekdays(wide_men$time)
wide_men[,c("YZ","doy")]<-lapply(wide_men[,c("YZ","doy")],factor)

wide_women <- read_dta("D:/Stata17/ado/personal/Data/wide_women.dta")
wide_women <- as.data.frame(wide_women)
wide_women$time = ymd(wide_women$time)
wide_women$doy <- weekdays(wide_women$time)
wide_women[,c("YZ","doy")]<-lapply(wide_women[,c("YZ","doy")],factor)

wide_19 <- read_dta("D:/Stata17/ado/personal/Data/wide_from_19_to_45_years.dta")
wide_19 <- as.data.frame(wide_19)
wide_19$time = ymd(wide_19$time)
wide_19$doy <- weekdays(wide_19$time)
wide_19[,c("gender","YZ","doy")]<-lapply(wide_19[,c("gender","YZ","doy")],factor)

wide_46 <- read_dta("D:/Stata17/ado/personal/Data/wide_from_46_to_59_years.dta")
wide_46 <- as.data.frame(wide_46)
wide_46$time = ymd(wide_46$time)
wide_46$doy <- weekdays(wide_46$time)
wide_46[,c("gender","YZ","doy")]<-lapply(wide_46[,c("gender","YZ","doy")],factor)

wide_60 <- read_dta("D:/Stata17/ado/personal/Data/wide_over_60_years.dta")
wide_60 <- as.data.frame(wide_60)
wide_60$time = ymd(wide_60$time)
wide_60$doy <- weekdays(wide_60$time)
wide_60[,c("gender","YZ","doy")]<-lapply(wide_60[,c("gender","YZ","doy")],factor)

wide_men_19 <- read_dta("D:/Stata17/ado/personal/Data/wide_men_from_19_to_45_years.dta")
wide_men_19 <- as.data.frame(wide_men_19)
wide_men_19$time = ymd(wide_men_19$time)
wide_men_19$doy <- weekdays(wide_men_19$time)
wide_men_19[,c("YZ","doy")]<-lapply(wide_men_19[,c("YZ","doy")],factor)

wide_women_19 <- read_dta("D:/Stata17/ado/personal/Data/wide_women_from_19_to_45_years.dta")
wide_women_19 <- as.data.frame(wide_women_19)
wide_women_19$time = ymd(wide_women_19$time)
wide_women_19$doy <- weekdays(wide_women_19$time)
wide_women_19[,c("YZ","doy")]<-lapply(wide_women_19[,c("YZ","doy")],factor)

wide_men_46 <- read_dta("D:/Stata17/ado/personal/Data/wide_men_from_46_to_59_years.dta")
wide_men_46 <- as.data.frame(wide_men_46)
wide_men_46$time = ymd(wide_men_46$time)
wide_men_46$doy <- weekdays(wide_men_46$time)
wide_men_46[,c("YZ","doy")]<-lapply(wide_men_46[,c("YZ","doy")],factor)

wide_women_46 <- read_dta("D:/Stata17/ado/personal/Data/wide_women_from_46_to_59_years.dta")
wide_women_46 <- as.data.frame(wide_women_46)
wide_women_46$time = ymd(wide_women_46$time)
wide_women_46$doy <- weekdays(wide_women_46$time)
wide_women_46[,c("YZ","doy")]<-lapply(wide_women_46[,c("YZ","doy")],factor)

wide_men_60 <- read_dta("D:/Stata17/ado/personal/Data/wide_men_over_60_years.dta")
wide_men_60 <- as.data.frame(wide_men_60)
wide_men_60$time = ymd(wide_men_60$time)
wide_men_60$doy <- weekdays(wide_men_60$time)
wide_men_60[,c("YZ","doy")]<-lapply(wide_men_60[,c("YZ","doy")],factor)

wide_women_60 <- read_dta("D:/Stata17/ado/personal/Data/wide_women_over_60_years.dta")
wide_women_60 <- as.data.frame(wide_women_60)
wide_women_60$time = ymd(wide_women_60$time)
wide_women_60$doy <- weekdays(wide_women_60$time)
wide_women_60[,c("YZ","doy")]<-lapply(wide_women_60[,c("YZ","doy")],factor)

#overall
#outdoor artificial
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")

#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")


#men
#outdoor artificial
alan_men <- wide_men[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                        "deng6","deng7","deng8","deng9","deng10","deng11",
                        "deng12","deng13","deng14","deng15","deng16","deng17",
                        "deng18","deng19","deng20","deng21","deng22","deng23",
                        "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men <- wide_men[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                        "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                        "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men) <- paste("lag", 0:30, sep="")

#PM10
pm10_men <- wide_men[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                        "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                        "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                        "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men <- wide_men[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                             "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                             "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                             "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                             "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men) <- paste("lag", 0:30, sep="")

#air pressure
ap_men <- wide_men[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                      "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                      "qiya12","qiya13","qiya14","qiya15","qiya16",
                      "qiya17","qiya18","qiya19","qiya20","qiya21",
                      "qiya22","qiya23","qiya24","qiya25",
                      "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men <- wide_men[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                      "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                      "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                      "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                      "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men) <- paste("lag", 0:30, sep="")


#women
#outdoor artificial
alan_women <- wide_women[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                            "deng6","deng7","deng8","deng9","deng10","deng11",
                            "deng12","deng13","deng14","deng15","deng16","deng17",
                            "deng18","deng19","deng20","deng21","deng22","deng23",
                            "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women <- wide_women[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                            "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                            "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women) <- paste("lag", 0:30, sep="")

#PM10
pm10_women <- wide_women[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                            "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                            "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                            "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women <- wide_women[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                 "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                 "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                 "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                 "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women) <- paste("lag", 0:30, sep="")

#air pressure
ap_women <- wide_women[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                          "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                          "qiya12","qiya13","qiya14","qiya15","qiya16",
                          "qiya17","qiya18","qiya19","qiya20","qiya21",
                          "qiya22","qiya23","qiya24","qiya25",
                          "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women <- wide_women[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                          "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                          "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                          "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                          "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women) <- paste("lag", 0:30, sep="")


#19~45 years old
#outdoor artificial
alan_19 <- wide_19[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                      "deng6","deng7","deng8","deng9","deng10","deng11",
                      "deng12","deng13","deng14","deng15","deng16","deng17",
                      "deng18","deng19","deng20","deng21","deng22","deng23",
                      "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_19) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_19 <- wide_19[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                      "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                      "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_19) <- paste("lag", 0:30, sep="")

#PM10
pm10_19 <- wide_19[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                      "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                      "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                      "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_19) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_19 <- wide_19[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                           "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                           "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                           "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                           "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_19) <- paste("lag", 0:30, sep="")

#air pressure
ap_19 <- wide_19[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                    "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                    "qiya12","qiya13","qiya14","qiya15","qiya16",
                    "qiya17","qiya18","qiya19","qiya20","qiya21",
                    "qiya22","qiya23","qiya24","qiya25",
                    "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_19) <- paste("lag", 0:30, sep="")

#relative humidity
rh_19 <- wide_19[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                    "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                    "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                    "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                    "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_19) <- paste("lag", 0:30, sep="")


#46~59 years old
#outdoor artificial
alan_46 <- wide_46[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                      "deng6","deng7","deng8","deng9","deng10","deng11",
                      "deng12","deng13","deng14","deng15","deng16","deng17",
                      "deng18","deng19","deng20","deng21","deng22","deng23",
                      "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_46) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_46 <- wide_46[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                      "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                      "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_46) <- paste("lag", 0:30, sep="")

#PM10
pm10_46 <- wide_46[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                      "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                      "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                      "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_46) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_46 <- wide_46[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                           "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                           "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                           "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                           "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_46) <- paste("lag", 0:30, sep="")

#air pressure
ap_46 <- wide_46[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                    "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                    "qiya12","qiya13","qiya14","qiya15","qiya16",
                    "qiya17","qiya18","qiya19","qiya20","qiya21",
                    "qiya22","qiya23","qiya24","qiya25",
                    "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_46) <- paste("lag", 0:30, sep="")

#relative humidity
rh_46 <- wide_46[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                    "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                    "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                    "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                    "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_46) <- paste("lag", 0:30, sep="")


# >60 years old
#outdoor artificial
alan_60 <- wide_60[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                      "deng6","deng7","deng8","deng9","deng10","deng11",
                      "deng12","deng13","deng14","deng15","deng16","deng17",
                      "deng18","deng19","deng20","deng21","deng22","deng23",
                      "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_60) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_60 <- wide_60[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                      "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                      "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_60) <- paste("lag", 0:30, sep="")

#PM10
pm10_60 <- wide_60[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                      "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                      "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                      "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_60) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_60 <- wide_60[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                           "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                           "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                           "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                           "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_60) <- paste("lag", 0:30, sep="")

#air pressure
ap_60 <- wide_60[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                    "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                    "qiya12","qiya13","qiya14","qiya15","qiya16",
                    "qiya17","qiya18","qiya19","qiya20","qiya21",
                    "qiya22","qiya23","qiya24","qiya25",
                    "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_60) <- paste("lag", 0:30, sep="")

#relative humidity
rh_60 <- wide_60[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                    "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                    "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                    "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                    "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_60) <- paste("lag", 0:30, sep="")


#men 19~45years
#outdoor artificial
alan_men_19 <- wide_men_19[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men_19) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men_19 <- wide_men_19[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men_19) <- paste("lag", 0:30, sep="")

#PM10
pm10_men_19 <- wide_men_19[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men_19) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men_19 <- wide_men_19[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men_19) <- paste("lag", 0:30, sep="")

#air pressure
ap_men_19 <- wide_men_19[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men_19) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men_19 <- wide_men_19[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men_19) <- paste("lag", 0:30, sep="")


#women 19~45years
#outdoor artificial
alan_women_19 <- wide_women_19[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                                  "deng6","deng7","deng8","deng9","deng10","deng11",
                                  "deng12","deng13","deng14","deng15","deng16","deng17",
                                  "deng18","deng19","deng20","deng21","deng22","deng23",
                                  "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women_19) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women_19 <- wide_women_19[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                                  "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                                  "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women_19) <- paste("lag", 0:30, sep="")

#PM10
pm10_women_19 <- wide_women_19[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                                  "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                                  "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                                  "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women_19) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women_19 <- wide_women_19[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                       "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                       "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                       "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                       "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women_19) <- paste("lag", 0:30, sep="")

#air pressure
ap_women_19 <- wide_women_19[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                                "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                                "qiya12","qiya13","qiya14","qiya15","qiya16",
                                "qiya17","qiya18","qiya19","qiya20","qiya21",
                                "qiya22","qiya23","qiya24","qiya25",
                                "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women_19) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women_19 <- wide_women_19[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                                "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                                "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                                "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                                "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women_19) <- paste("lag", 0:30, sep="")


#men 46~59years
#outdoor artificial
alan_men_46 <- wide_men_46[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men_46) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men_46 <- wide_men_46[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men_46) <- paste("lag", 0:30, sep="")

#PM10
pm10_men_46 <- wide_men_46[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men_46) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men_46 <- wide_men_46[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men_46) <- paste("lag", 0:30, sep="")

#air pressure
ap_men_46 <- wide_men_46[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men_46) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men_46 <- wide_men_46[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men_46) <- paste("lag", 0:30, sep="")


#women 46~59years
#outdoor artificial
alan_women_46 <- wide_women_46[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                                  "deng6","deng7","deng8","deng9","deng10","deng11",
                                  "deng12","deng13","deng14","deng15","deng16","deng17",
                                  "deng18","deng19","deng20","deng21","deng22","deng23",
                                  "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women_46) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women_46 <- wide_women_46[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                                  "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                                  "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women_46) <- paste("lag", 0:30, sep="")

#PM10
pm10_women_46 <- wide_women_46[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                                  "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                                  "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                                  "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women_46) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women_46 <- wide_women_46[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                       "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                       "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                       "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                       "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women_46) <- paste("lag", 0:30, sep="")

#air pressure
ap_women_46 <- wide_women_46[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                                "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                                "qiya12","qiya13","qiya14","qiya15","qiya16",
                                "qiya17","qiya18","qiya19","qiya20","qiya21",
                                "qiya22","qiya23","qiya24","qiya25",
                                "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women_46) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women_46 <- wide_women_46[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                                "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                                "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                                "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                                "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women_46) <- paste("lag", 0:30, sep="")


#men >60 years
#outdoor artificial
alan_men_60 <- wide_men_60[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men_60) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men_60 <- wide_men_60[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men_60) <- paste("lag", 0:30, sep="")

#PM10
pm10_men_60 <- wide_men_60[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men_60) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men_60 <- wide_men_60[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men_60) <- paste("lag", 0:30, sep="")

#air pressure
ap_men_60 <- wide_men_60[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men_60) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men_60 <- wide_men_60[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men_60) <- paste("lag", 0:30, sep="")


#women >60 years
#outdoor artificial
alan_women_60 <- wide_women_60[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                                  "deng6","deng7","deng8","deng9","deng10","deng11",
                                  "deng12","deng13","deng14","deng15","deng16","deng17",
                                  "deng18","deng19","deng20","deng21","deng22","deng23",
                                  "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women_60) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women_60 <- wide_women_60[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                                  "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                                  "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women_60) <- paste("lag", 0:30, sep="")

#PM10
pm10_women_60 <- wide_women_60[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                                  "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                                  "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                                  "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women_60) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women_60 <- wide_women_60[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                       "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                       "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                       "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                       "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women_60) <- paste("lag", 0:30, sep="")

#air pressure
ap_women_60 <- wide_women_60[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                                "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                                "qiya12","qiya13","qiya14","qiya15","qiya16",
                                "qiya17","qiya18","qiya19","qiya20","qiya21",
                                "qiya22","qiya23","qiya24","qiya25",
                                "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women_60) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women_60 <- wide_women_60[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                                "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                                "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                                "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                                "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women_60) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#overall
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#men
#outdoor ALAN
AL_men <- crossbasis(alan_men,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men <- crossbasis(pm25_men,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men <- crossbasis(pm10_men,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men,nk=2)
MeantTemp_men <- crossbasis(meanttemp_men,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men,nk=2)
AP_men <- crossbasis(ap_men,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men,nk=2)
RH_men <- crossbasis(rh_men,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#women
#outdoor ALAN
AL_women <- crossbasis(alan_women,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women <- crossbasis(pm25_women,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women <- crossbasis(pm10_women,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women,nk=2)
MeantTemp_women <- crossbasis(meanttemp_women,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women,nk=2)
AP_women <- crossbasis(ap_women,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women,nk=2)
RH_women <- crossbasis(rh_women,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#19~45 years
#outdoor ALAN
AL_19 <- crossbasis(alan_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_19 <- crossbasis(pm25_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_19 <- crossbasis(pm10_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_19,nk=2)
MeantTemp_19 <- crossbasis(meanttemp_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_19,nk=2)
AP_19 <- crossbasis(ap_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_19,nk=2)
RH_19 <- crossbasis(rh_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#46~59 years
#outdoor ALAN
AL_46 <- crossbasis(alan_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_46 <- crossbasis(pm25_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_46 <- crossbasis(pm10_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_46,nk=2)
MeantTemp_46 <- crossbasis(meanttemp_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_46,nk=2)
AP_46 <- crossbasis(ap_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_46,nk=2)
RH_46 <- crossbasis(rh_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#>60 years
#outdoor ALAN
AL_60 <- crossbasis(alan_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_60 <- crossbasis(pm25_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_60 <- crossbasis(pm10_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_60,nk=2)
MeantTemp_60 <- crossbasis(meanttemp_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_60,nk=2)
AP_60 <- crossbasis(ap_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_60,nk=2)
RH_60 <- crossbasis(rh_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#men 19~45
#outdoor ALAN
AL_men_19 <- crossbasis(alan_men_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men_19 <- crossbasis(pm25_men_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men_19 <- crossbasis(pm10_men_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men_19,nk=2)
MeantTemp_men_19 <- crossbasis(meanttemp_men_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men_19,nk=2)
AP_men_19 <- crossbasis(ap_men_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men_19,nk=2)
RH_men_19 <- crossbasis(rh_men_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#women 19~45
#outdoor ALAN
AL_women_19 <- crossbasis(alan_women_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women_19 <- crossbasis(pm25_women_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women_19 <- crossbasis(pm10_women_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women_19,nk=2)
MeantTemp_women_19 <- crossbasis(meanttemp_women_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women_19,nk=2)
AP_women_19 <- crossbasis(ap_women_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women_19,nk=2)
RH_women_19 <- crossbasis(rh_women_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#men 46~59
#outdoor ALAN
AL_men_46 <- crossbasis(alan_men_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men_46 <- crossbasis(pm25_men_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men_46 <- crossbasis(pm10_men_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men_46,nk=2)
MeantTemp_men_46 <- crossbasis(meanttemp_men_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men_46,nk=2)
AP_men_46 <- crossbasis(ap_men_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men_46,nk=2)
RH_men_46 <- crossbasis(rh_men_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#women 46~59
#outdoor ALAN
AL_women_46 <- crossbasis(alan_women_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women_46 <- crossbasis(pm25_women_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women_46 <- crossbasis(pm10_women_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women_46,nk=2)
MeantTemp_women_46 <- crossbasis(meanttemp_women_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women_46,nk=2)
AP_women_46 <- crossbasis(ap_women_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women_46,nk=2)
RH_women_46 <- crossbasis(rh_women_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#men >60
#outdoor ALAN
AL_men_60 <- crossbasis(alan_men_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men_60 <- crossbasis(pm25_men_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men_60 <- crossbasis(pm10_men_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men_60,nk=2)
MeantTemp_men_60 <- crossbasis(meanttemp_men_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men_60,nk=2)
AP_men_60 <- crossbasis(ap_men_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men_60,nk=2)
RH_men_60 <- crossbasis(rh_men_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#women >60
#outdoor ALAN
AL_women_60 <- crossbasis(alan_women_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women_60 <- crossbasis(pm25_women_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women_60 <- crossbasis(pm10_women_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women_60,nk=2)
MeantTemp_women_60 <- crossbasis(meanttemp_women_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women_60,nk=2)
AP_women_60 <- crossbasis(ap_women_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women_60,nk=2)
RH_women_60 <- crossbasis(rh_women_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+factor(doy),family=binomial,data=wide)

model_men <-glm(YZ ~ AL_men + PM25_men +PM10_men +AP_men +MeantTemp_men +RH_men +factor(doy),family=binomial,data=wide_men)

model_women <-glm(YZ ~ AL_women + PM25_women +PM10_women +AP_women +MeantTemp_women +RH_women +factor(doy),family=binomial,data=wide_women)

model_19 <-glm(YZ ~ AL_19 + PM25_19 +PM10_19 +AP_19 +MeantTemp_19 +RH_19 +factor(doy),family=binomial,data=wide_19)

model_46 <-glm(YZ ~ AL_46 + PM25_46 +PM10_46+AP_46+MeantTemp_46+RH_46+factor(doy),family=binomial,data=wide_46)

model_60 <-glm(YZ ~ AL_60 + PM25_60+PM10_60+AP_60+MeantTemp_60+RH_60+factor(doy),family=binomial,data=wide_60)

model_men_19 <-glm(YZ ~ AL_men_19 + PM25_men_19+PM10_men_19+AP_men_19+MeantTemp_men_19+RH_men_19+factor(doy),family=binomial,data=wide_men_19)

model_women_19 <-glm(YZ ~ AL_women_19 + PM25_women_19+PM10_women_19+AP_women_19+MeantTemp_women_19+RH_women_19+factor(doy),family=binomial,data=wide_women_19)

model_men_46 <-glm(YZ ~ AL_men_46 +PM25_men_46+PM10_men_46+AP_men_46+MeantTemp_men_46+RH_men_46+factor(doy),family=binomial,data=wide_men_46)

model_women_46 <-glm(YZ ~ AL_women_46 + PM25_women_46+PM10_women_46+AP_women_46+MeantTemp_women_46+RH_women_46+factor(doy),family=binomial,data=wide_women_46)

model_men_60 <-glm(YZ ~ AL_men_60 +PM25_men_60+PM10_men_60+AP_men_60+MeantTemp_men_60+RH_men_60+factor(doy),family=binomial,data=wide_men_60)

model_women_60 <-glm(YZ ~ AL_women_60 +PM25_women_60+PM10_women_60+AP_women_60+MeantTemp_women_60+RH_women_60+factor(doy),family=binomial,data=wide_women_60)
###############################predict##################################
pred <- crosspred(AL,model,cen=0,at=1:99,by=0.2,cumul = T)
pred_men <- crosspred(AL_men,model_men,cen=0,at=1:99,by=0.2,cumul = T)
pred_women <- crosspred(AL_women,model_women,cen=0,at=1:99,by=0.2,cumul = T)
pred_19 <- crosspred(AL_19,model_19,cen=0,at=1:99,by=0.2,cumul = T)
pred_46 <- crosspred(AL_46,model_46,cen=0,at=1:99,by=0.2,cumul = T)
pred_60 <- crosspred(AL_60,model_60,cen=0,at=1:99,by=0.2,cumul = T)
pred_men_19 <- crosspred(AL_men_19,model_men_19,cen=0,at=1:99,by=0.2,cumul = T)
pred_women_19 <- crosspred(AL_women_19,model_women_19,cen=0,at=1:99,by=0.2,cumul = T)
pred_men_46 <- crosspred(AL_men_46,model_men_46,cen=0,at=1:99,by=0.2,cumul = T)
pred_women_46 <- crosspred(AL_women_46,model_women_46,cen=0,at=1:99,by=0.2,cumul = T)
pred_men_60 <- crosspred(AL_men_60,model_men_60,cen=0,at=1:99,by=0.2,cumul = T)
pred_women_60 <- crosspred(AL_women_60,model_women_60,cen=0,at=1:99,by=0.2,cumul = T)

tiff(file = "Overall.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Men.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_men, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_men, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_men, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_men, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_men, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Women.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_women, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_women, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_women, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_women, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_women, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "19~45 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_19, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_19, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_19, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_19, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_19, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "46~59 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_46, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_46, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_46, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_46, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_46, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "≥60 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_60, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_60, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_60, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_60, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_60, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Men in 19~45 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_men_19, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_men_19, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_men_19, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_men_19, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_men_19, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Women in 19~45 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_women_19, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_women_19, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_women_19, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_women_19, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_women_19, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Men in 46~59 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_men_46, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_men_46, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_men_46, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_men_46, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_men_46, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Women in 46~59 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_women_46, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_women_46, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_women_46, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_women_46, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_women_46, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Men of ≥60 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_men_60, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_men_60, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_men_60, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_men_60, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_men_60, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()

tiff(file = "Women of ≥60 years.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))
plot(pred_women_60, "slices", var = 1,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P2.5",cex.main=1.5,cumul = T)
plot(pred_women_60, "slices", var = 11,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P25",cex.main=1.5,cumul = T)
plot(pred_women_60, "slices", var = 26,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P50",cex.main=1.5,cumul = T)
plot(pred_women_60, "slices", var = 38,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P75",cex.main=1.5,cumul = T)
plot(pred_women_60, "slices", var = 66,col=2,cex.axis=1.50,xlab="lag(days)",ylab="OR",main="P97.5",cex.main=1.5,cumul = T)
dev.off()


##############################################################################
##############################################################################
#####################Cumulative association#######################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)

wide_men <- read_dta("D:/Stata17/ado/personal/Data/wide_men.dta")
wide_men <- as.data.frame(wide_men)
wide_men$time = ymd(wide_men$time)
wide_men$doy <- weekdays(wide_men$time)
wide_men[,c("YZ","doy")]<-lapply(wide_men[,c("YZ","doy")],factor)

wide_women <- read_dta("D:/Stata17/ado/personal/Data/wide_women.dta")
wide_women <- as.data.frame(wide_women)
wide_women$time = ymd(wide_women$time)
wide_women$doy <- weekdays(wide_women$time)
wide_women[,c("YZ","doy")]<-lapply(wide_women[,c("YZ","doy")],factor)

wide_19 <- read_dta("D:/Stata17/ado/personal/Data/wide_from_19_to_45_years.dta")
wide_19 <- as.data.frame(wide_19)
wide_19$time = ymd(wide_19$time)
wide_19$doy <- weekdays(wide_19$time)
wide_19[,c("gender","YZ","doy")]<-lapply(wide_19[,c("gender","YZ","doy")],factor)

wide_46 <- read_dta("D:/Stata17/ado/personal/Data/wide_from_46_to_59_years.dta")
wide_46 <- as.data.frame(wide_46)
wide_46$time = ymd(wide_46$time)
wide_46$doy <- weekdays(wide_46$time)
wide_46[,c("gender","YZ","doy")]<-lapply(wide_46[,c("gender","YZ","doy")],factor)

wide_60 <- read_dta("D:/Stata17/ado/personal/Data/wide_over_60_years.dta")
wide_60 <- as.data.frame(wide_60)
wide_60$time = ymd(wide_60$time)
wide_60$doy <- weekdays(wide_60$time)
wide_60[,c("gender","YZ","doy")]<-lapply(wide_60[,c("gender","YZ","doy")],factor)

wide_men_19 <- read_dta("D:/Stata17/ado/personal/Data/wide_men_from_19_to_45_years.dta")
wide_men_19 <- as.data.frame(wide_men_19)
wide_men_19$time = ymd(wide_men_19$time)
wide_men_19$doy <- weekdays(wide_men_19$time)
wide_men_19[,c("YZ","doy")]<-lapply(wide_men_19[,c("YZ","doy")],factor)

wide_women_19 <- read_dta("D:/Stata17/ado/personal/Data/wide_women_from_19_to_45_years.dta")
wide_women_19 <- as.data.frame(wide_women_19)
wide_women_19$time = ymd(wide_women_19$time)
wide_women_19$doy <- weekdays(wide_women_19$time)
wide_women_19[,c("YZ","doy")]<-lapply(wide_women_19[,c("YZ","doy")],factor)

wide_men_46 <- read_dta("D:/Stata17/ado/personal/Data/wide_men_from_46_to_59_years.dta")
wide_men_46 <- as.data.frame(wide_men_46)
wide_men_46$time = ymd(wide_men_46$time)
wide_men_46$doy <- weekdays(wide_men_46$time)
wide_men_46[,c("YZ","doy")]<-lapply(wide_men_46[,c("YZ","doy")],factor)

wide_women_46 <- read_dta("D:/Stata17/ado/personal/Data/wide_women_from_46_to_59_years.dta")
wide_women_46 <- as.data.frame(wide_women_46)
wide_women_46$time = ymd(wide_women_46$time)
wide_women_46$doy <- weekdays(wide_women_46$time)
wide_women_46[,c("YZ","doy")]<-lapply(wide_women_46[,c("YZ","doy")],factor)

wide_men_60 <- read_dta("D:/Stata17/ado/personal/Data/wide_men_over_60_years.dta")
wide_men_60 <- as.data.frame(wide_men_60)
wide_men_60$time = ymd(wide_men_60$time)
wide_men_60$doy <- weekdays(wide_men_60$time)
wide_men_60[,c("YZ","doy")]<-lapply(wide_men_60[,c("YZ","doy")],factor)

wide_women_60 <- read_dta("D:/Stata17/ado/personal/Data/wide_women_over_60_years.dta")
wide_women_60 <- as.data.frame(wide_women_60)
wide_women_60$time = ymd(wide_women_60$time)
wide_women_60$doy <- weekdays(wide_women_60$time)
wide_women_60[,c("YZ","doy")]<-lapply(wide_women_60[,c("YZ","doy")],factor)

#overall
#outdoor artificial
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")

#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")


#men
#outdoor artificial
alan_men <- wide_men[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men <- wide_men[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men) <- paste("lag", 0:30, sep="")

#PM10
pm10_men <- wide_men[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men <- wide_men[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men) <- paste("lag", 0:30, sep="")

#air pressure
ap_men <- wide_men[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men <- wide_men[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men) <- paste("lag", 0:30, sep="")


#women
#outdoor artificial
alan_women <- wide_women[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                        "deng6","deng7","deng8","deng9","deng10","deng11",
                        "deng12","deng13","deng14","deng15","deng16","deng17",
                        "deng18","deng19","deng20","deng21","deng22","deng23",
                        "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women <- wide_women[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                        "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                        "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women) <- paste("lag", 0:30, sep="")

#PM10
pm10_women <- wide_women[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                        "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                        "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                        "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women <- wide_women[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                             "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                             "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                             "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                             "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women) <- paste("lag", 0:30, sep="")

#air pressure
ap_women <- wide_women[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                      "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                      "qiya12","qiya13","qiya14","qiya15","qiya16",
                      "qiya17","qiya18","qiya19","qiya20","qiya21",
                      "qiya22","qiya23","qiya24","qiya25",
                      "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women <- wide_women[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                      "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                      "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                      "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                      "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women) <- paste("lag", 0:30, sep="")


#19~45 years old
#outdoor artificial
alan_19 <- wide_19[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                            "deng6","deng7","deng8","deng9","deng10","deng11",
                            "deng12","deng13","deng14","deng15","deng16","deng17",
                            "deng18","deng19","deng20","deng21","deng22","deng23",
                            "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_19) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_19 <- wide_19[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                            "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                            "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_19) <- paste("lag", 0:30, sep="")

#PM10
pm10_19 <- wide_19[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                            "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                            "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                            "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_19) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_19 <- wide_19[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                 "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                 "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                 "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                 "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_19) <- paste("lag", 0:30, sep="")

#air pressure
ap_19 <- wide_19[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                          "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                          "qiya12","qiya13","qiya14","qiya15","qiya16",
                          "qiya17","qiya18","qiya19","qiya20","qiya21",
                          "qiya22","qiya23","qiya24","qiya25",
                          "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_19) <- paste("lag", 0:30, sep="")

#relative humidity
rh_19 <- wide_19[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                          "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                          "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                          "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                          "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_19) <- paste("lag", 0:30, sep="")


#46~59 years old
#outdoor artificial
alan_46 <- wide_46[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                      "deng6","deng7","deng8","deng9","deng10","deng11",
                      "deng12","deng13","deng14","deng15","deng16","deng17",
                      "deng18","deng19","deng20","deng21","deng22","deng23",
                      "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_46) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_46 <- wide_46[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                      "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                      "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_46) <- paste("lag", 0:30, sep="")

#PM10
pm10_46 <- wide_46[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                      "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                      "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                      "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_46) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_46 <- wide_46[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                           "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                           "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                           "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                           "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_46) <- paste("lag", 0:30, sep="")

#air pressure
ap_46 <- wide_46[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                    "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                    "qiya12","qiya13","qiya14","qiya15","qiya16",
                    "qiya17","qiya18","qiya19","qiya20","qiya21",
                    "qiya22","qiya23","qiya24","qiya25",
                    "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_46) <- paste("lag", 0:30, sep="")

#relative humidity
rh_46 <- wide_46[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                    "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                    "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                    "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                    "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_46) <- paste("lag", 0:30, sep="")


# >60 years old
#outdoor artificial
alan_60 <- wide_60[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                      "deng6","deng7","deng8","deng9","deng10","deng11",
                      "deng12","deng13","deng14","deng15","deng16","deng17",
                      "deng18","deng19","deng20","deng21","deng22","deng23",
                      "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_60) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_60 <- wide_60[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                      "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                      "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_60) <- paste("lag", 0:30, sep="")

#PM10
pm10_60 <- wide_60[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                      "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                      "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                      "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_60) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_60 <- wide_60[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                           "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                           "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                           "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                           "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_60) <- paste("lag", 0:30, sep="")

#air pressure
ap_60 <- wide_60[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                    "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                    "qiya12","qiya13","qiya14","qiya15","qiya16",
                    "qiya17","qiya18","qiya19","qiya20","qiya21",
                    "qiya22","qiya23","qiya24","qiya25",
                    "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_60) <- paste("lag", 0:30, sep="")

#relative humidity
rh_60 <- wide_60[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                    "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                    "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                    "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                    "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_60) <- paste("lag", 0:30, sep="")


#men 19~45years
#outdoor artificial
alan_men_19 <- wide_men_19[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                        "deng6","deng7","deng8","deng9","deng10","deng11",
                        "deng12","deng13","deng14","deng15","deng16","deng17",
                        "deng18","deng19","deng20","deng21","deng22","deng23",
                        "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men_19) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men_19 <- wide_men_19[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                        "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                        "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men_19) <- paste("lag", 0:30, sep="")

#PM10
pm10_men_19 <- wide_men_19[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                        "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                        "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                        "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men_19) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men_19 <- wide_men_19[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                             "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                             "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                             "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                             "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men_19) <- paste("lag", 0:30, sep="")

#air pressure
ap_men_19 <- wide_men_19[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                      "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                      "qiya12","qiya13","qiya14","qiya15","qiya16",
                      "qiya17","qiya18","qiya19","qiya20","qiya21",
                      "qiya22","qiya23","qiya24","qiya25",
                      "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men_19) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men_19 <- wide_men_19[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                      "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                      "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                      "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                      "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men_19) <- paste("lag", 0:30, sep="")


#women 19~45years
#outdoor artificial
alan_women_19 <- wide_women_19[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women_19) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women_19 <- wide_women_19[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women_19) <- paste("lag", 0:30, sep="")

#PM10
pm10_women_19 <- wide_women_19[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women_19) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women_19 <- wide_women_19[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women_19) <- paste("lag", 0:30, sep="")

#air pressure
ap_women_19 <- wide_women_19[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women_19) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women_19 <- wide_women_19[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women_19) <- paste("lag", 0:30, sep="")


#men 46~59years
#outdoor artificial
alan_men_46 <- wide_men_46[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men_46) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men_46 <- wide_men_46[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men_46) <- paste("lag", 0:30, sep="")

#PM10
pm10_men_46 <- wide_men_46[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men_46) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men_46 <- wide_men_46[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men_46) <- paste("lag", 0:30, sep="")

#air pressure
ap_men_46 <- wide_men_46[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men_46) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men_46 <- wide_men_46[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men_46) <- paste("lag", 0:30, sep="")


#women 46~59years
#outdoor artificial
alan_women_46 <- wide_women_46[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women_46) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women_46 <- wide_women_46[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women_46) <- paste("lag", 0:30, sep="")

#PM10
pm10_women_46 <- wide_women_46[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women_46) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women_46 <- wide_women_46[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women_46) <- paste("lag", 0:30, sep="")

#air pressure
ap_women_46 <- wide_women_46[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women_46) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women_46 <- wide_women_46[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women_46) <- paste("lag", 0:30, sep="")


#men >60 years
#outdoor artificial
alan_men_60 <- wide_men_60[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_men_60) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_men_60 <- wide_men_60[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_men_60) <- paste("lag", 0:30, sep="")

#PM10
pm10_men_60 <- wide_men_60[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_men_60) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_men_60 <- wide_men_60[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_men_60) <- paste("lag", 0:30, sep="")

#air pressure
ap_men_60 <- wide_men_60[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_men_60) <- paste("lag", 0:30, sep="")

#relative humidity
rh_men_60 <- wide_men_60[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_men_60) <- paste("lag", 0:30, sep="")


#women >60 years
#outdoor artificial
alan_women_60 <- wide_women_60[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                              "deng6","deng7","deng8","deng9","deng10","deng11",
                              "deng12","deng13","deng14","deng15","deng16","deng17",
                              "deng18","deng19","deng20","deng21","deng22","deng23",
                              "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan_women_60) <- paste("lag", 0:30, sep="")

#PM2.5
pm25_women_60 <- wide_women_60[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                              "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                              "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25_women_60) <- paste("lag", 0:30, sep="")

#PM10
pm10_women_60 <- wide_women_60[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                              "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                              "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                              "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10_women_60) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp_women_60 <- wide_women_60[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                                   "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                                   "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                                   "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                                   "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp_women_60) <- paste("lag", 0:30, sep="")

#air pressure
ap_women_60 <- wide_women_60[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
                            "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
                            "qiya12","qiya13","qiya14","qiya15","qiya16",
                            "qiya17","qiya18","qiya19","qiya20","qiya21",
                            "qiya22","qiya23","qiya24","qiya25",
                            "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap_women_60) <- paste("lag", 0:30, sep="")

#relative humidity
rh_women_60 <- wide_women_60[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
                            "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
                            "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
                            "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
                            "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh_women_60) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#overall
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#men
#outdoor ALAN
AL_men <- crossbasis(alan_men,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men <- crossbasis(pm25_men,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men <- crossbasis(pm10_men,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men,nk=2)
MeantTemp_men <- crossbasis(meanttemp_men,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men,nk=2)
AP_men <- crossbasis(ap_men,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men,nk=2)
RH_men <- crossbasis(rh_men,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#women
#outdoor ALAN
AL_women <- crossbasis(alan_women,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women <- crossbasis(pm25_women,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women <- crossbasis(pm10_women,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women,nk=2)
MeantTemp_women <- crossbasis(meanttemp_women,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women,nk=2)
AP_women <- crossbasis(ap_women,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women,nk=2)
RH_women <- crossbasis(rh_women,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#19~45 years
#outdoor ALAN
AL_19 <- crossbasis(alan_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_19 <- crossbasis(pm25_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_19 <- crossbasis(pm10_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_19,nk=2)
MeantTemp_19 <- crossbasis(meanttemp_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_19,nk=2)
AP_19 <- crossbasis(ap_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_19,nk=2)
RH_19 <- crossbasis(rh_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#46~59 years
#outdoor ALAN
AL_46 <- crossbasis(alan_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_46 <- crossbasis(pm25_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_46 <- crossbasis(pm10_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_46,nk=2)
MeantTemp_46 <- crossbasis(meanttemp_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_46,nk=2)
AP_46 <- crossbasis(ap_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_46,nk=2)
RH_46 <- crossbasis(rh_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#>60 years
#outdoor ALAN
AL_60 <- crossbasis(alan_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_60 <- crossbasis(pm25_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_60 <- crossbasis(pm10_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_60,nk=2)
MeantTemp_60 <- crossbasis(meanttemp_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_60,nk=2)
AP_60 <- crossbasis(ap_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_60,nk=2)
RH_60 <- crossbasis(rh_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#men 19~45
#outdoor ALAN
AL_men_19 <- crossbasis(alan_men_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men_19 <- crossbasis(pm25_men_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men_19 <- crossbasis(pm10_men_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men_19,nk=2)
MeantTemp_men_19 <- crossbasis(meanttemp_men_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men_19,nk=2)
AP_men_19 <- crossbasis(ap_men_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men_19,nk=2)
RH_men_19 <- crossbasis(rh_men_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#women 19~45
#outdoor ALAN
AL_women_19 <- crossbasis(alan_women_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women_19 <- crossbasis(pm25_women_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women_19 <- crossbasis(pm10_women_19,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women_19,nk=2)
MeantTemp_women_19 <- crossbasis(meanttemp_women_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women_19,nk=2)
AP_women_19 <- crossbasis(ap_women_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women_19,nk=2)
RH_women_19 <- crossbasis(rh_women_19,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#men 46~59
#outdoor ALAN
AL_men_46 <- crossbasis(alan_men_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men_46 <- crossbasis(pm25_men_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men_46 <- crossbasis(pm10_men_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men_46,nk=2)
MeantTemp_men_46 <- crossbasis(meanttemp_men_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men_46,nk=2)
AP_men_46 <- crossbasis(ap_men_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men_46,nk=2)
RH_men_46 <- crossbasis(rh_men_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#women 46~59
#outdoor ALAN
AL_women_46 <- crossbasis(alan_women_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women_46 <- crossbasis(pm25_women_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women_46 <- crossbasis(pm10_women_46,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women_46,nk=2)
MeantTemp_women_46 <- crossbasis(meanttemp_women_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women_46,nk=2)
AP_women_46 <- crossbasis(ap_women_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women_46,nk=2)
RH_women_46 <- crossbasis(rh_women_46,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

#men >60
#outdoor ALAN
AL_men_60 <- crossbasis(alan_men_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_men_60 <- crossbasis(pm25_men_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_men_60 <- crossbasis(pm10_men_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_men_60,nk=2)
MeantTemp_men_60 <- crossbasis(meanttemp_men_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_men_60,nk=2)
AP_men_60 <- crossbasis(ap_men_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_men_60,nk=2)
RH_men_60 <- crossbasis(rh_men_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))


#women >60
#outdoor ALAN
AL_women_60 <- crossbasis(alan_women_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25_women_60 <- crossbasis(pm25_women_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10_women_60 <- crossbasis(pm10_women_60,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp_women_60,nk=2)
MeantTemp_women_60 <- crossbasis(meanttemp_women_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap_women_60,nk=2)
AP_women_60 <- crossbasis(ap_women_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh_women_60,nk=2)
RH_women_60 <- crossbasis(rh_women_60,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+factor(doy),family=binomial,data=wide)

model_men <-glm(YZ ~ AL_men + PM25_men +PM10_men +AP_men +MeantTemp_men +RH_men +factor(doy),family=binomial,data=wide_men)

model_women <-glm(YZ ~ AL_women + PM25_women +PM10_women +AP_women +MeantTemp_women +RH_women +factor(doy),family=binomial,data=wide_women)

model_19 <-glm(YZ ~ AL_19 + PM25_19 +PM10_19 +AP_19 +MeantTemp_19 +RH_19 +factor(doy),family=binomial,data=wide_19)

model_46 <-glm(YZ ~ AL_46 + PM25_46 +PM10_46+AP_46+MeantTemp_46+RH_46+factor(doy),family=binomial,data=wide_46)

model_60 <-glm(YZ ~ AL_60 + PM25_60+PM10_60+AP_60+MeantTemp_60+RH_60+factor(doy),family=binomial,data=wide_60)

model_men_19 <-glm(YZ ~ AL_men_19 + PM25_men_19+PM10_men_19+AP_men_19+MeantTemp_men_19+RH_men_19+factor(doy),family=binomial,data=wide_men_19)

model_women_19 <-glm(YZ ~ AL_women_19 + PM25_women_19+PM10_women_19+AP_women_19+MeantTemp_women_19+RH_women_19+factor(doy),family=binomial,data=wide_women_19)

model_men_46 <-glm(YZ ~ AL_men_46 +PM25_men_46+PM10_men_46+AP_men_46+MeantTemp_men_46+RH_men_46+factor(doy),family=binomial,data=wide_men_46)

model_women_46 <-glm(YZ ~ AL_women_46 + PM25_women_46+PM10_women_46+AP_women_46+MeantTemp_women_46+RH_women_46+factor(doy),family=binomial,data=wide_women_46)

model_men_60 <-glm(YZ ~ AL_men_60 +PM25_men_60+PM10_men_60+AP_men_60+MeantTemp_men_60+RH_men_60+factor(doy),family=binomial,data=wide_men_60)

model_women_60 <-glm(YZ ~ AL_women_60 +PM25_women_60+PM10_women_60+AP_women_60+MeantTemp_women_60+RH_women_60+factor(doy),family=binomial,data=wide_women_60)
###############################predict##################################
cralla <- crossreduce(AL,model,cen=0)
crallb <- crossreduce(AL_men,model_men,cen=0)
crallc <- crossreduce(AL_women,model_women,cen=0)
cralld <- crossreduce(AL_19,model_19,cen=0)
cralle <- crossreduce(AL_46,model_46,cen=0)
crallf <- crossreduce(AL_60,model_60,cen=0)
crallg <- crossreduce(AL_men_19,model_men_19,cen=0)
crallh <- crossreduce(AL_women_19,model_women_19,cen=0)
cralli <- crossreduce(AL_men_46,model_men_46,cen=0)
crallj <- crossreduce(AL_women_46,model_women_46,cen=0)
crallk <- crossreduce(AL_men_60,model_men_60,cen=0)
cralll <- crossreduce(AL_women_60,model_women_60,cen=0)

###############################Overall图##################################
tiff(file = "Figure 2.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(2,3))

plot(cralla,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Overall",cex=1.0)

plot(crallb,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Men",cex=1.0)

plot(crallc,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Women",cex=1.0)

plot(cralli,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Men in 46~59 years",cex=1.0)

plot(crallj,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Women in 46~59 years",cex=1.0)
dev.off()


tiff(file = "Figure 2.tiff", width = 3000, height = 2500, res = 300)
par(mfrow=c(3,3))

plot(cralld,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="19~45 years",cex=1.0)

plot(cralle,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="46~59 years",cex=1.0)

plot(crallf,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="≥60 years",cex=1.0)

plot(crallg,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Men in 19~45 years",cex=1.0)

plot(crallh,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Women in 19~45 years",cex=1.0)

plot(crallk,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Men of ≥60 years",cex=1.0)

plot(cralll,xlab="outdoor ALAN",ylab="OR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Women of ≥60 years",cex=1.0)
dev.off()


##############################################################################
##############################################################################
###########Attributable fractions##################################################

########################attrdl##################################
### (c) Antonio Gasparrini 2015-2017###############################################
################################################################################
attrdl <- function(x,basis,cases,model=NULL,coef=NULL,vcov=NULL,model.link=NULL,
                   type="af",dir="back",tot=TRUE,cen,range=NULL,sim=FALSE,nsim=5000) {
  ################################################################################
  #
  # CHECK VERSION OF THE DLNM PACKAGE
  if(packageVersion("dlnm")<"2.2.0") 
    stop("update dlnm package to version >= 2.2.0")
  #
  # EXTRACT NAME AND CHECK type AND dir
  name <- deparse(substitute(basis))
  type <- match.arg(type,c("an","af"))
  dir <- match.arg(dir,c("back","forw"))
  #
  # DEFINE CENTERING
  if(missing(cen) && is.null(cen <- attr(basis,"argvar")$cen))
    stop("'cen' must be provided")
  if(!is.numeric(cen) && length(cen)>1L) stop("'cen' must be a numeric scalar")
  attributes(basis)$argvar$cen <- NULL
  #  
  # SELECT RANGE (FORCE TO CENTERING VALUE OTHERWISE, MEANING NULL RISK)
  if(!is.null(range)) x[x<range[1]|x>range[2]] <- cen
  #
  # COMPUTE THE MATRIX OF
  #   - LAGGED EXPOSURES IF dir="back"
  #   - CONSTANT EXPOSURES ALONG LAGS IF dir="forw"
  lag <- attr(basis,"lag")
  if(NCOL(x)==1L) {
    at <- if(dir=="back") tsModel:::Lag(x,seq(lag[1],lag[2])) else 
      matrix(rep(x,diff(lag)+1),length(x))
  } else {
    if(dir=="forw") stop("'x' must be a vector when dir='forw'")
    if(ncol(at <- x)!=diff(lag)+1) 
      stop("dimension of 'x' not compatible with 'basis'")
  }
  #
  # NUMBER USED FOR THE CONTRIBUTION AT EACH TIME IN FORWARD TYPE
  #   - IF cases PROVIDED AS A MATRIX, TAKE THE ROW AVERAGE
  #   - IF PROVIDED AS A TIME SERIES, COMPUTE THE FORWARD MOVING AVERAGE
  #   - THIS EXCLUDES MISSING ACCORDINGLY
  # ALSO COMPUTE THE DENOMINATOR TO BE USED BELOW
  if(NROW(cases)!=NROW(at)) stop("'x' and 'cases' not consistent")
  if(NCOL(cases)>1L) {
    if(dir=="back") stop("'cases' must be a vector if dir='back'")
    if(ncol(cases)!=diff(lag)+1) stop("dimension of 'cases' not compatible")
    den <- sum(rowMeans(cases,na.rm=TRUE),na.rm=TRUE)
    cases <- rowMeans(cases)
  } else {
    den <- sum(cases,na.rm=TRUE) 
    if(dir=="forw") 
      cases <- rowMeans(as.matrix(tsModel:::Lag(cases,-seq(lag[1],lag[2]))))
  }
  #
  ################################################################################
  #
  # EXTRACT COEF AND VCOV IF MODEL IS PROVIDED
  if(!is.null(model)) {
    cond <- paste0(name,"[[:print:]]*v[0-9]{1,2}\\.l[0-9]{1,2}")
    if(ncol(basis)==1L) cond <- name
    model.class <- class(model)
    coef <- dlnm:::getcoef(model,model.class)
    ind <- grep(cond,names(coef))
    coef <- coef[ind]
    vcov <- dlnm:::getvcov(model,model.class)[ind,ind,drop=FALSE]
    model.link <- dlnm:::getlink(model,model.class)
    if(!model.link %in% c("log","logit"))
      stop("'model' must have a log or logit link function")
  }
  #
  # IF REDUCED ESTIMATES ARE PROVIDED
  typebasis <- ifelse(length(coef)!=ncol(basis),"one","cb")
  #
  ################################################################################
  #
  # PREPARE THE ARGUMENTS FOR TH BASIS TRANSFORMATION
  predvar <- if(typebasis=="one") x else seq(NROW(at))
  predlag <- if(typebasis=="one") 0 else dlnm:::seqlag(lag)
  #  
  # CREATE THE MATRIX OF TRANSFORMED CENTRED VARIABLES (DEPENDENT ON typebasis)
  if(typebasis=="cb") {
    Xpred <- dlnm:::mkXpred(typebasis,basis,at,predvar,predlag,cen)
    Xpredall <- 0
    for (i in seq(length(predlag))) {
      ind <- seq(length(predvar))+length(predvar)*(i-1)
      Xpredall <- Xpredall + Xpred[ind,,drop=FALSE]
    }
  } else {
    basis <- do.call(onebasis,c(list(x=x),attr(basis,"argvar")))
    Xpredall <- dlnm:::mkXpred(typebasis,basis,x,predvar,predlag,cen)
  }
  #  
  # CHECK DIMENSIONS  
  if(length(coef)!=ncol(Xpredall))
    stop("arguments 'basis' do not match 'model' or 'coef'-'vcov'")
  if(any(dim(vcov)!=c(length(coef),length(coef)))) 
    stop("arguments 'coef' and 'vcov' do no match")
  if(typebasis=="one" && dir=="back")
    stop("only dir='forw' allowed for reduced estimates")
  #
  ################################################################################
  #
  # COMPUTE AF AND AN 
  af <- 1-exp(-drop(as.matrix(Xpredall%*%coef)))
  an <- af*cases
  #
  # TOTAL
  #   - SELECT NON-MISSING OBS CONTRIBUTING TO COMPUTATION
  #   - DERIVE TOTAL AF
  #   - COMPUTE TOTAL AN WITH ADJUSTED DENOMINATOR (OBSERVED TOTAL NUMBER)
  if(tot) {
    isna <- is.na(an)
    af <- sum(an[!isna])/sum(cases[!isna])
    an <- af*den
  }
  #
  ################################################################################
  #
  # EMPIRICAL CONFIDENCE INTERVALS
  if(!tot && sim) {
    sim <- FALSE
    warning("simulation samples only returned for tot=T")
  }
  if(sim) {
    # SAMPLE COEF
    k <- length(coef)
    eigen <- eigen(vcov)
    X <- matrix(rnorm(length(coef)*nsim),nsim)
    coefsim <- coef + eigen$vectors %*% diag(sqrt(eigen$values),k) %*% t(X)
    # RUN THE LOOP
    # pre_afsim <- (1 - exp(- Xpredall %*% coefsim)) * cases # a matrix
    # afsim <- colSums(pre_afsim,na.rm=TRUE) / sum(cases[!isna],na.rm=TRUE)
    afsim <- apply(coefsim,2, function(coefi) {
      ani <- (1-exp(-drop(Xpredall%*%coefi)))*cases
      sum(ani[!is.na(ani)])/sum(cases[!is.na(ani)])
    })
    ansim <- afsim*den
  }
  #
  ################################################################################
  #
  res <- if(sim) {
    if(type=="an") ansim else afsim
  } else {
    if(type=="an") an else af    
  }
  #
  return(res)
}

############################overall#################################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

################################Men################################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_men.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))


################################Women################################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_women.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################19~45 years#############################################
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_from_19_to_45_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################46~59 years#############################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_from_46_to_59_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################≥60 years#############################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_over_60_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################Men in 19~45 years#########################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_men_from_19_to_45_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################Women in 19~45 years#####################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_women_from_19_to_45_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################Men in 46~59 years#####################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_men_from_46_to_59_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################Women in 46~59 years#####################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_women_from_46_to_59_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################Men of ≥60 years#####################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_men_over_60_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#############################Women of ≥60 years#####################################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide_women_over_60_years.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
#################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")
alan <- as.matrix(alan)
#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

###########################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################

model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+ +factor(doy),family=binomial,data=wide)

Y <- wide[,7]
Y <- as.numeric(Y)
str(Y)
Y <- as.logical(Y)

pred <- crosspred(AL,model,cen=0,at=20,by=0.2)
with(pred,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=0)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=0),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=10,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=10)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=10),c(0.025,0.975))

pred1 <- crosspred(AL,model,cen=20,at=60,by=0.2)
with(pred1,cbind(allRRfit,allRRlow,allRRhigh))
attrdl(alan,AL,Y,model,cen=20)
quantile(attrdl(alan,AL,Y,model,sim=T,nsim=1000,cen=20),c(0.025,0.975))

#################################################################
#################################################################
####################Sensitivity analysis##############################
library(haven)
library(dlnm)
library(lubridate)
library(splines)
wide <- read_dta("D:/Stata17/ado/personal/Data/wide.dta")
wide <- as.data.frame(wide)
wide$time = ymd(wide$time)
wide$doy <- weekdays(wide$time)
wide[,c("gender","YZ","doy")]<-lapply(wide[,c("gender","YZ","doy")],factor)
################################################################
#outdoor ALAN
alan <- wide[,c("deng0","deng1","deng2","deng3","deng4","deng5",
                "deng6","deng7","deng8","deng9","deng10","deng11",
                "deng12","deng13","deng14","deng15","deng16","deng17",
                "deng18","deng19","deng20","deng21","deng22","deng23",
                "deng24","deng25","deng26","deng27","deng28","deng29","deng30")]
colnames(alan) <- paste("lag", 0:30, sep="")

#PM2.5
pm25 <- wide[,c("m0","m1","m2","m3","m4","m5","m6","m7","m8","m9","m10",
                "m11","m12","m13","m14","m15","m16","m17","m18","m19","m20",
                "m21","m22","m23","m24","m25","m26","m27","m28","m29","m30")]
colnames(pm25) <- paste("lag", 0:30, sep="")

#PM10
pm10 <- wide[,c("pm0","pm1","pm2","pm3","pm4","pm5","pm6","pm7","pm8",
                "pm9","pm10","pm11","pm12","pm13","pm14","pm15","pm16",
                "pm17","pm18","pm19","pm20","pm21","pm22","pm23","pm24",
                "pm25","pm26","pm27","pm28","pm29","pm30")]
colnames(pm10) <- paste("lag", 0:30, sep="")

#mean temperature
meanttemp <- wide[,c("qiwen0","qiwen1","qiwen2","qiwen3","qiwen4","qiwen5",
                     "qiwen6","qiwen7","qiwen8","qiwen9","qiwen10","qiwen11",
                     "qiwen12","qiwen13","qiwen14","qiwen15","qiwen16","qiwen17",
                     "qiwen18","qiwen19","qiwen20","qiwen21","qiwen22","qiwen23",
                     "qiwen24","qiwen25","qiwen26","qiwen27","qiwen28","qiwen29","qiwen30")]
colnames(meanttemp) <- paste("lag", 0:30, sep="")

#air pressure
ap <- wide[,c("qiya0","qiya1","qiya2","qiya3","qiya4","qiya5",
              "qiya6","qiya7","qiya8","qiya9","qiya10","qiya11",
              "qiya12","qiya13","qiya14","qiya15","qiya16",
              "qiya17","qiya18","qiya19","qiya20","qiya21",
              "qiya22","qiya23","qiya24","qiya25",
              "qiya26","qiya27","qiya28","qiya29","qiya30")]
colnames(ap) <- paste("lag", 0:30, sep="")

#relative humidity
rh <- wide[,c("shidu0","shidu1","shidu2","shidu3","shidu4","shidu5",
              "shidu6","shidu7","shidu8","shidu9","shidu10","shidu11",
              "shidu12","shidu13","shidu14","shidu15","shidu16","shidu17",
              "shidu18","shidu19","shidu20","shidu21","shidu22","shidu23",
              "shidu24","shidu25","shidu26","shidu27","shidu28","shidu29","shidu30")]
colnames(rh) <- paste("lag", 0:30, sep="")

############################crossbasis#######################################
#outdoor ALAN
AL <- crossbasis(alan,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM2.5
PM25 <- crossbasis(pm25,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#PM10
PM10 <- crossbasis(pm10,lag=c(0,30),argvar=list(fun="lin"),arglag=list(fun="poly", degree=4))
#mean temperature
lk = logknots(30,nk=3)
vk = equalknots(meanttemp,nk=2)
MeantTemp <- crossbasis(meanttemp,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#air pressure
vk = equalknots(ap,nk=2)
AP <- crossbasis(ap,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))
#relative humidity
vk = equalknots(rh,nk=2)
RH <- crossbasis(rh,lag=c(0,30),argvar=list(fun="bs",degree=2,knots=vk), arglag=list(knots=lk))

###############################model####################################
modela <-glm(YZ ~ AL + age +factor(gender) +factor(doy),family=binomial,data=wide)
modelb <-glm(YZ ~ AL + PM25+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelc <-glm(YZ ~ AL + PM10+age +factor(gender) +factor(doy),family=binomial,data=wide)
modeld <-glm(YZ ~ AL + MeantTemp+age +factor(gender) +factor(doy),family=binomial,data=wide)
modele <-glm(YZ ~ AL + AP+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelf <-glm(YZ ~ AL + RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelg <-glm(YZ ~ AL + PM25+PM10+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelh <-glm(YZ ~ AL + PM25+AP+age +factor(gender) +factor(doy),family=binomial,data=wide)
modeli <-glm(YZ ~ AL + PM25+MeantTemp+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelj <-glm(YZ ~ AL + PM25+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelk <-glm(YZ ~ AL + PM10+AP+age +factor(gender) +factor(doy),family=binomial,data=wide)
modell <-glm(YZ ~ AL + PM10+MeantTemp+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelm <-glm(YZ ~ AL + PM10+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
modeln <-glm(YZ ~ AL + AP+MeantTemp+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelo <-glm(YZ ~ AL + AP+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelp <-glm(YZ ~ AL + MeantTemp+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelq <-glm(YZ ~ AL + PM25+PM10+MeantTemp+AP+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelr <-glm(YZ ~ AL + PM25+PM10+MeantTemp+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
models <-glm(YZ ~ AL + PM25+PM10+AP+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
modelt <-glm(YZ ~ AL + PM10+AP+MeantTemp+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)
model <-glm(YZ ~ AL + PM25+PM10+AP+MeantTemp+RH+age +factor(gender) +factor(doy),family=binomial,data=wide)

###############################predict##################################
preda <- crosspred(AL,modela,cen=25,at=0:99,by=0.2,cumul = TRUE)
predb <- crosspred(AL,modelb,cen=25,at=0:99,by=0.2,cumul = TRUE)
predc <- crosspred(AL,modelc,cen=25,at=0:99,by=0.2,cumul = TRUE)
predd <- crosspred(AL,modeld,cen=25,at=0:99,by=0.2,cumul = TRUE)
prede <- crosspred(AL,modele,cen=25,at=0:99,by=0.2,cumul = TRUE)
predf <- crosspred(AL,modelf,cen=25,at=0:99,by=0.2,cumul = TRUE)
predg <- crosspred(AL,modelg,cen=25,at=0:99,by=0.2,cumul = TRUE)
predh <- crosspred(AL,modelh,cen=25,at=0:99,by=0.2,cumul = TRUE)
predi <- crosspred(AL,modeli,cen=25,at=0:99,by=0.2,cumul = TRUE)
predj <- crosspred(AL,modelj,cen=25,at=0:99,by=0.2,cumul = TRUE)
predk <- crosspred(AL,modelk,cen=25,at=0:99,by=0.2,cumul = TRUE)
predl <- crosspred(AL,modell,cen=25,at=0:99,by=0.2,cumul = TRUE)
predm <- crosspred(AL,modelm,cen=25,at=0:99,by=0.2,cumul = TRUE)
predn <- crosspred(AL,modeln,cen=25,at=0:99,by=0.2,cumul = TRUE)
predo <- crosspred(AL,modelo,cen=25,at=0:99,by=0.2,cumul = TRUE)
predp <- crosspred(AL,modelp,cen=25,at=0:99,by=0.2,cumul = TRUE)
predq <- crosspred(AL,modelq,cen=25,at=0:99,by=0.2,cumul = TRUE)
predr <- crosspred(AL,modelr,cen=25,at=0:99,by=0.2,cumul = TRUE)
preds <- crosspred(AL,models,cen=25,at=0:99,by=0.2,cumul = TRUE)
predt <- crosspred(AL,modelt,cen=25,at=0:99,by=0.2,cumul = TRUE)
pred <- crosspred(AL,model,cen=25,at=0:99,by=0.2,cumul = TRUE)

###############################plot##################################
tiff(file = "a.tiff", width = 3000, height = 2500, res = 300)
plot(preda,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "b.tiff", width = 3000, height = 2500, res = 300)
plot(predb,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "c.tiff", width = 3000, height = 2500, res = 300)
plot(predc,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "d.tiff", width = 3000, height = 2500, res = 300)
plot(predd,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "e.tiff", width = 3000, height = 2500, res = 300)
plot(prede,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "f.tiff", width = 3000, height = 2500, res = 300)
plot(predf,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()


tiff(file = "g.tiff", width = 3000, height = 2500, res = 300)
plot(predg,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "h.tiff", width = 3000, height = 2500, res = 300)
plot(predh,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "i.tiff", width = 3000, height = 2500, res = 300)
plot(predi,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "j.tiff", width = 3000, height = 2500, res = 300)
plot(predj,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "k.tiff", width = 3000, height = 2500, res = 300)
plot(predk,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "l.tiff", width = 3000, height = 2500, res = 300)
plot(predl,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "m.tiff", width = 3000, height = 2500, res = 300)
plot(predm,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "n.tiff", width = 3000, height = 2500, res = 300)
plot(predn,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "o.tiff", width = 3000, height = 2500, res = 300)
plot(predo,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "p.tiff", width = 3000, height = 2500, res = 300)
plot(predp,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "q.tiff", width = 3000, height = 2500, res = 300)
plot(predq,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "r.tiff", width = 3000, height = 2500, res = 300)
plot(predr,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "s.tiff", width = 3000, height = 2500, res = 300)
plot(preds,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "t.tiff", width = 3000, height = 2500, res = 300)
plot(predt,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

tiff(file = "all.tiff", width = 3000, height = 2500, res = 300)
plot(pred,"contour", xlab="ALAN", key.title=title("OR"),cex=5,cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title("Contour plot",xlab="ALAN",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()