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