I am using a GAM in R to compare time series data from 3 countries. The data sets are of hourly measurements for one year. The main aim here is to show when at which time of the day and day of the year are the data well matched to the mean series. The script below shows my attempt:
## Ozone measurements for three countries in Europe ## Find similarities between time series require(plyr) require(lattice) require(mgcv) TopFolder <- list("http://www.nilu.no/projects/ccc/onlinedata/ozone/CZ03_2009.dat" ,"http://www.nilu.no/projects/ccc/onlinedata/ozone/CY02_2009.dat" ,"http://www.nilu.no/projects/ccc/onlinedata/ozone/BE35_2009.dat" ) ## create variable for the data data = ldply(TopFolder, header = TRUE, read.table, sep = "", skip = 3) ## define Ozone levels Ozone <- data$Value Ozone[Ozone==-999] <- NA Ozone <- data.frame(Ozone) ## define Datetime - need to concatenate arrays DateTime <- paste(data$Date,data$Hour, sep = " ") Date <- as.POSIXct(DateTime, format = "%d.%m.%Y %H:%M") ## define countries Countries <- c("Czech","Cyprus","Belgium") Country <- data.frame(Country = rep(Countries, each = 8760)) ## bind together Dat <- cbind(Ozone, Country = Country) Dat <- transform(Dat, Doy = as.numeric(format(Date,format = "%j")), Tod = as.numeric(format(Date,format = "%H")), DecTime = rep(seq(1,365, length = 8760),by = 3)) ## plot the Ozone data xyplot(Ozone~DecTime | Country, data = Dat, type = "l", col = 1, strip = function(bg = 'white',...)strip.default(bg = 'white',...)) ## generalised additive model mod1 <- gam(Ozone ~ Country + s(Doy,bs = "cc", k = 20) + s(Doy, by = Country, bs = "cc", k = 20) + s(Tod,bs = "cc", k=7) + s(Tod,by = Country,bs = "cr", k=7), data = Dat, method = "ML") plot(mod1,pages=1,scale=0,shade=TRUE) xyplot(resid(mod1) ~ Doy | Country, data = Dat, type = c("l","smooth")) mod2 <- gamm(Ozone ~ Country + s(Doy,bs = "cc", k = 20) + s(Doy, by = Country, bs = "cc", k = 20) + s(Tod,bs = "cc", k=7) + s(Tod,by = Country,bs = "cr", k=7), data = Dat, method = "ML",correlation = corAR1(form = ~ DecTime | Country))
One problem of this is that I am using an additive model with correlated errors i.e. in each time series of hourly measurements, a high concentration is followed by another high concentration. This can be addressed by adding a correlation matrix into the GAM:
However, R throws out an error as it is unable to get any more RAM from the OS. So, to tackle this I need to: (1) from mod1 - fit a time series model (mod3) to the residuals (2) obtain the acf(1) from mod3 (2) use the acf(1) to build a covariance matrix that we can then pass into mod2 thus simplifying the correlation matrix.
Any advice on how to to complete some or all of the three steps above would be highly appreciated. Although I think the first point that I an unsure of is how to extract the residuals for the different countries in mod1?