library(MASS,T)
library(splines)
library(chron)

###############################################################################
##  Author:  Bill Cox, USEPA,OAQPS,AQAD.  Dec 1, 2007
##   R code to model daily maximum 8-hour ozone as a function of daily 
##   meteorological variables plus terms to account for calendar day,
##   dow-of-week and a term "year" to account for trends over time.
##   This functions plots the raw and meteorologically adjusted trends
##   over the availabe years with data.  Assumes an R data file is present
##   that identifies the urban area along with the daily ozone and 
##   meteorological variables used in the prediction.  An object is saved
##   that contains the fitted object, the raw and meteorological adjusted
##   seasonal average ozone for each year, a vector of F-statistics for
##   each of the variables in the model and the R-square statistics for fit. 
##   R software is free and available at http://cran.r-project.org/
###############################################################################

###############################################################################
## Function to create circular splines
###############################################################################
bc <- function(x, period = 24., degree = 3., nknots = 6.)
 {
if(any(x < 0. | x >= period))
stop("x out of bound")
if(degree < 1.)
stop("degree must be a positive integer")
knots.x <- seq(0., period, length = nknots + 1.)
delta <- diff(knots.x[1.:2.])
knots.left <- seq(from =  - delta, by =  - delta, length = degree)
knots.right <- seq(from = period + delta, by = delta, length = degree)
knots <- c(knots.left, knots.x, knots.right)
X <- spline.des(x = x, knots = knots, ord = degree + 1.)$design
X[, 1.:degree] <- X[, 1.:degree] + X[, ncol(X) - (degree - 1.):0.]
X <- X[, 1.:(ncol(X) - degree)]
return(X)
 }

###############################################################################
##   Read in the testdata set.
##   A complete description of the variables can be found in the publication
##   "The effects of meteorology on ozone in urban areas and their use in
##   assessing ozone trends", by Louise Camalier, et.al.  Atmospheric
##   Environment 41 (2007) 7127-7137
###############################################################################

cc <- c("character",rep("numeric",10),"factor","Date","numeric",
        "numeric","factor","numeric","character")
testdata <- read.table(file="C:/Ddrive/Ozone2006/testdata.txt",
                       header=TRUE,colClasses=cc)
#  Note:  User must specify directory location containing "testdata.txt"
###############################################################################
## Start of the function that performs all analysis
###############################################################################

fitmod <- function(mnam,name,cval)
{
###############################################################################
##  Subset the data by selecting an urban areas (mnam), and the variables
##  to be used in the fit.  Create fitted objects using all of the 
##  predictor variables and one with only the year factor present (raw)
###############################################################################

  mnam <- paste(substitute(mnam))
  
  varlist <-  c("m8max","tmax","rhavgmid","wsavgam","wsavgpm",
                "dt925","devt850","trandis","trandir","rainhrs",
                "yrf","sdate","jday","year","dowf","reg","state")
 dat  <- na.omit(testdata[testdata$mnam==mnam & testdata$year >= 1997,varlist])
 levels(dat$dowf) <- c("sun","mon","tue","wed","thu","fri","sat")
 
 state <- unique(dat$state)
 reg   <- unique(dat$reg)

 form    <- formula(m8max ~ ns(tmax,3)    + ns(rhavgmid,3)
                          + ns(wsavgam,3) + ns(wsavgpm,3) 
                          + ns(dt925,3)   + ns(devt850,3) + ns(rainhrs,3)
                          + ns(trandis,3) + bc(trandir,period=360,nknots=4)
                          + ns(jday,3)    + dowf + yrf)
 forr    <- formula(m8max ~ yrf)
 
 fitm <- glm(form,data=dat,x=F,y=T,model=F,
              na.action=na.omit,quasi(link=log,variance=mu))
 fitr <- glm(forr,data=dat,x=F,y=T,
              na.action=na.omit,quasi(link=log,variance=mu))
              
###############################################################################
## Create object with the urban names and the R-square statistics
###############################################################################
  
  cor    <-  round(100*(cor(fitted(fitm),fitm$y)^2),1)
  corout <- data.frame(mnam=mnam,cor=cor)
  
###############################################################################
## function that uses the fitted object to calculate the meteorologicaly 
## adjusted and unadjusted seasonal mean 8-hour ozone for each year
###############################################################################

predval <- function(mnam,fitm,fitr)
 {
  const <- mean(fitm$y)
  tempf <- predict(fitm, type = "terms")
  tempr <- predict(fitr, type = "terms")
  pvf   <- const*exp(tempf[,"yrf"])
  pvr   <- const*exp(tempr[,"yrf"])
  dat2 <- cbind(year=dat$year,pvf=pvf,pvr=pvr)
  pboth <- aggregate(dat2[,c("pvf","pvr")],list(year=dat2[,"year"]),mean)
  pboth$year <- as.numeric(as.character(pboth$year))
  zz <- data.frame(mnam=mnam,pboth)
  zz
 }
pval <- predval(mnam,fitm,fitr)

###############################################################################
##  Tabulate the number of daily observations by year
##  and append to the metadjusted and unadjsuted means for each year 
###############################################################################

cnt  <- aggregate(dat$m8max,list(year=dat$year),function(x) sum(!is.na(x)))
pval$cnt <- cnt$x

###############################################################################
##  Use the dropterm function from the MASS package to determine the F value
##  associated with droping each term from the regression model.  Create object
##  Fval to be output later
###############################################################################

aovw  <- dropterm(fitm,test="F")
 nams <- rownames(aovw)
 nams <- nams[-1]
 t1 <- gsub(" .*","",nams)
 t2 <- gsub("\\(","",t1)
 t3 <- gsub(",","",t2)
 t4 <- gsub("ns","",t3)
 vnam <- gsub("bc","",t4)
 Fval   <- data.frame(mnam=mnam,vnam=vnam,Fw=round(aovw[,"F value"][-1]))

###############################################################################
## Plot trends lines of the adjusted and unadjusted trends
###############################################################################

 ylim <- c(cval - 50,cval + 30)
 plot(pval$year,pval$pvf, las = 1.,cex = 1.4,cex.axis=1.4,cex.lab=1.4,
  ylim = ylim, pch = " ", lwd = 1.5,type="n",
  xlab="Year",ylab="Seasonal Average Ozone (ppb)")
  lines(pval$year,pval$pvf,type="b",lty=1,lwd=3,pch=19)
  lines(pval$year,pval$pvr,type="b",lty=3,lwd=3,pch=21)
  legend(1998,28,xjust=0,cex=1,bty="o",pch=c(19,21),
     c("Adjusted for Weather","Unadjusted for Weather"),
      lty=c(1,3),lwd=c(2,2),y.intersp = 1,text.width=4,box.lwd=1.5)
  mtext(cex=1.6,name,outer=T,side=3,line=-3)
 
###############################################################################
## create and output a list with the R-square, predicted seasonal averages 
## (adjusted and unadjusted), the fitted objects, the Fval and
## vector of names of variables used in the fitting.
###############################################################################

out <- list(cor=corout,pval=pval,fitm=fitm,
            fitr=fitr,Fval=Fval,vnam=vnam)
out
}
###  End of the analysis function

###############################################################################
## Call the function to fit the model, plot the trends and output statistics
###############################################################################

cinctest <- fitmod("cinc","Cincinnati, OH 8-hr Ozone Trends",65)
cinctest    #  Prints the contents of the object to the screen

##  End of Computer Code