# [R-sig-ME] results lme unstructured covariance matrix, again

John Fox j|ox @end|ng |rom mcm@@ter@c@
Mon Jul 18 19:32:04 CEST 2022

Dear Ben,

My apologies---I paid attention to the inessential part of your message,
which was the parametrization of the model rather than the
variance-covariance structure.

This is a little awkward because of the LaTeX (and as far as I know, I
can't attach a PDF, but maybe you can paste this into a LaTeX editor):

The model is $Y_{ij} = \alpha + \beta x_{ij} + \delta_{1i} + \delta_{2i} x_{ij} + \varepsilon_{ij}$ for $i = 1, \ldots, n$ and $j = 1,2$, where
$x_{i1} = 1$ and $x_{i2} = 0$.

The variances of the random effects are $V(\delta_{1i}) = \psi_1^2$ and
$V(\delta_{2i}) = \psi_2^2$, and their covariance is $C(\delta_{1i}, \delta_{2i}) = \psi_{12}$. The observation-level error variance is
$V(\varepsilon_{ij}) = \sigma^2$.

Then the composite error is $\zeta_{ij} = \delta_{1i} + \delta_{2i} x_{ij} + \varepsilon_{ij}$ with variance $V(\zeta_{ij}) = \psi_1^2 + x^2_{ij} \psi_2^2 + 2 x_{ij} \psi_{12} + \sigma^2$, which is
$V(\zeta_{i1}) = \psi_1^2 + \psi_2^2 + 2\psi_{12} + \sigma^2$ for $j = 1$ and  $V(\zeta_{i2}) = \psi_1^2 + \sigma^2$ for $j = 2$, and
covariance $C(\zeta_{i1}, \zeta_{i2}) = \psi_1^2 + \psi_2^2$.

There are, as you say, 4 variance-covariance components, $\psi_1^2, \psi_2^2, \psi_{12}, \sigma^2$, and just 3 variances and covariances
among the composite disturbances, $V(\zeta_{i1}), V(\zeta_{i2}), C(\zeta_{i1}, \zeta_{i2})$, and so---as you and several others have
noted (but I missed)---the variance-covariance components are
underidentified.

lmer(), but not lme(), appears to use a heuristic, counting the number
of random effects and observations, to detect underidentification. I
don't know whether that's bullet-proof, but it works in this case. Maybe
Ben Bolker, who has already responded, can comment further.

Best,
John

John Fox, Professor Emeritus
McMaster University
web: https://socialsciences.mcmaster.ca/jfox/

On 2022-07-18 6:39 a.m., ben pelzer wrote:
> Dear John,
>
> Thank you for answering my question and your nice example with the
> alternative model formulations!!. But still I'm in doubt about the
> results. The fact that there are only three observed (co)variances which
> are being reproduced by four parameters estimates of lme leaves me
> confused. Actually, I would think that there is no unique solution
> without some constraint being applied. But I could not find something
> about such a constraint in lme documentation. The random effects of t1
> and t2, or of the intercept and t2 in your alternative model, should be
> sufficient to reproduce the observed (co)variances, so the residual
> variance is "unnecassary". I guess that is the reason that lmer is not
> able to estimate the model.
>
> library(lme4)
> m3 <- lmer(y ~ 1+t2+(1+t2|person), data=da)
>
> Error: number of observations (=20) <= number of random effects (=20)
> for term (1 + t2 | person); the random-effects parameters and the
> residual variance (or scale parameter) are probably unidentifiable
>
>
> It's interesting to notice that the two model specifications (t1+t2 and
> 1+t2) lead to different residual variances estimated by lme. What do
> these residual variances mean? And of course also: what do the random
> effect variances estimated by both model formulations actually mean or
> stand for? Do you have any ideas? Kind regards,
>
> Ben.
>
> On Sat, 16 Jul 2022 at 16:50, John Fox <jfox using mcmaster.ca
> <mailto:jfox using mcmaster.ca>> wrote:
>
>     Dear Ben,
>
>     First, I'll make this into a reproducible example:
>
>       > set.seed(123)
>       > t1 <- c(rep(1, 10), rep(0, 10))
>       > t2 <- 1 - t1
>       > person <- rep(1:10, 2)
>       > y <- t2 + rnorm(20)
>       > da <- data.frame(y, t1, t2, person)
>
>       > library(nlme)
>
>     Then note that the random-effect specification 0 + t1 + t2 is simply a
>     reparametrization of 1 + t2 (i.e., 1 = t1 + t2), which produces the
>     same
>     fit to the data (same fixed effects, same restricted log-likelihood):
>
>       > m1 <- lme(y ~ 1 + t2, random = ~ 0 + t1 + t2 | person, data=da)
>       > m2 <- lme(y ~ 1 + t2, random = ~ 1 + t2 | person, data=da)
>       > m1
>     Linear mixed-effects model fit by REML
>         Data: da
>         Log-restricted-likelihood: -25.92726
>         Fixed: y ~ 1 + t2
>     (Intercept)          t2
>        0.07462564  1.13399632
>
>     Random effects:
>        Formula: ~0 + t1 + t2 | person
>        Structure: General positive-definite, Log-Cholesky parametrization
>                StdDev    Corr
>     t1       0.8964136 t1
>     t2       0.9856215 0.647
>     Residual 0.3258015
>
>     Number of Observations: 20
>     Number of Groups: 10
>
>       > m2
>     Linear mixed-effects model fit by REML
>         Data: da
>         Log-restricted-likelihood: -25.92726
>         Fixed: y ~ 1 + t2
>     (Intercept)          t2
>        0.07462564  1.13399632
>
>     Random effects:
>        Formula: ~1 + t2 | person
>        Structure: General positive-definite, Log-Cholesky parametrization
>                   StdDev    Corr
>     (Intercept) 0.8787887 (Intr)
>     t2          0.7540826 -0.302
>     Residual    0.3707215
>
>     Number of Observations: 20
>     Number of Groups: 10
>
>     Finally, it's unnecessary to supply the intercept 1 in the model
>     formula
>     since the intercept is implied if it's not explicitly excluded:
>
>       > m3 <- lme(y ~ t2, random = ~ t2 | person, data=da)
>       > m3
>     Linear mixed-effects model fit by REML
>         Data: da
>         Log-restricted-likelihood: -25.92726
>         Fixed: y ~ t2
>     (Intercept)          t2
>        0.07462564  1.13399632
>
>     Random effects:
>        Formula: ~t2 | person
>        Structure: General positive-definite, Log-Cholesky parametrization
>                   StdDev    Corr
>     (Intercept) 0.8787887 (Intr)
>     t2          0.7540826 -0.302
>     Residual    0.3707215
>
>     Number of Observations: 20
>     Number of Groups: 10
>
>     I hope this helps,
>        John
>
>     On 2022-07-16 10:05 a.m., ben pelzer wrote:
>      > Sorry, my previous mailed contained another question which is
>     irrelevant...
>      > I deleted that now.
>      >
>      >
>      > Hi all,
>      >
>      > I have a question about results from lme of package nlme.
>      >
>      > Suppose the data consists of repeated measures at two fixed time
>     points.
>      >
>      > I used the following equation:
>      >
>      >
>      >
>      > Model1 <- lme ( y ~ 1+t2 , random = ~ 0 + t1+t2|person, data=da)
>      >
>      >
>      >
>      > y is the dependent, t1 and t2 are binary dummy variables, valued
>     0 or 1,
>      > indicating the time point.  Model1 is estimated without any
>     convergence
>      > problems and the reproduced (co)variances found with
>      >
>      >
>      >
>      > getVarCov(Model1, type=”marginal”, indivual=”1”)
>      >
>      >
>      >
>      > are identical to the observed (co)variances.
>      >
>      >
>      > My question is:  how can lme estimate 4 (co)variances with only 3
>     known
>      > (co)variances?
>      >
>      >
>      >
>      > The 4 estimates concern:
>      >
>      > -          std. deviation of the random effect of dummy t1
>      >
>      > -          std. deviation of the random effect of dummy t2
>      >
>      > -          covariance of the random effects of the dummies t1 and
>     t2 t1
>      >
>      > -          residual std. error
>      >
>      >
>      >
>      > Related to the question above: how can the variances of the
>     random effects
>      > and the residual std. error be interpreted?
>      >
>      >
>      >
>      > Thanks for any help,
>      >
>      >
>      >
>      > Ben.
>      >
>      >       [[alternative HTML version deleted]]
>      >
>      > _______________________________________________
>      > R-sig-mixed-models using r-project.org
>     <mailto:R-sig-mixed-models using r-project.org> mailing list
>      > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>     <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>     --
>     John Fox, Professor Emeritus
>     McMaster University
>     Hamilton, Ontario, Canada
>     web: https://socialsciences.mcmaster.ca/jfox/
>     <https://socialsciences.mcmaster.ca/jfox/>
>