混合效应模型中lmer估计Omega与真实值不符的问题求助

我在拟合混合效应模型时遇到了一个棘手的问题：我原本认为可以从lmer的输出结果中获取Omega矩阵（3×3对角矩阵）和sigma，但得到的估计值和我生成数据时用的真实值差异极大，这绝对不是普通的估计误差能解释的。

当前拟合结果与真实值对比 #

我的lmer拟合结果如下：

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

按照之前看到的帖子，Omega作为3×3对角矩阵，其对角元素应该对应上面subject组下的三个Variance值，而Residual行的Variance值1.453就是sigma²，另一帖子也验证了这个逻辑。但我生成数据时用的真实参数是：

Omega
[,1] [,2] [,3]
[1,] 0.6181442 0 0
[2,] 0.0000000 0 0
[3,] 0.0000000 0 0
真实sigma=1，两者的差异完全无法用估计误差来解释。

尝试的另一种提取方法及异常结果 #

之后我按照另一帖子的方法尝试获取Omega，代码如下：

M1 <- lmer(ym ~ 0+XC1 + (0+X2 | subject))
rr <- ranef(M1,condVar=TRUE)
pv <- attr(rr[[1]],"postVar")
str(pv)

得到的pv是一个维度为3×3×20的数组，每个pv[,,i]都是类似下面的矩阵：

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
这个结果既和lmer输出的随机效应方差不符，也和真实的Omega矩阵相差甚远。恳请各位大佬帮忙分析一下问题出在哪里！

附：数据生成与拟合的完整代码 #

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

内容的提问来源于stack exchange，提问作者Matata