library(openxlsx)
b<-read.xlsx("data18-22IQR.xlsx")
b$date<-as.Date(b$date,origin="1899-12-30")
b$dow<-ifelse(b$dow=="Saturday"|b$dow=="Sunday",0,1)
b$dow<-factor(b$dow,levels = c(0,1),labels = c("weekend","weekday"))
b$holiday<-factor(b$holiday,levels = c(0,1),labels = c("nonholiday","holiday"))
library(mgcv)
library(splines)
library(ggplot2)
library(dlnm)
library(reshape2)

####
cb.pm2.5<-crossbasis(b$PM25,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
cb.pm10<-crossbasis(b$PM10,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
cb.SO2<-crossbasis(b$SO2,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
cb.CO<-crossbasis(b$CO,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
cb.NO2<-crossbasis(b$NO2,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
cb.TM<-crossbasis(b$temp,lag = 7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))

#####PM2.5 and SO2
quantile(b$PM25,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.pm2.5+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.pm2.5+cb.SO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.pm2.5<-crosspred(cb.pm2.5,fit0,cen=0.5467174,at=0.5467174:1.5467174,bylag =1,cumul = T)
pred.pm2.5<-crosspred(cb.pm2.5,fit1,cen=0.5467174,at=0.5467174:1.5467174,bylag =1,cumul = T)

#####
E1<-pred.pm2.5$cumRRfit["1.5467174",]
E2<-pred.pm2.5$cumRRlow["1.5467174",]
E3<-pred.pm2.5$cumRRhigh["1.5467174",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")

A1<-pred0.pm2.5$cumfit["1.5467174",3]
A2<-pred0.pm2.5$cumse["1.5467174",3]
A3<-pred.pm2.5$cumfit["1.5467174",3]
A4<-pred.pm2.5$cumse["1.5467174",3]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A)

#####PM2.5 and NO2
quantile(b$PM25,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.pm2.5+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.pm2.5+cb.NO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.pm2.5<-crosspred(cb.pm2.5,fit0,cen=0.5467174,at=0.5467174:1.5467174,bylag =1,cumul = T)
pred.pm2.5<-crosspred(cb.pm2.5,fit1,cen=0.5467174,at=0.5467174:1.5467174,bylag =1,cumul = T)
E1<-pred.pm2.5$cumRRfit["1.5467174",]
E2<-pred.pm2.5$cumRRlow["1.5467174",]
E3<-pred.pm2.5$cumRRhigh["1.5467174",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")

A1<-pred0.pm2.5$cumfit["1.5467174",3]
A2<-pred0.pm2.5$cumse["1.5467174",3]
A3<-pred.pm2.5$cumfit["1.5467174",3]
A4<-pred.pm2.5$cumse["1.5467174",3]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A) 


#####PM10 and SO2
quantile(b$PM10,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.pm10+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.pm10+cb.SO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.pm10<-crosspred(cb.pm10,fit0,cen=0.7487055,at=0.7487055:1.7487055,bylag =1,cumul = T)
pred.pm10<-crosspred(cb.pm10,fit1,cen=0.7487055,at=0.7487055:1.7487055,bylag =1,cumul = T)

E1<-pred.pm10$cumRRfit["1.7487055",]
E2<-pred.pm10$cumRRlow["1.7487055",]
E3<-pred.pm10$cumRRhigh["1.7487055",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")

A1<-pred0.pm10$cumfit["1.7487055",3]
A2<-pred0.pm10$cumse["1.7487055",3]
A3<-pred.pm10$cumfit["1.7487055",3]
A4<-pred.pm10$cumse["1.7487055",3]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A) 

####PM10 and NO2
quantile(b$PM10,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.pm10+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.pm10+cb.NO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.pm10<-crosspred(cb.pm10,fit0,cen=0.7487055,at=0.7487055:1.7487055,bylag =1,cumul = T)
pred.pm10<-crosspred(cb.pm10,fit1,cen=0.7487055,at=0.7487055:1.7487055,bylag =1,cumul = T)

E1<-pred.pm10$cumRRfit["1.7487055",]
E2<-pred.pm10$cumRRlow["1.7487055",]
E3<-pred.pm10$cumRRhigh["1.7487055",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")

A1<-pred0.pm10$cumfit["1.7487055",3]
A2<-pred0.pm10$cumse["1.7487055",3]
A3<-pred.pm10$cumfit["1.7487055",3]
A4<-pred.pm10$cumse["1.7487055",3]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A) 

#####PM10 and CO
quantile(b$PM10,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.pm10+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.pm10+cb.CO+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.pm10<-crosspred(cb.pm10,fit0,cen=0.7487055,at=0.7487055:1.7487055,bylag =1,cumul = T)
pred.pm10<-crosspred(cb.pm10,fit1,cen=0.7487055,at=0.7487055:1.7487055,bylag =1,cumul = T)

E1<-pred.pm10$cumRRfit["1.7487055",]
E2<-pred.pm10$cumRRlow["1.7487055",]
E3<-pred.pm10$cumRRhigh["1.7487055",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")

A1<-pred0.pm10$cumfit["1.7487055",3]
A2<-pred0.pm10$cumse["1.7487055",3]
A3<-pred.pm10$cumfit["1.7487055",3]
A4<-pred.pm10$cumse["1.7487055",3]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A) 


#####NO2 and CO
quantile(b$NO2,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.NO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.NO2+cb.CO+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.NO2<-crosspred(cb.NO2,fit0,cen=0.5748137,bylag = 1,at=c(0.5748137:1.5748137),cumul = T)
pred.NO2<-crosspred(cb.NO2,fit1,cen=0.5748137,bylag = 1,at=c(0.5748137:1.5748137),cumul = T)

E1<-pred.NO2$cumRRfit["1.5748137",]
E2<-pred.NO2$cumRRlow["1.5748137",]
E3<-pred.NO2$cumRRhigh["1.5748137",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")

A1<-pred0.NO2$cumfit["1.5748137",6]
A2<-pred0.NO2$cumse["1.5748137",6]
A3<-pred.NO2$cumfit["1.5748137",6]
A4<-pred.NO2$cumse["1.5748137",6]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A)

####NO2 and SO2
quantile(b$NO2,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.NO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.NO2+cb.SO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.NO2<-crosspred(cb.NO2,fit0,cen=0.5748137,bylag = 1,at=c(0.5748137:1.5748137),cumul = T)
pred.NO2<-crosspred(cb.NO2,fit1,cen=0.5748137,bylag = 1,at=c(0.5748137:1.5748137),cumul = T)
E1<-pred.NO2$cumRRfit["1.5748137",]
E2<-pred.NO2$cumRRlow["1.5748137",]
E3<-pred.NO2$cumRRhigh["1.5748137",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")


A1<-pred0.NO2$cumfit["1.5748137",6]
A2<-pred0.NO2$cumse["1.5748137",6]
A3<-pred.NO2$cumfit["1.5748137",6]
A4<-pred.NO2$cumse["1.5748137",6]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A) 


#####CO and SO2
quantile(b$CO,probs = c(0,0.25,0.5,0.75,1))
fit0<-glm(b$death~ cb.CO+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
fit1<-glm(b$death~ cb.CO+cb.SO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred0.CO<-crosspred(cb.CO,fit0,cen = 1.7721274,bylag = 1,at=c(1.7721274:2.7721274),cumul = T)
pred.CO<-crosspred(cb.CO,fit1,cen = 1.7721274,bylag = 1,at=c(1.7721274:2.7721274),cumul = T)
E1<-pred.CO$cumRRfit["2.7721274",]
E2<-pred.CO$cumRRlow["2.7721274",]
E3<-pred.CO$cumRRhigh["2.7721274",]
E<-rbind(E1,E2,E3)
library(tidyr)
E4<-t(E)
E4<-as.data.frame(E4)
E4$E1<-round(E4$E1,3)
E4$E2<-round(E4$E2,3)
E4$E3<-round(E4$E3,3)
E4$a<-paste(E4$E1,E4$E2,sep = "(")
E4$b<-paste(E4$a,E4$E3,sep = "-")
E4$C<-paste(E4$b,"",")")
write.csv(E4,file = "E4.CSV")


A1<-pred0.CO$cumfit["2.7721274",3] 
A2<-pred0.CO$cumse["2.7721274",3]
A3<-pred.CO$cumfit["2.7721274",3]
A4<-pred.CO$cumse["2.7721274",3]
A<-rbind(A1,A2,A3,A4)
library(tidyr)
A5<-t(A) 

####
z<-read.xlsx("XXX.xlsx")
z$z<-(z$β2-z$β1)/(z$se2^2+z$se1^2)^0.5
z$p<-pnorm(-abs(z$z))*2