mydata<-read.table("database.csv",header=T, sep=",") head(mydata,3) mydata$airpollution<-factor(mydata$airpollution,order=T) library(dlnm) cb1.pm1<-crossbasis(mydata$airpullition,lag=7,argvar=list(fun="lin"),arglag=list(fun="poly",degree=3)) cb1.pm2<-crossbasis(mydata$PM2.5,lag=7,argvar=list(fun="lin"),arglag=list(fun="poly",degree=3)) cb1.pm3<-crossbasis(mydata$O3,lag=7,argvar=list(fun="lin"),arglag=list(fun="poly",degree=3)) library(splines) library(mgcv) library(nlme) model1<-gam(morbidity~cb1.pm1+cb1.pm2+cb1.pm3+ns(precipitation,3)+ns(temperature,3)+ns(airpressure,3)+ns(time,7)+as.factor(week)+as.factor(holiday),family=quasipoisson(),data=mydata) pred1.pm1<-crosspred(cb1.pm1, model1, at=c(1,2,3,4), by=1,cumul=TRUE,cen=0) plot(pred1.pm1, "slices", var=1, col=3, ylab="RR", ci.arg=list(density=15,lwd=2), main="Association with mild air pollution") plot(pred1.pm1, "slices", var=1, cumul=TRUE, ylab="Cumulative RR", main="Cumulative association with mild air pollution") pred1.pm2<-crosspred(cb1.pm2, model1, at=c(0:304), by=1,cumul=TRUE,cen=75) plot(pred1.pm2, "slices", var=85, col=3, ylab="RR", ci.arg=list(density=15,lwd=2), main="Association with a 10-unit increase in PM2.5") plot(pred1.pm2, "slices", var=85, cumul=TRUE, ylab="Cumulative RR", main="Cumulative association with a 10-unit increase in PM2.5") cbind(pred1.pm2$matRRfit)["85",] cbind(pred1.pm2$matRRlow, pred1.pm2$matRRhigh)["85",] cbind(pred1.pm2$cumRRfit)["85",] cbind(pred1.pm2$cumRRlow, pred1.pm2$cumRRhigh)["85",] pred1.pm3<-crosspred(cb1.pm3, model1, at=c(0:254), by=1,cumul=TRUE,cen=200) plot(pred1.pm7, "slices", var=210, col=3, ylab="RR", ci.arg=list(density=15,lwd=2), main="Association with a 10-unit increase in O3") plot(pred1.pm7, "slices", var=210, cumul=TRUE, ylab="Cumulative RR", main="Cumulative association with a 10-unit increase in O3") cbind(pred1.pm3$matRRfit)["210",] cbind(pred1.pm3$matRRlow, pred1.pm3$matRRhigh)["210",] cbind(pred1.pm3$cumRRfit)["210",] cbind(pred1.pm3$cumRRlow, pred1.pm3$cumRRhigh)["210",]