################################################### ### chunk number 1: ################################################### options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) ################################################### ### chunk number 2: load ################################################### library("dlnm") ################################################### ### chunk number 3: vignette eval=FALSE ################################################### ## vignette("dlnmOverview") ################################################### ### chunk number 4: internal.mkbasis ################################################### mkbasis(1:5, type="bs", df=4, degree=2, cenvalue=3) ################################################### ### chunk number 5: internal.mklagbasis ################################################### mklagbasis(maxlag=5, type="strata", knots=c(2,4)) ################################################### ### chunk number 6: crossbasis ################################################### basis.o3 <- crossbasis(chicagoNMMAPS$o3, vartype="hthr", varknots=40.3, lagtype="strata", lagknots=c(2,6), maxlag=10) basis.temp <- crossbasis(chicagoNMMAPS$temp, vartype="bs", vardegree=3, vardf=6, cenvalue=25, lagdf=5, maxlag=30) ################################################### ### chunk number 7: summary ################################################### summary(basis.temp) ################################################### ### chunk number 8: model ################################################### library("splines") model <- glm(death ~ basis.temp + basis.o3 + ns(time,7*14) + dow, family=quasipoisson(), chicagoNMMAPS) ################################################### ### chunk number 9: crosspred ################################################### pred.o3 <- crosspred(basis.o3, model, at=c(0:65,40.3,50.3)) pred.temp <- crosspred(basis.temp, model, by=2) ################################################### ### chunk number 10: o3effect ################################################### pred.o3$allRRfit["50.3"] cbind(pred.o3$allRRlow,pred.o3$allRRhigh)["50.3",] ################################################### ### chunk number 11: o3slice ################################################### plot(pred.o3, "slices", type="p", pch=19, cex=1.5, var=50.3, ci="bars", ylab="RR",main="Lag-specific effects") ################################################### ### chunk number 12: o3overall ################################################### plot(pred.o3, "overall", ci="lines", ylim=c(0.95,1.25), lwd=2, col=4, xlab="Ozone", ylab="RR", main="Overall effect") ################################################### ### chunk number 13: temp3d ################################################### plot(pred.temp, xlab="Temperature", theta=240, phi=40, ltheta=-185, zlab="RR", main="3D graph") ################################################### ### chunk number 14: tempcontour ################################################### plot(pred.temp, "contour", plot.title=title(xlab="Temperature", ylab="Lag", main="Contour graph"), key.title=title("RR")) ################################################### ### chunk number 15: tempslices1 ################################################### plot(pred.temp, "slices", var=-20, ci="n",ylim=c(0.95,1.22), lwd=1.5) for(i in 1:2) lines(pred.temp, "slices", var=c(0,32)[i], col=i+2, lwd=1.5) legend("topright",paste("Temperature =",c(-20,0,32)), col=2:4, lwd=1.5) ################################################### ### chunk number 16: tempslices2 ################################################### plot(pred.temp, "slices", var=c(-20,0,32), lag=c(0,5,20), ci.level=0.99, xlab="Temperature",ci.arg=list(density=20,col=grey(0.7))) ################################################### ### chunk number 17: selectionprep ################################################### basis.temp2 <- crossbasis(chicagoNMMAPS$temp, vartype="poly", vardegree=6, cenvalue=25, lagdf=5, maxlag=30) model2 <- update(model, .~. - basis.temp + basis.temp2) pred.temp2 <- crosspred(basis.temp2, model2, by=2) basis.temp3 <- crossbasis(chicagoNMMAPS$temp, vartype="dthr", varknots=25, lagdf=5, maxlag=30) model3 <- update(model, .~. - basis.temp + basis.temp3) pred.temp3 <- crosspred(basis.temp3, model3, by=2) ################################################### ### chunk number 18: selectionplot1 ################################################### plot(pred.temp, "overall", ylim=c(0.5,2.5), ci="n", lwd=1.5, main="Overall effect") lines(pred.temp2, "overall", col=3, lty=2, lwd=2) lines(pred.temp3, "overall", col=4, lty=4, lwd=2) legend("top", c("natural spline","polynomial","double threshold"), col=2:4, lty=c(1:2,4), lwd=1.5, inset=0.1, cex=0.8) ################################################### ### chunk number 19: selectionplot2 ################################################### plot(pred.temp, "slices", var=32, ylim=c(0.95,1.22), ci="n", lwd=1.5, main="Lag-specific effect") lines(pred.temp2, "slices", var=32, col=3, lty=2, lwd=2) lines(pred.temp3, "slices", var=32, col=4, lty=4, lwd=2) legend("top", c("natural spline","polynomial","double threshold"), col=2:4, lty=c(1:2,4), inset=0.1, cex=0.8)