### seasonal analysis
### space of predictor: linear effect above 40.3 microgr/m3 for O3
### space of predictor: linear effects below 15C and above 25C for temperature
### lag function: integer lag parameterization (unconstrained) for O3 up to lag5
### lag function: strata intervals at lag 0-1, 2-5 and 6-10 for temperature
# SELECT SUMMER MONTHS OF THE SERIES
chicagoNMMAPSseas <- subset(chicagoNMMAPS, month>5 & month<10)
# CREATE THE CROSS-BASIS FOR EACH PREDICTOR, SPECIFYING THE GROUP VARIABLE
cb2.o3 <- crossbasis(chicagoNMMAPSseas$o3, lag=5, argvar=list(type="hthr",
knots=40.3), arglag=list(type="integer"), group=chicagoNMMAPSseas$year)
cb2.temp <- crossbasis(chicagoNMMAPSseas$temp, lag=10,
argvar=list(type="dthr",knots=c(15,25)), arglag=list(type="strata",
knots=c(2,6)), group=chicagoNMMAPSseas$year)
summary(cb2.o3)
summary(cb2.temp)
# RUN THE MODEL AND GET THE PREDICTION FOR 03
library(splines)
model2 <- glm(death ~ cb2.o3 + cb2.temp + ns(doy, 4) + dow,
family=quasipoisson(), chicagoNMMAPSseas)
pred2.o3 <- crosspred(cb2.o3, model2, at=c(0:65,40.3,50.3))
# PLOT THE LINEAR ASSOCIATION OF 03 ALONG LAGS (WITH 80%CI)
plot(pred2.o3, "slices", var=50.3, ci="bars", type="p", pch=19, ci.level=0.80,
main="Effects of 10-unit increase above the threshold (80CI)")
# PLOT THE OVERALL ASSOCIATION
plot(pred2.o3,"overall",xlab="Ozone", ci="lines", ylim=c(0.9,1.3), lwd=2,
ci.arg=list(col=1,lty=3), main="Overall effects over 5 days of lag")
# GET THE FIGURES FOR THE OVERALL ASSOCIATION ABOVE, WITH CI
pred2.o3$allRRfit["50.3"]
cbind(pred2.o3$allRRlow, pred2.o3$allRRhigh)["50.3",]
# 3-D PLOT WITH DEFAULT AND USER-DEFINED PERSPECTIVE
plot(pred2.o3, xlab="Ozone", main="3D: default perspective")
plot(pred2.o3, xlab="Ozone", main="3D: different perspective", theta=250, phi=40)
### See the vignette 'dlnmOverview' for a detailed explanation of this example
Run the code above in your browser using DataLab