La varianza de la suma de los valores de predicción de un modelo de efectos mixtos en un unicc

Question

Tengo un modelo de efectos mixtos (de hecho, un aditivo generalizado (modelo mixto) que me da predicciones para un unicc. Para contrarrestar la autocorrelación, yo uso un corCAR1 modelo, dado el hecho de que tengo los datos que faltan. Los datos se supone que me da un total de la carga, así que tengo que suma más de todo el intervalo de predicción. Pero también debo obtener una estimación del error estándar en que la carga total.

Si todas las predicciones sería independiente, esto puede ser fácilmente resuelto por :

$Var(\sum^{n}_{i=1}E[X_i]) = \sum^{n}_{i=1}Var(E[X_i])$ con $Var(E[X_i]) = E(E[X_i])^2$

El problema es, la predicción de los valores vienen de un modelo, y el original de los datos ha de autocorrelación. Todo el problema conduce a las siguientes preguntas :

1. Estoy en lo cierto al suponer que la SE en el cálculo de las predicciones puede ser interpretado como la raíz de la varianza en el valor esperado de esa predicción? Tiendo a interpretar las predicciones como "la media de las predicciones", y por lo tanto la suma de un conjunto de medios.
2. ¿Cómo puedo incorporar la autocorrelación en este problema, o se puede asumir con seguridad que no influyen en los resultados demasiado?

Este es un ejemplo en el R. Mi verdadero conjunto de datos tiene cerca de 34.000 mediciones, por lo que escalabilidad es un problema. Esa es la razón por la que el modelo de la autocorrelación dentro de cada mes, de lo contrario, los cálculos no son posibles más. 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 :

La lección a aprender : primero ir a través de todas las muestras en todos los archivos de ayuda antes de que el pánico. En los archivos de ayuda de predecir.gam, me pueden 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)

Que parece estar cerca de lo que quiero hacer. Con esto no me diga exactamente cómo se hace. Yo podría llegar tan lejos como el hecho de que se basa en el predictor lineal de la matriz. Cualquier conocimiento son bienvenidos.

Answer 1

En notación matricial de un modelo mixto puede ser representado como

y = X*beta + Z*u + epsilon

donde X y Z son conocidos diseño de matrices relacionadas con los efectos fijos y de efectos aleatorios, observaciones, respectivamente.

Me gustaría aplicar una sencilla y adecuada (pero no el mejor) de transformación para la corrección de auto-correlación, que implica la pérdida de la primera observación, y reemplazar el vector columna de [y1, y2,...yn] con un menor por parte de una observación vector columna, a saber: [y2 - rho*y1, y3 - rho*y2,..., yn - rho*y(n-1)], donde rho es su valor estimado para la serie de auto-correlación.

Esto se puede realizar mediante la multiplicación por una matriz T, formando T*y, donde la 1ª fila de T se compone como sigue: [ -rho, 1, 0, 0,....], la 2ª fila: [ 0, -rho, 1, 0, 0, ...], etc. Del mismo modo, el diseño de las matrices se cambian a T*X y T*Z. También, la varianza-covarianza de la matriz de los términos de error es alterado, ahora con independiente de los términos de error.

Ahora, acaba de calcular la solución con el nuevo diseño de las matrices.