Eu tenho um modelo de efeito misto (de fato, um modelo misto aditivo generalizado) que me fornece previsões para séries temporais. Para combater a autocorrelação, eu uso um modelo corCAR1, pois tenho dados ausentes. Os dados devem fornecer uma carga total, portanto, preciso somar todo o intervalo de previsão. Mas também devo obter uma estimativa do erro padrão nessa carga total.
Se todas as previsões fossem independentes, isso poderia ser facilmente resolvido por:
com
O problema é que os valores previstos são provenientes de um modelo e os dados originais têm autocorrelação. Todo o problema leva às seguintes perguntas:
- Estou correto ao supor que o SE nas previsões calculadas pode ser interpretado como a raiz da variação no valor esperado dessa previsão? Costumo interpretar as previsões como "previsões médias" e, portanto, somar um conjunto inteiro de meios.
- Como incorporar a autocorrelação nesse problema ou posso assumir com segurança que isso não influenciaria muito os resultados?
Este é um exemplo em R. Meu conjunto de dados real tem cerca de 34.000 medições, então a escalabilidade é um problema. Essa é a razão pela qual eu modelo a autocorrelação dentro de cada mês, caso contrário, os cálculos não serão mais possíveis. Não é a solução mais correta, mas a mais correta não é viável.
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:
Lição a aprender: primeiro passe por todas as amostras em todos os arquivos de ajuda antes de entrar em pânico. Nos arquivos de ajuda do predict.gam, posso 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)
O que parece estar próximo do que eu quero fazer. Isso ainda não me diz exatamente como é feito. Eu poderia chegar até o fato de que é baseado na matriz preditora linear. Quaisquer insights ainda são bem-vindos.
fonte
Respostas:
Na notação matricial, um modelo misto pode ser representado como
y = X * beta + Z * u + epsilon
onde X e Z são matrizes de projeto conhecidas relacionadas aos efeitos fixos e observações de efeitos aleatórios, respectivamente.
Eu aplicaria uma transformação simples e adequada (mas não a melhor) para corrigir a correlação automática que envolve a perda da primeira observação e substituir o vetor da coluna [y1, y2, ... yn] por um menor por um vetor da coluna de observação, a saber: [y2 - rho * y1, y3 - rho * y2, ..., yn - rho * y (n-1)], em que rho é o valor estimado para a correlação automática serial.
Isso pode ser realizado multiplicando por uma matriz T, formando T * y, em que a primeira linha de T é composta da seguinte forma: [-rho, 1, 0, 0, ....], a segunda linha: [0, -rho, 1, 0, 0, ...] etc. Da mesma forma, as outras matrizes de design são alteradas para T * X e T * Z. Além disso, a matriz de variância-covariância dos termos de erro também é alterada, agora com termos de erro independentes.
Agora, apenas calcule a solução com as novas matrizes de design.
fonte