daily pastebin goal
62%
SHARE
TWEET

Untitled

a guest Feb 19th, 2019 63 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. Random effects:
  2. Groups   Name Variance  Std.Dev. Corr      
  3. subject  X21  8.558e+00 2.925380            
  4.          X22  2.117e-03 0.046011 -1.00      
  5.          X23  2.532e-05 0.005032  1.00 -1.00
  6. Residual      1.453e+00 1.205402            
  7. Number of obs: 100, groups:  subject, 20
  8.      
  9. > Omega
  10.           [,1] [,2] [,3]
  11. [1,] 0.6181442    0    0
  12. [2,] 0.0000000    0    0
  13. [3,] 0.0000000    0    0
  14.      
  15. M1 <- lmer(ym ~ 0+XC1 + (0+X2 | subject))
  16. rr <- ranef(M1,condVar=TRUE)
  17.  
  18. pv <- attr(rr[[1]],"postVar")
  19. str(pv)
  20.      
  21. > pv[,,1]
  22.               [,1]          [,2]          [,3]
  23. [1,]  0.2913922395 -4.588735e-03  5.017337e-04
  24. [2,] -0.0045887347  7.523883e-05 -8.101042e-06
  25. [3,]  0.0005017337 -8.101042e-06  8.773362e-07
  26.      
  27. N=20          #number of subject
  28. n=5           #number of observations per subject
  29.  
  30. sigma2 =1
  31.  
  32. #equally space point not include the start point
  33. time <- function(from, to, length.out) {
  34.     length.out <- length.out + 1
  35.     result <- seq(from, to, length.out = length.out)
  36.     result <- result[-1]
  37.     return(result)
  38. }
  39.  
  40. subject = matrix(0,nrow=N*n,ncol=1)
  41. for(i in 1:N){
  42.     for(j in (n*(i-1)+1):(n*i)){
  43.         subject[j]=i
  44.     }
  45. }
  46.  
  47.  
  48. X = array(0, dim = c(N, n, 3))#each X[i,,] is a nx3 matrix
  49. for (i in 1:N){
  50.     for (j in 1:n){
  51.     X[i,j,] <-c(1,time(0,10,n)[j],(time(0,10,n)[j])^2)
  52.     }
  53. }
  54.  
  55.  
  56. y = array(0, dim = c(N, n, 1))
  57.  
  58.  
  59. Omega <- matrix(0,nrow=3,ncol=3)
  60. Omega[1,1] = runif(1,0.01,1.01)#only omega1^2 is not equal to 0
  61.  
  62. beta <-rep(0,5)
  63. beta[1]= rnorm(1,mean=0.01,sd=1) #mu0
  64. beta[2]= rnorm(1,mean=0.005,sd=1) #mu1
  65. beta[3]= rnorm(1,mean=0.0025,sd=1) #mu2
  66.  
  67.  
  68. C1 = array(0, dim = c(N, 3, 5))
  69. for(i in 1:N){
  70.     C1[i,1,1]=C1[i,2,2]=C1[i,3,3]=1
  71. }
  72.  
  73.  
  74. muy = array(0, dim = c(N, n, 1))     #store the expextation of each yi
  75. Cov = array(0, dim = c(N, n, n))     #store the covariance matrix of y
  76.  
  77. for (i in 1:N){
  78.     muy[i,,] <- X[i,,]%*%C1[i,,]%*%beta
  79.     Cov[i,,] <- X[i,,]%*%Omega%*%t(X[i,,])+ sigma2*diag(n)
  80.     y[i,,] <- mvrnorm(n = 1, muy[i,,], Cov[i,,])
  81. }
  82.  
  83. ym <- as.vector(y[,,1])
  84.  
  85.  
  86. #change X into X2, which is in a matrix format, easy for compitation later
  87. X2 <- rbind(X[1,,],X[2,,])
  88. for(i in 2:(N-1)){
  89.     X2 = rbind(X2,X[i+1,,])
  90. }
  91.  
  92.  
  93. XC1=matrix(0,nrow=N*n,ncol=5)
  94. for(i in 1:N){
  95. XC1[(n*(i-1)+1):(i*n),]=X[i,,]%*%C1[i,,]
  96. }
  97.  
  98.  
  99.  
  100. M1<-lmer(ym ~ 0+XC1+(0+X2|subject))
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top