5
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ I would like to run mod2 in order to see how much RAM is actually needed. But if I define DecTime=(as.numeric(Date)-min(as.numeric(Date)))/3600 there is a Covariate must have unique values within groups for corAR1 objects error. Also, are you familiar with residuals and its method residuals.gam? $\endgroup$ – Roland Jul 26 '12 at 9:23
  • $\begingroup$ When trying to run mod2 R says 'Error: cannot allocate vector of size 530.2 Mb'. I am fairly familiar with residuals, e.g. I can plot the residuals for each country with 'xyplot(resid(mod1) ~ Doy | Country, data = Dat, type = c("l","smooth"))' but I don't know how I would go about fitting a model to the residuals of each country $\endgroup$ – KatyB Jul 26 '12 at 9:40
  • $\begingroup$ note from the plot of the residuals, there still appears to be some patterns present, hence the need for fitting the correlation matrix in gamm $\endgroup$ – KatyB Jul 26 '12 at 9:42
  • $\begingroup$ I know that you don't have enough memory, but maybe the 16 GB in total on my machine would suffice. But you didn't include how you defined DecTime in your code and my attempt to define it myself gave me an error that there are sometime more than one value per date/hour/country. You can subset the residuals to fit models per country, but I don't know if that would help you ultimately. $\endgroup$ – Roland Jul 26 '12 at 9:47
  • $\begingroup$ DecTime is now added to the data.frame. $\endgroup$ – KatyB Jul 26 '12 at 10:19
5
$\begingroup$

Coincidentally I have been pondering this problem for a little while and resorted to emailing Simon Wood about this only the other day.

His advice to me was (amongst other things but I don't think using bam()'s rho argument is applicable when you have by terms in the model) to check the example at the end of ?magic. That example I reproduce below with outputs:

## Now a correlated data example ... 
library(nlme)
## simulate truth
set.seed(1)
n <- 400
sig <- 2
x <- 0:(n-1)/(n-1)
f <- 0.2 * x^11 * (10 * (1-x))^6 + 10 * (10*x)^3 * (1-x)^10
## produce scaled covariance matrix for AR1 errors...
V <- corMatrix(Initialize(corAR1(.6),data.frame(x=x)))
## note here that V is what you would estimate from the model residuals
## but here we are generating the data from known correlation.
## You would plug your estimate of the correlation parameter in to corAR1(X)
##
## for the Cholesky factorisation of V
Cv <- chol(V)  # t(Cv)%*%Cv=V
## Simulate AR1 errors ...
e <- t(Cv) %*% rnorm(n, 0, sig) # so cov(e) = V * sig^2
## Observe truth + AR1 errors
y <- f + e 

## GAM ignoring correlation
b <- gam(y ~ s(x, k = 20))

## Fit smooth, taking account of *known* correlation...
## form a weight matrix
w <- solve(t(Cv)) # V^{-1} = w'w
## Use `gam' to set up model for fitting...
G <- gam(y ~ s(x, k = 20), fit=FALSE)
## fit using magic, with weight *matrix*
mgfit <- magic(G$y, G$X, G$sp, G$S, G$off, rank=G$rank, C=G$C, w=w)
## Modify previous gam object using new fit, for plotting...    
mg.stuff <- magic.post.proc(G$X, mgfit, w)
b2 <- b ## copy b
b2$edf <- mg.stuff$edf
b2$Vp <- mg.stuff$Vb
b2$coefficients <- mgfit$b 


## compare fits
layout(matrix(1:2, ncol = 2))
plot(b, main = "Ignoring correlation")
lines(x, f-mean(f), col=2)
plot(b2, main = "Known correlation")
lines(x, f-mean(f), col=2)
layout(1)

The resulting plot of the two fits looks like this:

enter image description here

You have the additional problem that your $V$ will not be a simple matrix like the one here. As you are fitting by Country your $V$ will be a block diagonal (I think that is the right term) matrix formed by first generating each $V_{country}$ and then stacking them in a diagonal fashion such that observations within a single Country are correlated with one another (so their entries in the matrix are none-zero, being $\rho^{|s|}$ where $\rho$ is the estiated AR(1) parameter for that Country and s is the separation in time points of the two residuals) but are uncorrelated with the observations from other countries (those entries in $V$ are zero). At least that is what the gamm() code you show would produce IIRC.

Although I haven't tried this, there are several options to form block diagonal matrices in R and you could use the function posted to R-Help some years ago to stick your individual $V_{country}$ together. Other options include converting each $V$ into a sparse matrix that the Matrix package understands, use its bdiag() function and then coerce back to a dense matrix using as.matrix(). See the help for the Matrix package for details.

As for obtaining the AR(1) terms for each countries residuals, fit your model via gam() but to help, place all the data in time order within country first before fitting. If you have any NA in the data, make sure to add na.action = na.exclude in the call to gam() as that will place NA back into the resdiuals that you extract via resid() from the model. You can split the vector of residuals by country like this:

res <- resid(mod)
res.spl <- split(res, Dat$Country)

Then you could simply sapply() the function acf() or pacf() to each component of res.spl to extract the lag 1 coefficient for each country. Those are the values to plug into corAR1() in the above example code.

Amazingly, your example code looks very similar to some that of my PhD student (although they are working on lakes) sent me the other day. Do you know someone working on hi-res lake data? If not perhaps they saw your post here?

$\endgroup$
  • $\begingroup$ many thanks for you very detailed answer. I used a similar script some time ago when I was investigating air temperature variations for a number of countries in Europe and I posted that script on r.789695.n4.nabble.com so there may be a link there. I cant say that I know of any lake scientsts, although some researches in this department are looking at NAO related processes in aquatic systems? $\endgroup$ – KatyB Jul 27 '12 at 13:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.