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.


 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
 (eCSV <- mvrnorm(1, mu = rep(0,15), Sigma = sd %*% as.matrix(corMatrixCS) %*% sd))
 yCSV <- ( dm %*% beta ) + eCSV

Plot of the simulated data (y values vs days). enter image description here

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

Output table:

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.

enter image description here

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


If you mean you get same estimates with different compound symmetric correlation matrices, that is because GLSE reduces to LSE under compound symmetric covariance structure.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.