Build a covariance matrix for GAM

Question

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?

Answer 1

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:

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?