# sample analysis code using Richard Smith's "fitTS" analysis function

# first load data and separate into 8 air basins

# if necessary, put the full file path in the following command
X=read.csv('for-analysis.csv',header=T)
X1=X[X$basin=='mountain-counties',]
X2=X[X$basin=='sacramento-valley',]
X3=X[X$basin=='salton-sea',]
X4=X[X$basin=='san-diego',]
X5=X[X$basin=='san-francisco-bay',]
X6=X[X$basin=='san-joaquin-valley',]
X7=X[X$basin=='south-central-coast',]
X8=X[X$basin=='south-coast',]

# Example of fitTS applied to deaths from South Coast (all nonaccidental deaths aged 65+), using tmax, tmin and MAXRH as
# meteorological covariates, ozone (in ppm) as air pollution covariate. The program produces a plot of ozone risk at
# levels from 0.02 to 0.15 (ppm) at an interval of 0.001, using 0.075 as the reference level for which the relative
# risk is 1. The program is based on lag 0, uses 7 df per year for the long-term trend function, 6 df for both
# the current day and lagged days of meteorology, and uses ndf3=1 to denote a linear C-R curve (ndf3>1 would indicate
# a nonlinear curve with ndf3 degrees of freedom

f1=fitTS(X8$alldeathnoac65,cbind(X8$tmax,X8$tmin,X8$MAXRH),0.001*X8$o3DMOL8N,0.02,0.15,0.001,0.075,0,7,6,6,1)

# show results of LRT for meteorological variables
f1$lrt
# result: 0.000000e+00 4.589908e-07 2.453967e-04 2.350494e-05 1.846576e-01 1.519901e-10

# show estimate, standard error, t-statistic and two-sided p-value for main air pollution coefficient
round(summary(f1$glmfit)$coef[2,],4)
#  Estimate Std. Error    t value   Pr(>|t|) 
#    0.0870     0.1135     0.7668     0.4432 

# show number of non-missing observations both for met-only analysis and for met+pollution analysis 
# (these two numbers are the same for this particular dataset)

f1$nsub1
# [1] 4746
f1$nsub2
# [1] 4746

# plot results
plotTS(f1)

# The result of the plot can be seen by clicking on the following link:
# http://www.stat.unc.edu/faculty/rs/CApollution/PlotExample.pdf


# now do the same thing again but with nonlinear fit with 6 df and for the distributed lag model for ozone based
# on lags 0, 1, 2, 3

f2=fitTS(X8$alldeathnoac65,cbind(X8$tmax,X8$tmin,X8$MAXRH),0.001*X8$o3DMOL8N,0.02,0.15,0.001,0.075,c(0,1,2,3),7,6,6,6)

# plot results
plotTS(f2)

# The result of the plot can be seen by clicking on the following link:
# http://www.stat.unc.edu/faculty/rs/CApollution/PlotExample2.pdf