/
simdatx.R
118 lines (90 loc) · 3.25 KB
/
simdatx.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
######## simulation for x grouping ######
### time is simulated from a uniform distribution Unif(0,1)
### yih = xietai + mui(tih) + eih
## nx is the number of x
## etag group parameter for x ## nx \times ng matrix, each column is for one group
##### parametric part is the same for the same individual ####
simdatx = function(xlist = list(nx = 2, meanx = 0, sdx = 1, etag),
sig2, lamj, mvec = c(10,20), ncl = 50,
funlist, eigenlist, seed = 2228)
{
set.seed(seed)
ngroup = length(funlist)
nobs = ncl*ngroup
if(mvec[1]== mvec[2])
{
nsub = rep(mvec[1], nobs)
}else{
nsub = sample(mvec[1]:mvec[2], nobs, replace = TRUE)
}
ntotal = sum(nsub)
group = rep(1:ngroup, each = ncl)
xmat = matrix(rnorm(nobs*xlist$nx,xlist$meanx, xlist$sdx), nrow = nobs)
dat = data.frame(group = rep(group, nsub),
ind = rep(1:(ngroup*ncl),nsub),
time = runif(ntotal),
x = xmat[rep(1:nobs, nsub),],
obs = rep(0, ntotal),
mean = rep(0, ntotal),
meanx = rep(0, ntotal))
### mean is the mean curve function #####
colnames(dat)[4:(xlist$nx+3)] = paste0("x",1:xlist$nx)
rande = sqrt(sig2) * rnorm(ntotal)
for(j in 1:ngroup)
{
dat$mean[dat$group==j] = funlist[[j]](dat$time[dat$group==j])
dat$meanx[dat$group==j] = as.matrix(dat[dat$group==j,4:(xlist$nx+3)]) %*% xlist$etag[,j]
}
vi = rep(0, ntotal)
for(i in 1:(ncl*ngroup))
{
timei = dat$time[dat$ind==i]
vi[dat$ind ==i] = sqrt(lamj[1]) * rnorm(1) * eigenlist[[1]](timei) +
sqrt(lamj[2]) * rnorm(1) * eigenlist[[2]](timei)
}
dat$obs = dat$meanx + dat$mean + vi + rande
return(dat)
}
####### parametric part is not the same for the same individual ####
### yih = xih etai + mui(tih) + eih
simdatx2 = function(xlist = list(nx = 2, meanx = 0, sdx = 1, etag),
sig2, lamj, mvec = c(10,20), ncl = 50,
funlist, eigenlist, seed = 2228)
{
set.seed(seed)
ngroup = length(funlist)
nobs = ncl*ngroup
if(mvec[1]== mvec[2])
{
nsub = rep(mvec[1], nobs)
}else{
nsub = sample(mvec[1]:mvec[2], nobs, replace = TRUE)
}
ntotal = sum(nsub)
group = rep(1:ngroup, each = ncl)
xmat = matrix(rnorm(ntotal*xlist$nx,xlist$meanx, xlist$sdx), nrow = ntotal)
dat = data.frame(group = rep(group, nsub),
ind = rep(1:(ngroup*ncl),nsub),
time = runif(ntotal),
x = xmat,
obs = rep(0, ntotal),
mean = rep(0, ntotal),
meanx = rep(0, ntotal))
### mean is the mean curve function #####
colnames(dat)[4:(xlist$nx+3)] = paste0("x",1:xlist$nx)
rande = sqrt(sig2) * rnorm(ntotal)
for(j in 1:ngroup)
{
dat$mean[dat$group==j] = funlist[[j]](dat$time[dat$group==j])
dat$meanx[dat$group==j] = as.matrix(dat[dat$group==j,4:(xlist$nx+3)]) %*% xlist$etag[,j]
}
vi = rep(0, ntotal)
for(i in 1:(ncl*ngroup))
{
timei = dat$time[dat$ind==i]
vi[dat$ind ==i] = sqrt(lamj[1]) * rnorm(1) * eigenlist[[1]](timei) +
sqrt(lamj[2]) * rnorm(1) * eigenlist[[2]](timei)
}
dat$obs = dat$meanx + dat$mean + vi + rande
return(dat)
}