Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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
- [,1] [,2] [,3]
- [1,] 0.6181442 0 0
- [2,] 0.0000000 0 0
- [3,] 0.0000000 0 0
- M1 <- lmer(ym ~ 0+XC1 + (0+X2 | subject))
- rr <- ranef(M1,condVar=TRUE)
- pv <- attr(rr[[1]],"postVar")
- str(pv)
- > 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
- 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))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement