Nothing
## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(season)
library(MASS) # for mvrnorm
library(survival) # for coxph
library(ggplot2)
library(dlnm) # for splines
## -----------------------------------------------------------------------------
data(CVDdaily)
ggplot(CVDdaily, aes(x=factor(month), y=cvd))+
geom_boxplot()+
ylab('Daily number of CVD deaths')+
xlab('Month')+
theme_bw()
## ---- fig.width=7-------------------------------------------------------------
years = 1987:2001
Januarys = as.numeric(as.Date(paste(years, '-01-01', sep=''), origin='1970-01-01'))
ggplot(CVDdaily, aes(x=as.numeric(date), y=cvd))+
geom_line()+
scale_x_continuous(breaks=Januarys, labels=years)+
ylab('Daily number of CVD deaths')+
xlab('Time')+
theme_bw()+
theme(panel.grid.minor = element_blank())
## ---- fig.width=6-------------------------------------------------------------
ggplot(CVDdaily, aes(x=tmpd, y=cvd))+
geom_point()+
ylab('Daily number of CVD deaths')+
xlab('Temperature (degrees F)')+
theme_bw()
## ---- warning=FALSE, message=FALSE--------------------------------------------
# make a spline basis that has a lag and is non-linear
tmpd.basis = crossbasis(CVDdaily$tmpd, lag=14, # 14 day lag
arglag=list(fun='ns', df=3), # 3 degrees of freedom for lag; ns = natural spline
argvar=list(fun='ns', knots=c(60, 75))) # knots at 65 and 75 degrees
# add the spline basis variables to the data
CVDdaily = cbind(CVDdaily, tmpd.basis[1:nrow(CVDdaily), ])
# create the regression formula
spline.names = colnames(tmpd.basis)
formula = paste('cvd ~', paste(spline.names, collapse = ' + '), '+ Mon + Tue + Wed + Thu + Fri + Sat')
model = casecross(as.formula(formula), data=CVDdaily)
summary(model)
## -----------------------------------------------------------------------------
confint(model$c.model)
## ---- fig.height=5, fig.width=5-----------------------------------------------
# extract the coefficients and variance--covariance matrix for the spline terms
coef = coefficients(model$c.model)
index = names(coef) %in% spline.names
coef = coef[index]
vcov = vcov(model$c.model)[index, index]
for.plot = crosspred(basis=tmpd.basis, coef=coef, vcov=vcov, at=seq(45, 85, 1), model.link = 'log', cen=70)
par(mai=c(0.2, 0, 0, 0)) # reduce plot margins
plot(for.plot, xlab='Temperature (degrees F)', zlab='Odds ratio', ylab='Lag (days)')
## ---- fig.width=5-------------------------------------------------------------
to.plot = data.frame(temperature = for.plot$predvar,
mean = for.plot$allRRfit,
lower = for.plot$allRRlow,
upper = for.plot$allRRhigh)
ggplot(data=to.plot, aes(x=temperature, y=mean, ymin=lower, ymax=upper))+
geom_hline(lty=2, yintercept = 1)+ # horizontal reference line at no change in odds
geom_ribbon(alpha=0.3)+
geom_line()+
xlab('Temperature (degrees F)')+
ylab('Odds ratio')+
theme_bw()
## ---- fig.width=7, warning=FALSE----------------------------------------------
set.seed(1234) # set the random seed to give repeatable results
data(CVD)
f = c(12) # a single twelve month cycle
tau = c(10, 50) # achieved via trial-and-error; small tau -> less variability
ns.season = nscosinor(data=CVD, response='adj', cycles=f, niters=2000,
burnin=500, tau=tau, div=1000)
summary(ns.season)
plot(ns.season)
## -----------------------------------------------------------------------------
cat('Mean phase = ', round(invyrfraction(0.6952055/(2*pi), type='monthly', text=F), 2), ' months.\n', sep='')
cat('Lower 95% interval = ', round(invyrfraction(0.5732958/(2*pi), type='monthly', text=F), 2), ' months.\n', sep='')
cat('Upper 95% interval = ', round(invyrfraction(0.8216251/(2*pi), type='monthly', text=F), 2), ' months.\n', sep='')
## -----------------------------------------------------------------------------
set.seed(1234) # set the random seed to give repeatable results
ntest.residuals = nonlintest(ns.season$residuals, n.lag=12, n.boot=500)
ntest.residuals
## ---- warnings=FALSE, fig.width=5, fig.height=5-------------------------------
plot = plot(ntest.residuals, plot=FALSE)
plot + scale_x_continuous(breaks = 0:12) +
scale_y_continuous(breaks = 0:12) +
theme_bw()+
xlab('Lag 1 (months)')+
ylab('Lag 2 (months)')+
theme(panel.grid.minor = element_blank())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.