For repeated observation analysis, what's the difference between fitting GLS with compound symmetry and random-intercept LMM with CS?

Question

I will refer to R and SAS examples, but this is a general statistical question.

I don't fully understand two things.

Let's assume I want to fit a model for observations repeated over time and compare if the means differ between men and women.

In SAS you can use the PROC MIXED procedure for it in 2 ways:

• using a mixed model, by employing random effects - "RANDOM" statement.
• using a Generalized Least Square linear model corrected for dependency and heteroscedasticity, using the "REPEATED" statement.

In R you can do the same using, for example:

• nlme::lme() for a mixed model
• nlme::gls() for the GLS estimated LM.

Now, in both, SAS and R cases I can define the residual covariance structure. Let's say it will be "compound symmetry".

And here I'm confused.

What does it mean to have, for example, a mixed model with "random intercept" AND select the CS (where does it apply? to random effects? to the residuals)?

from

Just having a GLS model with the CS (here I know it goes to the residuals).

What's the sense, the point to define both random effect relationships AND the residual covariance?

Exemplary code taken from the internet:

lme(weight ~ Time * Diet, BodyWeight,
                   random = ~1|Rat,
             corr=corCompSymm(form=~1|Rat))

It specifies both random intercept (for the Rat being the cluster, experimental subject) AND also the corCompSymm - CS for the residuals.

...and I even saw there are structures for the random effects itself, for example pdCompSymm(). Same in SAS - the covariance can be set for both residuals and random effects.

What's the point to specify both relationships between random effects and residuals?

Answer 1

To the best of my knowledge you're exactly right that these terms are redundant in this particular case. As an experimental demonstration (not "proof"):

First, lme is notoriously quiet about unidentifiable models. It fits this model without complaint:

library(nlme)
## lme with both Rat REs and cs
m1 <- lme(weight ~ Time * Diet, BodyWeight,
    random = ~1|Rat,
    corr=corCompSymm(form=~1|Rat))

But trying to compute intervals(m1) gives

Error in intervals.lme(m1) : cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance

and the compound-symmetric correlation parameter Rho is estimated as 1.2e-17, i.e. effectively zero (this is not a reliable indicator, however, as there are infinitely many combinations of Rho and among-rat variation that would give the same answer).

Fitting the two special-case models (rat RE without R-side correlation, corCompSymm in a gls model):

m2 <- lme(weight ~ Time * Diet, BodyWeight,
          random = ~1|Rat)
m3 <- gls(weight ~ Time * Diet, BodyWeight,
          corr=corCompSymm(form=~1|Rat))

Comparing the log-likelihoods of these models shows that they are all the same (up to floating-point precision: the relative difference is on the order of 1e-15) — again, not mathematical proof, but convincing evidence.

• The one case where it's important to distinguish between the two special cases (R-side comp symm and rat-level intercept RE) is that if there is negative compound-symmetric correlation among observations within a group, that can only be modeled by the R-side/gls model, not by the rat-level RE (unless you are using a computational framework that allows negative RE variances ...)
• Not all RE-level/correlation structures are redundant/jointly unidentifiable. For example, it could make sense to use an autoregressive model for the within-rat correlations (corAR1 in R) ...