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

[Ben Bolker]

   For what it's worth you *can* fit this overparameterized model in 
lmer as well, you just have to override some controls

m4 <- lmer(y ~ 1 + t2 + (0 + t1 + t2 | person), data  = dat,
            control = lmerControl(check.nobs.vs.nRE = "ignore"))


preds <- cbind(sapply(list(m1, m2, m3), predict, level = 0),
       predict(m4, re.form = NA))

  The predictions are all the same **if you predict at the population 
level**; if you try to make predictions taking the random effects into 
account, the results differ because what is counted as "residual error" 
vs "random effects variation" differs among the models.

On 2022-07-18 12:05 PM, James Pustejovsky wrote:
> Hi Ben,
>
> I agree that this model is over-parameterized, so the distinction between
> the random effects variance components and the level-1 variance is
> arbitrary. You can see this by re-estimating the model with a fixed sigma
> that is less than or equal to the MLE. Here's an example, modified a bit
> from Dr. Fox's:
>
> library(nlme)
> library(mvtnorm)
> set.seed(123)
> N <- 1000
> t1 <- c(rep(1, N), rep(0, N))
> t2 <- 1 - t1
> person <- rep(1:N, 2)
> e <- rmvnorm(N, mean = c(0,1), sigma = matrix(c(1, 0.4, 0.4, 0.5), 2, 2))
> y <- as.vector(e)
> dat <- data.frame(y, t1, t2, person)
>
> m1 <- lme(y ~ 1 + t2, random = ~ 0 + t1 + t2 | person, data=dat)
> m2 <- lme(y ~ 1 + t2, random = ~ 1 + t2 | person, data=dat)
> m3 <- lme(y ~ 1 + t2, random = ~ 0 + t1 + t2 | person, data=dat,
>            control = lmeControl(sigma = 0.1))
>
> The marginal variance-covariance is identical across all three models:
> getVarCov(m1, type="marginal", individual="1")
> getVarCov(m2, type="marginal", individual="1")
> getVarCov(m3, type="marginal", individual="1")
>
> I would guess that the fact that lme() returns a result at all is due to
> how the model is estimated. lme() profiles out the sigma estimator and then
> solves for the MLEs of the random effects variance-covariance matrix. As
> Wolfgang noted, the result might depend on the optimizer. One peculiarity
> is that the log likelihood of m3 differs from that of the other two:
>> logLik(m1)
> 'log Lik.' -2301.013 (df=6)
>> logLik(m2)
> 'log Lik.' -2301.013 (df=6)
>> logLik(m3)
> 'log Lik.' -2323.984 (df=5)
>
> James
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models