# Example 2. See crossbasis and plot.crosspred for other examples
### 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
basis.o3 <- crossbasis(chicagoNMMAPSseas$o3, group=chicagoNMMAPSseas$year,
vartype="hthr", varknots=40.3, lagtype="integer", maxlag=5)
basis.temp <- crossbasis(chicagoNMMAPSseas$temp, group=chicagoNMMAPSseas$year,
vartype="dthr", varknots=c(15,25), lagtype="strata",
lagknots=c(2,6), maxlag=10)
summary(basis.o3)
summary(basis.temp)
# RUN THE MODEL AND GET THE PREDICTION FOR 03
library(splines)
model <- glm(death ~ basis.o3 + basis.temp + ns(doy, 4) + dow,
family=quasipoisson(), chicagoNMMAPSseas)
pred.o3 <- crosspred(basis.o3, model, at=c(0:65,40.3,50.3))
# PLOT THE LINEAR EFFECTS OF 03 ABOVE THE THRESHOLD (WITH 80%CI)
plot(pred.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(pred.o3,"overall",xlab="Ozone", ci="lines", ylim=c(0.9,1.3), lwd=2,
ci.arg=list(col=1,lty=3), main="Overall effect over 5 days of lag")
# GET THE FIGURES FOR THE OVERALL EFFECT ABOVE, WITH CI
pred.o3$allRRfit["50.3"]
cbind(pred.o3$allRRlow, pred.o3$allRRhigh)["50.3",]
# 3-D PLOT WITH DEFAULT AND USER-DEFINED PERSPECTIVE
plot(pred.o3, xlab="Ozone", main="3D graph of ozone effect")
plot(pred.o3, xlab="Ozone", main="3D graph of ozone effect", theta=250, phi=40)
### See the vignette 'dlnmOverview' for a detailed explanation of this exampleRun the code above in your browser using DataLab