R: Ковариационная матрица для случайного эффекта в модели смешанных эффектов

Question

Согласно этот столбец, матрицы Omega и 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 является диагональной матрицей 3x3, то три числа в Variance должны быть элементами в диагонали Omega, аномер на Residual & Variance т.е. 1.453 должен быть моим sigma^2.

Вот другой пост также подтверждает это.

Однако,Omega Я использовал для генерации y сильно отличается от него.Вот мое истинное Omega:

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

И мое истинное sigma=1.Я не верю, что это связано с оценочной ошибкой, поскольку эти два числа так сильно различаются.

Тогда Согласно этой публикации, принятый ответ дал другой способ получить ковариационную матрицу Omega as

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

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

Но когда я использую этот метод для получения Omega, я получаю

> 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

Итак, pv - это массив с размерностью 3x3x20,т. е. каждая pv[,,i] - это та же матрица 3x3 с элементами, показанными выше.

Если я возьму любой pv[,,i] в качестве Omega, это не будет совпадать с Residuals в моем lmer результаты, и ни один из них не близок к моему истинному Omega.

Любая помощь с этим будет принята с благодарностью!

Редактировать: вот мой код для генерации данных и выполнения lmer примерка.

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