2
$\begingroup$

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?

$\endgroup$

1 Answer 1

3
$\begingroup$

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) ...
$\endgroup$
2
  • $\begingroup$ Thank you very much, Sir! It clears my doubts entirely. $\endgroup$ Commented Nov 2, 2022 at 22:24
  • $\begingroup$ It is in general unidentifiable if we specify both the random effect and correlation structure with the same grouping. I actually tried to estimate a model with both autocorrelation in the residual variance and random effect in the intercept grouped by the subject, using another set of data with repeated measurements. Same situation happened: the estimated standard deviation of the random effect was almost zero, and/or the correlation coefficient at its boundary (-1 or 1). $\endgroup$
    – DrJerryTAO
    Commented Sep 24, 2023 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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