Covariance matrix for the random effect in mixed effects model

Question

According to this post, matrix Omega and sigma are in the results of lmer when we fitting the mixed effect model. And here is my result.

Random effects:
Groups   Name Variance  Std.Dev. Corr       
subject  X21  8.558e+00 2.925380            
         X22  2.117e-03 0.046011 -1.00      
         X23  2.532e-05 0.005032  1.00 -1.00
Residual      1.453e+00 1.205402            
Number of obs: 100, groups:  subject, 20

Since my Omega is a 3x3 diagonal matrix, so the three numbers in the Variance should be the elements in the diagonal of Omega and the number on the Residual&Variance i.e. 1.453 should be my sigma^2.

Here is another post also confirms this.

However, the Omega I used to generate y is quite different from it. Here is my true Omega:

> Omega
          [,1] [,2] [,3]
[1,] 0.6181442    0    0
[2,] 0.0000000    0    0
[3,] 0.0000000    0    0

And my true sigma=1. I do not believe the difference is due to the estimated error since these two numbers differ so much.

Then According to this post, the accepted answer gave another way to get the covariance matrix Omega as

M1 <- lmer(ym ~ 0+XC1 + (0+X2 | subject))
rr <- ranef(M1,condVar=TRUE)

pv <- attr(rr[[1]],"postVar")
str(pv)

But when I use this method to get my Omega, I get

> pv[,,1]
              [,1]          [,2]          [,3]
[1,]  0.2913922395 -4.588735e-03  5.017337e-04
[2,] -0.0045887347  7.523883e-05 -8.101042e-06
[3,]  0.0005017337 -8.101042e-06  8.773362e-07

So, pv is an array with dimension 3x3x20, i.e. each pv[,,i] is the same 3x3 matrix with elements shows above.

If I take any pv[,,i] as my Omega, this will not the same as the Residuals in my lmer results, and none of these close to my true Omega.

Any help with this will be greatly appreciated!

Edit: here is my code to generate data and do the lmer fitting.

N=20          #number of subject
n=5           #number of observations per subject

sigma2 =1

#equally space point not include the start point
time <- function(from, to, length.out) {
    length.out <- length.out + 1
    result <- seq(from, to, length.out = length.out)
    result <- result[-1]
    return(result)
}

subject = matrix(0,nrow=N*n,ncol=1)
for(i in 1:N){
    for(j in (n*(i-1)+1):(n*i)){
        subject[j]=i
    }
}


X = array(0, dim = c(N, n, 3))#each X[i,,] is a nx3 matrix
for (i in 1:N){
    for (j in 1:n){
    X[i,j,] <-c(1,time(0,10,n)[j],(time(0,10,n)[j])^2)
    }
} 


y = array(0, dim = c(N, n, 1))


Omega <- matrix(0,nrow=3,ncol=3)
Omega[1,1] = runif(1,0.01,1.01)#only omega1^2 is not equal to 0

beta <-rep(0,5)
beta[1]= rnorm(1,mean=0.01,sd=1) #mu0
beta[2]= rnorm(1,mean=0.005,sd=1) #mu1
beta[3]= rnorm(1,mean=0.0025,sd=1) #mu2


C1 = array(0, dim = c(N, 3, 5))
for(i in 1:N){
    C1[i,1,1]=C1[i,2,2]=C1[i,3,3]=1
}


muy = array(0, dim = c(N, n, 1))     #store the expextation of each yi
Cov = array(0, dim = c(N, n, n))     #store the covariance matrix of y

for (i in 1:N){
    muy[i,,] <- X[i,,]%*%C1[i,,]%*%beta
    Cov[i,,] <- X[i,,]%*%Omega%*%t(X[i,,])+ sigma2*diag(n)
    y[i,,] <- mvrnorm(n = 1, muy[i,,], Cov[i,,])
}

ym <- as.vector(y[,,1])


#change X into X2, which is in a matrix format, easy for compitation later
X2 <- rbind(X[1,,],X[2,,])
for(i in 2:(N-1)){
    X2 = rbind(X2,X[i+1,,])
}


XC1=matrix(0,nrow=N*n,ncol=5)
for(i in 1:N){
XC1[(n*(i-1)+1):(i*n),]=X[i,,]%*%C1[i,,]
}



M1<-lmer(ym ~ 0+XC1+(0+X2|subject))

Edit: According to the posted answer, I know why the VarCorr() and ranef() gave different solutions. But still, I am confused since VarCorr() should be the estimate of my $\Omega$ then, but they are quite different.

Answer 1

What you get from ranef() is not an estimate of the covariance matrix $\Omega$. For a more detailed explanation of the differences, check this post.

Moreover, you seem to be fitting a rather overly complicated model for the data you have simulated. Note that the optimization algorithms under the postulated model work under the premise that the covariance matrix $\Omega$ needs to be positive definite, meaning that the variances need to be positive. When the variance is actually zero, the shape of the log-likelihood can be rather unusual near leading to problems during the optimization. In your case, you have two variances that are zero making it a hard problem.

Finally, for another way to simulate from a mixed model check the code provided in this example - you will need to change the last line and use rnorm() instead of rbinom(), and specify as sd the $\sigma$ in your formulation.

Answer 2

I am confused since VarCorr() should be the estimate of my $\Omega$ then, but they are quite different.

There seems to be a typo at ym.

    ym <- as.vector(y[,,1])

    # compared with:
    # ym <- as.vector(t(y[,,1]))

In addition, if the Omega is zeros everywhere except at Omega[1,1], only the first column of X has random effect. As shown in the equivalent simulation below, only alpha[1] is non-zero. Thus it is unlikely to see meaningful random effects from columns 2 and 3 of X.

for (i in 1:N){
    muy[i,,] <- X[i,,]%*%C1[i,,]%*%beta
    Cov[i,,] <- X[i,,]%*%Omega%*%t(X[i,,])+ sigma2*diag(n)
    # y[i,,] <- mvrnorm(n = 1, muy[i,,], Cov[i,,])
    
    # is equivalent to:
    alpha = mvrnorm(n = 1, rep(0,n), Omega)
    y[i,,] <- muy[i,,] +
              rnorm(n,sd=sqrt(sigma2)) +
              alpha[1] * X[i,,1] +
              alpha[2] * X[i,,2] +
              alpha[3] * X[i,,3] 
}