I am running simulation on bivariate longitudinal data.
I will describe my code which I hope to reduce variation on csse for group 0.
So, there is tow classes Y= (0,1) and for each class I generate linear mixed effects model. There are 4 random parameters (intercept and slopes ) for each group, that is mean the covariance matrix is 4x4 and there are correlated between random effects.
Also, the residuals are matrix 2x2 (each variable has error) and there is a correlation . I assume that there are 4 fixed effects parameters for each variable in each group.
Therefore, my code is:
###~~~ simulate two correlated responses variables in longitudinal data ~~~###
## set means of random effects (intercepts and slopes)
m.0 = c(8,8,4,5)
m.1 = c(7,6,3,0)
## correlation matrix of random effects
cor.0 = matrix(c(1.00,0.20,0.55,0.29,
0.20,1.00,0.13,0.18,
0.55,0.13,1.00,0.65,
0.29,0.18,0.65,1.00),nrow=4)
cor.1 = matrix(c(1.00,0.30,0.40,0.56,
0.30,1.00,0.51,0.65,
0.40,0.51,1.00,0.81,
0.56,0.65,0.81,1.00),nrow=4)
### set correlation matrix of residuals
cor.R0 = matrix(c(1.00,0.15,
0.15,1.00),nrow=2)
cor.R1 = matrix(c(1.00,0.42,
0.42,1.00),nrow=2)
## generate covariance matrices
set.seed(1)
sds <- rnorm(4)^2
S0 = cor.0 * sds * rep(sds, each = nrow(cor.0))
sds <- rnorm(4)^2
S1 = cor.1 * sds * rep(sds, each = nrow(cor.1))
sds <- rnorm(2)^2
R0 = cor.R0 * sds * rep(sds, each = nrow(cor.R0))
sds <- rnorm(2)^2
R1 = cor.R1 * sds * rep(sds, each = nrow(cor.R1))
n = 200
## generate intercepte and slpos
library(mixAK)
B0 = rMVN(n,m.0,S0)$x
B1 = rMVN(n,m.1,S1)$x
E0 = rMVN(n,c(0,0),R0)$x
E1 = rMVN(n,c(0,0),R1)$x
Time = rep(c(0,3,6,9),times=n/2)
B0Time = B0[,c(2,4)]*Time[1:n]
B1Time = B1[,c(2,4)]*Time[1:n]
id = rep(1:n/2,each=4)
## set fix intercepts for each variable and groups
B0Fix1 = rep(0.5,n*2) #cs0
B1Fix1 = rep(40,n*2) #cs1
B0Fix2 = rep(40,n*2) #va0
B1Fix2 = rep(20,n*2) #va1
## get the equation
Y1.0 = B0Fix1 + B0[,1] + B0Time[,1] + E0[,1]
Y2.0 = B0Fix2 + B0[,3] + B0Time[,2] + E0[,2]
Y1.1 = B1Fix1 + B1[,1] + B1Time[,1] + E1[,1]
Y2.1 = B1Fix2 + B1[,3] + B1Time[,2] + E1[,2]
csse = c(Y1.0,Y1.1)
vase = c(Y2.0,Y2.1)
Y = as.factor(c(rep(0,400),rep(1,400)))
data.Sim = data.frame(id,Time,csse,vase ,Y)
The summaries are:
summary(data.Sim)
summary(data.Sim[data.Sim$Y==1,])
summary(data.Sim[data.Sim$Y==0,])
> summary(data.Sim[data.Sim$Y==1,])
id Time csse vase Y
Min. : 50.50 Min. :0.00 Min. : 15.65 Min. :-47.71 0: 0
1st Qu.: 62.88 1st Qu.:2.25 1st Qu.: 55.83 1st Qu.: 11.38 1:400
Median : 75.25 Median :4.50 Median : 73.44 Median : 23.33
Mean : 75.25 Mean :4.50 Mean : 76.43 Mean : 21.97
3rd Qu.: 87.62 3rd Qu.:6.75 3rd Qu.: 95.21 3rd Qu.: 32.56
Max. :100.00 Max. :9.00 Max. :155.59 Max. : 98.10
> summary(data.Sim[data.Sim$Y==0,])
id Time csse vase Y
Min. : 0.50 Min. :0.00 Min. :-503.75 Min. : 27.17 0:400
1st Qu.:12.88 1st Qu.:2.25 1st Qu.: -18.55 1st Qu.: 51.49 1: 0
Median :25.25 Median :4.50 Median : 10.88 Median : 67.10
Mean :25.25 Mean :4.50 Mean : 49.01 Mean : 67.13
3rd Qu.:37.62 3rd Qu.:6.75 3rd Qu.: 122.14 3rd Qu.: 83.80
Max. :50.00 Max. :9.00 Max. : 564.92 Max. :125.68
Just focus on the last summaries, which my question is how to reduce the variation on csse on group 0 Y==0
that not includes the range of group 1 .
May any one has answer, please?
Thanks