Tengo un modelo de efectos mixtos (de hecho, un modelo mixto aditivo generalizado) que me da predicciones para una serie temporal. Para contrarrestar la autocorrelación, utilizo un modelo corCAR1, dado que me faltan datos. Se supone que los datos me dan una carga total, por lo que necesito sumar todo el intervalo de predicción. Pero también debería obtener una estimación del error estándar en esa carga total.
Si todas las predicciones fueran independientes, esto podría resolverse fácilmente mediante:
con
El problema es que los valores pronosticados provienen de un modelo y los datos originales tienen autocorrelación. Todo el problema lleva a las siguientes preguntas:
- ¿Estoy en lo cierto al suponer que el SE en las predicciones calculadas puede interpretarse como la raíz de la varianza en el valor esperado de esa predicción? Tiendo a interpretar las predicciones como "predicciones medias" y, por lo tanto, sumar un conjunto completo de medias.
- ¿Cómo incorporo la autocorrelación en este problema, o puedo suponer con seguridad que no influirá demasiado en los resultados?
Este es un ejemplo en R. Mi conjunto de datos real tiene aproximadamente 34,000 mediciones, por lo que la escalabilidad es un problema. Esa es la razón por la que modelo la autocorrelación dentro de cada mes, de lo contrario, los cálculos ya no son posibles. No es la solución más correcta, pero la más correcta no es factible.
set.seed(12)
require(mgcv)
Data <- data.frame(
dates = seq(as.Date("2011-1-1"),as.Date("2011-12-31"),by="day")
)
Data <- within(Data,{
X <- abs(rnorm(nrow(Data),3))
Y <- 2*X + X^2 + scale(Data$dates)^2
month <- as.POSIXlt(dates)$mon+1
mday <- as.POSIXlt(dates)$mday
})
model <- gamm(Y~s(X)+s(as.numeric(dates)),correlation=corCAR1(form=~mday|month),data=Data)
preds <- predict(model$gam,se=T)
Total <- sum(preds$fit)
Editar:
Lección para aprender: primero revise todas las muestras en todos los archivos de ayuda antes de entrar en pánico. En los archivos de ayuda de predic.gam, puedo encontrar:
#########################################################
## now get variance of sum of predictions using lpmatrix
#########################################################
Xp <- predict(b,newd,type="lpmatrix")
## Xp %*% coef(b) yields vector of predictions
a <- rep(1,31)
Xs <- t(a) %*% Xp ## Xs %*% coef(b) gives sum of predictions
var.sum <- Xs %*% b$Vp %*% t(Xs)
Lo que parece estar cerca de lo que quiero hacer. Esto todavía no me dice exactamente cómo se hace. Podría llegar al hecho de que se basa en la matriz predictiva lineal. Cualquier idea aún es bienvenida.
fuente
Respuestas:
En notación matricial, un modelo mixto se puede representar como
y = X * beta + Z * u + épsilon
donde X y Z son matrices de diseño conocidas relacionadas con los efectos fijos y las observaciones de efectos aleatorios, respectivamente.
Aplicaría una transformación simple y adecuada (pero no la mejor) para corregir la autocorrelación que implica la pérdida de la primera observación, y reemplazar el vector de columna de [y1, y2, ... yn] por uno más pequeño por uno vector de columna de observación, a saber: [y2 - rho * y1, y3 - rho * y2, ..., yn - rho * y (n-1)], donde rho es su valor estimado para la autocorrelación en serie.
Esto se puede realizar multiplicando por una matriz T, formando T * y, donde la primera fila de T se compone de la siguiente manera: [-rho, 1, 0, 0, ....], la segunda fila: [0, -rho, 1, 0, 0, ...], etc. De manera similar, las otras matrices de diseño se cambian a T * X y T * Z. Además, la matriz de varianza-covarianza de los términos de error también se altera, ahora con términos de error independientes.
Ahora, solo calcule la solución con las nuevas matrices de diseño.
fuente