[R-sig-ME] results lme unstructured covariance matrix, again
John Fox
j|ox @end|ng |rom mcm@@ter@c@
Tue Jul 19 04:15:53 CEST 2022
Dear Ben,
On 2022-07-18 1:32 p.m., John Fox wrote:
> 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):
Rolf Turner pointed out to me that PDF attachments *are* acceptable, and
so I've attached a PDF with the LaTeX part of my message.
Best,
John
>
> 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
> Hamilton, Ontario, Canada
> 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/>
>>
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> 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/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mixed-model.pdf
Type: application/pdf
Size: 44332 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20220718/71e792e8/attachment-0001.pdf>
More information about the R-sig-mixed-models
mailing list