library(openxlsx)
b<-read.xlsx("data1822.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)
library(qcc)
####possion
qcc.overdispersion.test(b$death, type = 'poisson')
####
quantile(b$death,probs = c(0,0.1,0.25,0.5,0.75,0.9,1))
quantile(b$PM25,probs = c(0,0.1,0.25,0.5,0.75,0.9,1))
quantile(b$PM10,probs = c(0,0.1,0.25,0.5,0.75,0.9,1))
quantile(b$NO2,probs = c(0,0.1,0.25,0.5,0.75,0.9,1))
quantile(b$CO,probs = c(0,0.1,0.25,0.5,0.75,0.9,1))
quantile(b$SO2,probs = c(0,0.1,0.25,0.5,0.75,0.9,1))
mean(b$death)
sd(b$death)
####qaicbic
qaicbic <- function(model1) {
  phi <- summary(model1)$dispersion
  logll <- sum(dpois(ceiling(model1$y), lambda=exp(predict(model1)), log=TRUE))
  cbind((-2*logll + 2*summary(model1)$df[3]*phi),
        (-2*logll + log(length(resid(model1)))*phi*summary(model1)$df[3]))}
tqba <- data.frame()
for (ii in 1:21) {
  cb.pm2.5<-crossbasis(b$PM25,lag=ii,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
  cb.TM<-crossbasis(b$temp,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
  model1<-glm(b$death~ cb.pm2.5+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
  summary(model1)
  q <- cbind(ii,qaicbic(model1))
  tqba <- rbind(tqba,q)
}
print(tqba)

####PM2.5
cb.TM<-crossbasis(b$temp,lag = 7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))
cb.pm2.5<-crossbasis(b$PM25,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
model1<-glm(b$death~ cb.pm2.5+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.pm2.5<-crosspred(cb.pm2.5,model1,cen=3.724206,at=3.724206:13.724206,bylag =1,cumul = T)
####
A1<-pred.pm2.5$matRRfit["13.724206",]
A2<-pred.pm2.5$matRRlow["13.724206",]
A3<-pred.pm2.5$matRRhigh["13.724206",]
A<-rbind(A1,A2,A3)
library(tidyr)
A4<-t(A)
A4<-as.data.frame(A4)
A4$A1<-round(A4$A1,3)
A4$A2<-round(A4$A2,3)
A4$A3<-round(A4$A3,3)
A4$a<-paste(A4$A1,A4$A2,sep = "(")
A4$b<-paste(A4$a,A4$A3,sep = "-")
A4$C<-paste(A4$b,"",")")
write.csv(A4,file = "A4.CSV")
####
A1<-pred.pm2.5$cumRRfit["13.724206",]
A2<-pred.pm2.5$cumRRlow["13.724206",]
A3<-pred.pm2.5$cumRRhigh["13.724206",]
A<-rbind(A1,A2,A3)
library(tidyr)
A4<-t(A)
A4<-as.data.frame(A4)
A4$A1<-round(A4$A1,3)
A4$A2<-round(A4$A2,3)
A4$A3<-round(A4$A3,3)
A4$a<-paste(A4$A1,A4$A2,sep = "(")
A4$b<-paste(A4$a,A4$A3,sep = "-")
A4$C<-paste(A4$b,"",")")
write.csv(A4,file = "B4.CSV")
####PM10
cb.pm10<-crossbasis(b$PM10,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
model2<-glm(b$death~ cb.pm10+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.pm10<-crosspred(cb.pm10,model2,cen=8.694444,at=8.694444:18.694444,bylag =1,cumul = T)
####
A1<-pred.pm10$matRRfit["18.694444",]
A2<-pred.pm10$matRRlow["18.694444",]
A3<-pred.pm10$matRRhigh["18.694444",]
A<-rbind(A1,A2,A3)
library(tidyr)
A4<-t(A)
A4<-as.data.frame(A4)
A4$A1<-round(A4$A1,3)
A4$A2<-round(A4$A2,3)
A4$A3<-round(A4$A3,3)
A4$a<-paste(A4$A1,A4$A2,sep = "(")
A4$b<-paste(A4$a,A4$A3,sep = "-")
A4$C<-paste(A4$b,"",")")
write.csv(A4,file = "A4.CSV")
####
A1<-pred.pm10$cumRRfit["18.694444",]
A2<-pred.pm10$cumRRlow["18.694444",]
A3<-pred.pm10$cumRRhigh["18.694444",]
A<-rbind(A1,A2,A3)
library(tidyr)
A4<-t(A)
A4<-as.data.frame(A4)
A4$A1<-round(A4$A1,3)
A4$A2<-round(A4$A2,3)
A4$A3<-round(A4$A3,3)
A4$a<-paste(A4$A1,A4$A2,sep = "(")
A4$b<-paste(A4$a,A4$A3,sep = "-")
A4$C<-paste(A4$b,"",")")
write.csv(A4,file = "B4.CSV")
####SO2
cb.SO2<-crossbasis(b$SO2,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
model3<-glm(b$death~ cb.SO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.SO2<-crosspred(cb.SO2,model3,cen = 2.4,bylag = 1,at=c(2.4:12.4),cumul = T)
A1<-pred.SO2$matRRfit["12.4",]
A2<-pred.SO2$matRRlow["12.4",]
A3<-pred.SO2$matRRhigh["12.4",]
A<-rbind(A1,A2,A3)
library(tidyr)
A4<-t(A)
A4<-as.data.frame(A4)
A4$A1<-round(A4$A1,3)
A4$A2<-round(A4$A2,3)
A4$A3<-round(A4$A3,3)
A4$a<-paste(A4$A1,A4$A2,sep = "(")
A4$b<-paste(A4$a,A4$A3,sep = "-")
A4$C<-paste(A4$b,"",")")
write.csv(A4,file = "A4.CSV")
B1<-pred.SO2$cumRRfit["12.4",]
B2<-pred.SO2$cumRRlow["12.4",]
B3<-pred.SO2$cumRRhigh["12.4",]
B<-rbind(B1,B2,B3)
B4<-t(B)
B4<-as.data.frame(B4)
B4$B1<-round(B4$B1,3)
B4$B2<-round(B4$B2,3)
B4$B3<-round(B4$B3,3)
B4$a<-paste(B4$B1,B4$B2,sep = "(")
B4$b<-paste(B4$a,B4$B3,sep = "-")
B4$C<-paste(B4$b,"",")")
write.csv(B4,file = "B4.CSV")
####CO
cb.CO<-crossbasis(b$CO,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
model4<-glm(b$death~ cb.CO+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.CO<-crosspred(cb.CO,model4,cen = 0.2208333,bylag = 1,at=c(0.2208333:1.2208333),cumul = T)
A1<-pred.CO$matRRfit["1.2208333",]
A2<-pred.CO$matRRlow["1.2208333",]
A3<-pred.CO$matRRhigh["1.2208333",]
A<-rbind(A1,A2,A3)
library(tidyr)
A4<-t(A)
A4<-as.data.frame(A4)
A4$A1<-round(A4$A1,3)
A4$A2<-round(A4$A2,3)
A4$A3<-round(A4$A3,3)
A4$a<-paste(A4$A1,A4$A2,sep = "(")
A4$b<-paste(A4$a,A4$A3,sep = "-")
A4$C<-paste(A4$b,"",")")
write.csv(A4,file = "A4.CSV")
B1<-pred.CO$cumRRfit["1.2208333",]
B2<-pred.CO$cumRRlow["1.2208333",]
B3<-pred.CO$cumRRhigh["1.2208333",]
B<-rbind(B1,B2,B3)
B4<-t(B)
B4<-as.data.frame(B4)
B4$B1<-round(B4$B1,3)
B4$B2<-round(B4$B2,3)
B4$B3<-round(B4$B3,3)
B4$a<-paste(B4$B1,B4$B2,sep = "(")
B4$b<-paste(B4$a,B4$B3,sep = "-")
B4$C<-paste(B4$b,"",")")
write.csv(B4,file = "B4.CSV")
####NO2
cb.NO2<-crossbasis(b$NO2,lag=7,argvar = list(fun="lin"),arglag = list(fun="ns",df=3))
model5<-glm(b$death~ cb.NO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.NO2<-crosspred(cb.NO2,model5,cen=4.65625  ,bylag = 1,at=c(4.65625  :14.65625  ),cumul = T)
A1<-pred.NO2$matRRfit["14.65625",]
A2<-pred.NO2$matRRlow["14.65625",]
A3<-pred.NO2$matRRhigh["14.65625",]
A<-rbind(A1,A2,A3)
library(tidyr)
A4<-t(A)
A4<-as.data.frame(A4)
A4$A1<-round(A4$A1,3)
A4$A2<-round(A4$A2,3)
A4$A3<-round(A4$A3,3)
A4$a<-paste(A4$A1,A4$A2,sep = "(")
A4$b<-paste(A4$a,A4$A3,sep = "-")
A4$C<-paste(A4$b,"",")")
write.csv(A4,file = "A4.CSV")
B1<-pred.NO2$cumRRfit["14.65625",]
B2<-pred.NO2$cumRRlow["14.65625",]
B3<-pred.NO2$cumRRhigh["14.65625",]
B<-rbind(B1,B2,B3)
B4<-t(B)
B4<-as.data.frame(B4)
B4$B1<-round(B4$B1,3)
B4$B2<-round(B4$B2,3)
B4$B3<-round(B4$B3,3)
B4$a<-paste(B4$B1,B4$B2,sep = "(")
B4$b<-paste(B4$a,B4$B3,sep = "-")
B4$C<-paste(B4$b,"",")")
write.csv(B4,file = "B4.CSV")
####
library(ggplot2)
library(eoffice)
rm(list=ls())
b<-read.xlsx("XXXX.xlsx")
b$x<-rep(c("0","1","2","3","4","5","6","7","01","02","03","04","05","06","07"),6)
b$x<-factor(b$x,levels = c("0","1","2","3","4","5","6","7","01","02","03","04","05","06","07"),ordered=T)
b$group3<-paste(b$group1,b$group2,sep = "")
names(b)
b<- within(b, group3 <- factor(group3, levels = c("PM2.5single","PM10single","PM2.5cum","PM10cum","SO2single","NO2single","SO2cum","NO2cum",
                                                  "COsingle","xxsingle","COcum","xxcum")))
with(b, levels(group3))
p<-ggplot(b,aes(x=x,y=(RR-1)*100,color=group2))+
  geom_point(position=position_dodge(width=0.5))+
  geom_errorbar(aes(min=(lower-1)*100,max=(higher-1)*100),width=.2,position=position_dodge(0.5),lwd=0.8)+
  labs(x="Lag Days", y= 'Percentage change')+ 
  geom_hline(yintercept = 0,lty='dashed')+
  theme(text= element_text(size=18),axis.text = element_text(colour = "black"))+
  scale_color_manual(values = c('red', "blue"))+
  theme(panel.background=element_rect(fill = "white"))+
  theme(axis.line.x = element_line(color="black", linewidth = 0.8),axis.line.y = element_line(color="black", size = 0.8))+
  facet_wrap(group3~.,scales = "free")+
  theme(legend.position = "bottom")
p
topptx(filename = file.path("XXXX.pptx"),width = 10,
       height = 7,units="in")
dev.off()
####
cb.TM<-crossbasis(b$temp,lag = 7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))
cb.pm2.5<-crossbasis(b$PM25,lag=7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))
model1<-glm(b$death~ cb.pm2.5+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.pm2.5<-crosspred(cb.pm2.5,model1,cen=3.724206,at=3.724206:241.111111 ,bylag =1,cumul = T)
cb.pm10<-crossbasis(b$PM10,lag=7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))
model2<-glm(b$death~ cb.pm10+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.pm10<-crosspred(cb.pm10,model2,cen=8.694444,at=8.694444:1017.135417 ,bylag =1,cumul = T)
cb.SO2<-crossbasis(b$SO2,lag=7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))
model3<-glm(b$death~ cb.SO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.SO2<-crosspred(cb.SO2,model3,cen = 2.4,bylag = 1,at=c(2.4:49.763889 ),cumul = T)
cb.CO<-crossbasis(b$CO,lag=7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))
model4<-glm(b$death~ cb.CO+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.CO<-crosspred(cb.CO,model4,cen = 0.2208333,bylag = 1,at=c(0.2208333:2.8111111 ),cumul = T)
cb.NO2<-crossbasis(b$NO2,lag=7,argvar = list(fun="ns",df=3),arglag = list(fun="ns",df=3))
model5<-glm(b$death~ cb.NO2+cb.TM+ns(b$time,df=7*5)+b$holiday+b$dow,family=quasipoisson,data=b)
pred.NO2<-crosspred(cb.NO2,model5,cen=4.65625  ,bylag = 1,at=c(4.65625  :108.30556   ),cumul = T)
###3Dͼ
plot(pred.pm2.5,ticktype="detailed",border="#3366FF",xlab="PM2.5(ug/m3)",ylab="Lag",zlab="RR",
     col="#99FFCC",shade=0.1,cex.lab=1.3,cex.axis=1.3,lwd=1,theta=41,phi=25,ltheta=-35)
topptx(filename = file.path("D:XXXX","3dpm2.5.pptx"),width = 10,
       height = 8,units="in")
dev.off()
plot(pred.pm10,ticktype="detailed",border="#3366FF",xlab="PM10(ug/m3)",ylab="Lag",zlab="RR",
     col="#99FFCC",shade=0.1,cex.lab=1.3,cex.axis=1.3,lwd=1,theta=41,phi=25,ltheta=-35)
plot(pred.SO2,ticktype="detailed",border="#3366FF",xlab="SO2(ug/m3)",ylab="Lag",zlab="RR",
     col="#99FFCC",shade=0.1,cex.lab=1.3,cex.axis=1.3,lwd=1,theta=41,phi=25,ltheta=-35)
plot(pred.CO,ticktype="detailed",border="#3366FF",xlab="CO(mg/m3)",ylab="Lag",zlab="RR",
     col="#99FFCC",shade=0.1,cex.lab=1.3,cex.axis=1.3,lwd=1,theta=41,phi=25,ltheta=-35)
plot(pred.NO2,ticktype="detailed",border="#3366FF",xlab="NO2(ug/m3)",ylab="Lag",zlab="RR",
     col="#99FFCC",shade=0.1,cex.lab=1.3,cex.axis=1.3,lwd=1,theta=41,phi=25,ltheta=-35)
plot(pred.O3,ticktype="detailed",border="#3366FF",xlab="O3(ug/m3)",ylab="Lag",zlab="RR",
     col="#99FFCC",shade=0.1,cex.lab=1.3,cex.axis=1.3,lwd=1,theta=41,phi=25,ltheta=-35)
####?ȸ???ͼ
plot(pred.pm2.5,"contour",xlab="PM2.5",key.title=title("RR"),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(xlab="PM2.5(ug/m3)",ylab="Lag",cex.main=2,cex.lab=1.5))
topptx(filename = file.path("XXX,pm2.5.pptx"),width = 10,
       height = 8,units="in")
dev.off()
plot(pred.pm10,"contour",xlab="PM10",key.title=title("RR"),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(xlab="PM10(ug/m3)",ylab="Lag",cex.main=2,cex.lab=1.5))
#overallͼ
library(eoffice)
par(mfrow=c(2,3))
plot(pred.pm2.5,"overall",col="#DC0000E5",
     xlab="PM2.5",ylab=expression("Effect "),
     main="Overall Effect", lwd=2)
plot(pred.pm10,"overall",col="#DC0000E5",
     xlab="PM10",ylab=expression("Effect "),
     main="Overall Effect", lwd=2)
plot(pred.SO2,"overall",col="#DC0000E5",
     xlab="SO2",ylab=expression("Effect "),
     main="Overall Effect", lwd=2)
plot(pred.NO2,"overall",col="#DC0000E5",
     xlab="NO2",ylab=expression("Effect "),
     main="Overall Effect", lwd=2)
plot(pred.CO,"overall",col="#DC0000E5",
     xlab="CO",ylab=expression("Effect "),
     main="Overall Effect", lwd=2)
topptx(filename = file.path("D:\\XXX","ns.pptx"),width = 12,
       height = 8,units="in")
dev.off()
####
library(eoffice)
par(mfrow=c(2,3))
crall1<- crossreduce(cb.pm2.5,model1,cen=34,type="overall",lag=c(0,5))
plot(crall1,xlab="PM2.5",ylab="Relative risk",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="",cex=0.89)
crall2<- crossreduce(cb.pm10,model2,cen=54,type="overall",lag=c(0,5))
plot(crall2,xlab="PM10",ylab="Relative risk",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="",cex=0.89)
crall3<- crossreduce(cb.SO2,model3,cen=10,type="overall",lag=c(0,5))
plot(crall3,xlab="SO2",ylab="Relative risk",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="",cex=0.89)
crall4<- crossreduce(cb.NO2,model5,cen=36,type="overall",lag=c(0,5))
plot(crall4,xlab="NO2",ylab="Relative risk",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="",cex=0.89)
crall5<- crossreduce(cb.CO,model4,cen=0.7,type="overall",lag=c(0,5))
plot(crall5,xlab="CO",ylab="Relative risk",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="",cex=0.89)
crall6<- crossreduce(cb.O3,model6,cen=29,type="overall",lag=c(0,5))
plot(crall6,xlab="O3",ylab="Relative risk",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="",cex=0.89)
topptx(filename = file.path("D:\\XXXX","xx.pptx"),width = 8,
       height = 5,units="in")
dev.off()