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

[James Pustejovsky]

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]]