Construir una matriz de covarianza para GAM

Question

Estoy usando un GAM en R para comparar datos de series temporales de 3 países. Los conjuntos de datos son de mediciones por hora durante un año. El objetivo principal aquí es mostrar cuándo y en qué momento del día y día del año los datos coinciden bien con la serie media. El script a continuación muestra mi intento:

## Mediciones de ozono para tres países en Europa
## Encontrar similitudes entre series temporales
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"
)

## crear variable para los datos
data = ldply(TopFolder, header = TRUE, read.table, sep = "", skip = 3)

## definir niveles de ozono
Ozone <- data$Value
Ozone[Ozone==-999] <- NA
Ozone <- data.frame(Ozone)

## definir FechaHora - necesario concatenar arrays
DateTime <- paste(data$Date,data$Hour, sep = " ")
Date <- as.POSIXct(DateTime, format = "%d.%m.%Y %H:%M")

## definir países
Countries <- c("Czech","Cyprus","Belgium")
Country <- data.frame(Country = rep(Countries, each = 8760))

## unir
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))

## plotear los datos de ozono
xyplot(Ozone~DecTime | Country, data = Dat, type  = "l", col = 1,
       strip = function(bg = 'white',...)strip.default(bg = 'white',...))

## modelo aditivo generalizado
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))

Uno de los problemas es que estoy usando un modelo aditivo con errores correlacionados es decir, en cada serie temporal de mediciones por hora, una alta concentración es seguida por otra alta concentración. Esto se puede abordar agregando una matriz de correlación al GAM:

Sin embargo, R muestra un error ya que no puede obtener más RAM del sistema operativo. Por lo tanto, para abordar esto necesito: (1) de mod1 - ajustar un modelo de series temporales (mod3) a los residuos (2) obtener el acf(1) de mod3 (2) usar el acf(1) para construir una matriz de covarianza que luego podemos pasar a mod2 simplificando así la matriz de correlación.

Cualquier consejo sobre cómo completar alguno o todos los tres pasos anteriores sería muy apreciado. Aunque creo que el primer punto del que no estoy seguro es cómo extraer los residuos para los diferentes países en mod1?

Answer 1

Casualmente he estado reflexionando sobre este problema durante un tiempo y recurrí a enviar un correo electrónico a Simon Wood sobre esto solo el otro día.

Su consejo para mí fue (entre otras cosas, pero no creo que usar el argumento rho de bam() sea aplicable cuando tienes términos by en el modelo) revisar el ejemplo al final de ?magic. Ese ejemplo lo reproduzco a continuación con salidas:

## Ahora un ejemplo de datos correlacionados ...
library(nlme)
## simular la verdad
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
## producir una matriz de covarianza escalada para errores AR1...
V <- corMatrix(Initialize(corAR1(.6),data.frame(x=x)))
## nota aquí que V es lo que estimarías de los residuos del modelo
## pero aquí estamos generando los datos de una correlación conocida.
## Insertarías tu estimación del parámetro de correlación en corAR1(X)
##
## para la factorización de Cholesky de V
Cv <- chol(V)  # t(Cv)%*%Cv=V
## Simular errores AR1 ...
e <- t(Cv) %*% rnorm(n, 0, sig) # entonces cov(e) = V * sig^2
## Observar verdad + errores AR1
y <- f + e 

## GAM ignorando correlación
b <- gam(y ~ s(x, k = 20))

## Ajustar suavizado, teniendo en cuenta la correlación *conocida*...
## formar una matriz de pesos
w <- solve(t(Cv)) # V^{-1} = w'w
## Usar 'gam' para configurar el modelo para el ajuste...
G <- gam(y ~ s(x, k = 20), fit=FALSE)
## ajustar usando magic, con peso *matriz*
mgfit <- magic(G$y, G$X, G$sp, G$S, G$off, rank=G$rank, C=G$C, w=w)
## Modificar objeto gam anterior usando nuevo ajuste, para trazar...    
mg.stuff <- magic.post.proc(G$X, mgfit, w)
b2 <- b ## copiar b
b2$edf <- mg.stuff$edf
b2$Vp <- mg.stuff$Vb
b2$coefficients <- mgfit$b 

## comparar ajustes
layout(matrix(1:2, ncol = 2))
plot(b, main = "Ignorando correlación")
lines(x, f-mean(f), col=2)
plot(b2, main = "Correlación conocida")
lines(x, f-mean(f), col=2)
layout(1)

La gráfica resultante de los dos ajustes se ve así:

Tienes el problema adicional de que tu $V$ no será una matriz simple como la de aquí. Como estás ajustando por Country, tu $V$ será una matriz diagonal por bloques (creo que es el término correcto) formada primero generando cada $V_{country}$ y luego apilándolos de manera diagonal de modo que las observaciones dentro de un solo Country estén correlacionadas entre sí (por lo que sus entradas en la matriz no son cero, siendo $\rho^{|s|}$ donde $\rho$ es el parámetro AR(1) estimado para ese Country y s es la separación en puntos temporales de los dos residuos) pero no están correlacionadas con las observaciones de otros países (esas entradas en $V$ son cero). Al menos es lo que el código de gamm() que muestras produciría, si mal no recuerdo.

Aunque no he probado esto, hay varias opciones para formar matrices diagonales por bloques en R y podrías usar la función publicada en R-Help algunos años atrás para unir tus $V_{country}$ individuales. Otras opciones incluyen convertir cada $V$ en una matriz dispersa que el paquete Matrix entiende, usar su función bdiag() y luego convertir de nuevo a una matriz densa usando as.matrix(). Consulta la ayuda del paquete Matrix para más detalles.

En cuanto a obtener los términos AR(1) para los residuos de cada país, ajusta tu modelo vía gam() pero para ayudar, coloca todos los datos en orden de tiempo dentro de cada país antes de ajustar. Si tienes algún NA en los datos, asegúrate de agregar na.action = na.exclude en la llamada a gam() ya que eso colocará los NA de vuelta en los residuos que extraes mediante resid() del modelo. Puedes dividir el vector de residuos por país de esta manera:

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

Luego podrías simplemente aplicar sapply() la función acf() o pacf() a cada componente de res.spl para extraer el coeficiente de retardo 1 para cada país. Esos son los valores para insertar en corAR1() en el código de ejemplo anterior.

Sorprendentemente, tu código de ejemplo se parece mucho a algunos que mi estudiante de doctorado (aunque están trabajando en lagos) me envió el otro día. ¿Conoces a alguien que esté trabajando en datos de lagos de alta resolución? Si no, ¿quizás vieron tu publicación aquí?