• API
• FAQ
• Tools
• Archive
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.

Top