I'm currently struggling with linear autocorrelated models. I have correctly simulated some datasets to understand how homoscedasticity and heteroscedasticity work and now I'm focused on correlated residuals. Therefore I do simulate a dataset with correlated residuals.
The problem: in the end I compare the estimates of the models that do not or do include a specification to model the correlation structure, expecting the latter one to obtain better results and/or estimates. However, all the estimates are the same (?)
My simulated dataset:
- 5 subjects
- 3 days
- correlation structure = compound symmetry. Rho value = 0.85 (I've also tried AR1, and happens the same).
- Variance function: the variance is depending on the mean. Therefore the model is heterocedastic. More specifically, the variance of the residuals = mean of the y values.
# DATA SIMULATION days <- gl(3,1,15) # 3 days. subject <- gl(5,3,15) # subjects dm <- model.matrix(~ 1+ days, data.frame(days)) # design matrix beta <- c(10,20,30) mu <- dm%*%beta # variance-covariance matrix CS <- corCompSymm(value = 0.85, form = ~ days|subject, fixed = FALSE) CS <- Initialize(CS, data =data.frame(days)) corMatrixCS<- bdiag(corMatrix(CS)) sd <- diag(sqrt(as.vector(mu))) # residuals set.seed(1234) (eCSV <- mvrnorm(1, mu = rep(0,15), Sigma = sd %*% as.matrix(corMatrixCS) %*% sd)) yCSV <- ( dm %*% beta ) + eCSV
Input arguments of the gls function.
mB <- gls( model = y ~ 1, data = dataCorrVar) # offset model mB2 <- gls( model = y ~ 1 + days, data = dataCorrVar) # No correlation or variance applied mB3 <- gls(model = y ~ 1 + days, data = dataCorrVar, correlation = corCompSymm(value = 0.85, form=~ days|subject) ) mB4 <- gls( model = y ~ 1 + days, data = dataCorrVar, weights = varPower(fixed = 1)) mB5 <- gls(model = y ~ 1 + days, data = dataCorrVar, correlation = corCompSymm(value = 0.85, form=~ days|subject), weights = varPower(fixed = 1))
1st model: offset model. Only the intercept is calculated.
2nd model: "raw model". No variance or correlation modeling.
3rd model: Only correlation is modeled.
4th model: Only variance is modeled.
5th model: Both variance and correlation structure are modeled.
As I said before, I expect the 4th model to give the best results. And it does in terms of AIC or BIC values. However, to me, it is extremely rare that all the estimates are exactly the same and only the standard errors are modified (std error between parenthesis).