Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I have a mixed effect model (in fact a generalized additive mixed model) that gives me predictions for a timeseries. To counter the autocorrelation, I use a corCAR1 model, given the fact I have missing data. The data is supposed to give me a total load, so I need to sum over the whole prediction interval. But I should also get an estimate of the standard error on that total load.

If all predictions would be independent, this could be easily solved by :

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

Problem is, the predicted values are coming from a model, and the original data has autocorrelation. The whole problem leads to following questions :

  1. Am I correct in assuming that the SE on the calculated predictions can be interpreted as the root of the variance on the expected value of that prediction? I tend to interprete the predictions as "mean predictions", and hence sum a whole set of means.
  2. How do I incorporate the autocorrelation in this problem, or can I safely assume that it wouldn't influence the results too much?

This is an example in R. My real dataset has about 34.000 measurements, so scaleability is a problem. That's the reason why I model the autocorrelation within each month, otherwise the calculations aren't possible any more. It's not the most correct solution, but the most correct one isn't feasible.

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)

Edit :

Lesson to learn : first go through all the samples in all the help files before panicking. In the help files of predict.gam, I can find :

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

Which seems to be close to what I want to do. This still doesn't tell me exactly how it is done. I could get as far as the fact that it's based on the linear predictor matrix. Any insights are still welcome.

share|improve this question
6  
Im not sure what the r program is doing but we have $$var(\sum_iE[X_i])=a^Tvar(E[X])a$$ Where $a$ is a column vector of ones and $var(E[X])$ is the covariance matrix for $E[X]=(E[X_1],\dots,E[X_n])^T$. Does this help? –  probabilityislogic Apr 4 '12 at 22:36
    
@probabilityislogic That's basically what the r program is doing. Thx for the math –  Joris Meys Apr 27 '12 at 8:16
1  
@probabilityislogic If you can wrap that into an answer, you can grab my +50 bounty. ;) –  e-sushi Nov 11 '13 at 12:08
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.