-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsimulation code in R
76 lines (69 loc) · 4.42 KB
/
simulation code in 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
# targeted h2 = 0.5, c2 = 0.35, e2 = 0.15
library(MASS); library(lme4); library(nlme)
nr <- 100 # number of simulations per case
rDZ.t <- rep(0, nr); rMZ.t <- rep(0, nr) # rMZ & rDZ for trial-level modeling
rDZ.c <- rep(0, nr); rMZ.c <- rep(0, nr) # rMZ & rDZ for condition-level modeling through aggregation
rDZ.1 <- rep(0, nr); rMZ.1 <- rep(0, nr) # rMZ & rDZ for condition-level modeling w adjustment 1
rDZ.2 <- rep(0, nr); rMZ.2 <- rep(0, nr) # rMZ & rDZ for condition-level modeling w adjustment 2
rho.DZ <- 0.5; rho.MZ <- 0.375
gamma1 <- 0.8; gamma2 <- 0.5
tauMZ2 <- 0.7; tauDZ2 <- 0.5 # cross-family variance for MZ & DZ under each condition
# covariance matrix
Sigma <- matrix(c(1, gamma1, gamma2, rho.DZ, 0,0,0,0, # DZ; condition 1, twin 1 rDZ: (gamma1-rho.DZ)/(1-gamma2)
gamma1, 1, rho.DZ, gamma2, 0,0,0,0, # DZ; condition 1, twin 2 targeted rDZ: 0.6
gamma2, rho.DZ, 1, gamma1, 0,0,0,0, # DZ; condition 2, twin 1
rho.DZ, gamma2, gamma1, 1, 0,0,0,0, # DZ; condition 2, twin 2
0,0,0,0, 1, gamma1, gamma2, rho.MZ, # MZ; condition 1, twin 1 rMZ: (gamma1-rho.MZ)/(1-gamm2)
0,0,0,0, gamma1, 1, rho.MZ, gamma2, # MZ; condition 1, twin 2 targeted rMZ: 0.85
0,0,0,0, gamma2, rho.MZ, 1, 0.8, # MZ; condition 2, twin 1
0,0,0,0, rho.MZ, gamma2, 0.8, 1), 8, 8) # MZ; condition 2, twin 2
Sigma <- rbind(sweep(Sigma[1:4,], MARGIN=2, rep(tauDZ2,4), `*`),
sweep(Sigma[5:8,], MARGIN=2, rep(tauMZ2,4), `*`))
for(Rv in c(1,3,10)) # within-family cross-trial Rv ratio relative to cross-family Rv
for(ns in c(50,250,500)) { # number of family: 2*ns
s1 <- 2*ns; s2 <- 4*ns # s1: number of families; s2 <- number of twins
for(nt in c(50,100,500)) { # number of trials
dat <- data.frame(
fam = rep(rep(paste0('f',1:s1), each=4), nt),
mem = rep(c('m1','m2'), s2*nt),
cond = rep(c(-0.5,-0.5,0.5,0.5),s1*nt),
g1 = rep(c(rep(1,4), rep(0,4)),ns*nt),
g2 = rep(c(rep(0,4), rep(1,4)),ns*nt),
gg = rep(c(rep('g1',4), rep('g2',4)),ns*nt),
y = 0)
for(ii in 1:nr) {
dd <- mvrnorm(n=ns, mu=rep(0,8), Sigma=Sigma)
dat$y <- c(t(dd)) + rnorm(2*s2*nt, 0, Rv)
# trial-level modeling
try(m1 <- lmer(y ~ 0+g1+g2+g1:cond+g2:cond+
(0+g1|fam)+(0+g2|fam)+(0+g1:cond|fam)+(0+g2:cond|fam)+
(0+g1|fam:mem)+(0+g2|fam:mem)+(0+g1:cond|fam:mem)+(0+g2:cond|fam:mem), data=dat))
if(!is.null(m1)) {
rDZ.t[ii] <- VarCorr(m1)$fam.1[1]/(VarCorr(m1)$fam.1[1]+VarCorr(m1)$`fam.mem.1`[1])
rMZ.t[ii] <- VarCorr(m1)$fam[1]/(VarCorr(m1)$fam[1]+VarCorr(m1)$`fam.mem`[1])
}
# condition-level modeling via aggregation
d2 <- aggregate(y~fam+mem+gg+g1+g2, with(dat, data.frame(fam=fam, g1=g1, g2=g2,
mem=mem, gg=gg, y=ifelse(cond==0.5, 1, -1)*y)), function(x) 2*mean(x))
# assume heterogeneity of inter-individual variability
try(fm <- lme(y~0+g1+g2,random=list(fam=~0+g1,fam=~0+g2), data=d2, weights=varIdent(form=~1|gg)))
if(!is.null(fm)) {
v <- VarCorr(fm)[c('g1','g2','Residual'),"Variance"]
w <- coef(fm$modelStruct$varStruct, uncons=F, allCoef=T)
t <- as.numeric(v[c('g1','g2')])
res <- t/(t+as.numeric(v['Residual'])*w[c('g1','g2')]^2)
rDZ.c[ii] <- res[1]
rMZ.c[ii] <- res[2]
rDZ.1[ii] <- t[1]/(t[1] + max(as.numeric(v['Residual'])*w['g1']^2 - 2*Rv^2/nt,0))
rMZ.1[ii] <- t[2]/(t[2] + max(as.numeric(v['Residual'])*w['g2']^2 - 2*Rv^2/nt,0))
ee <- aggregate(y~fam+mem, data=dat, var)
rDZ.2[ii] <- t[1]/(t[1] + max(as.numeric(v['Residual'])*w['g1']^2 - 2*(nt-1)/(s2*nt-1)*sum(ee$y)/nt,0)) # s2: #individuals/twins
rMZ.2[ii] <- t[2]/(t[2] + max(as.numeric(v['Residual'])*w['g2']^2 - 2*(nt-1)/(s2*nt-1)*sum(ee$y)/nt,0))
}
}
tmp <- data.frame(cbind(Rv, ns, nt, rDZ.t, rMZ.t, rDZ.c, rMZ.c, rDZ.1, rMZ.1, rDZ.2, rMZ.2))
names(tmp) <- c('Rv', 'ns', 'nt', 'rDZ.t', 'rMZ.t', 'rDZ.c', 'rMZ.c', 'rDZ.1', 'rMZ.1', 'rDZ.2', 'rMZ.2')
write.table(tmp, file = 'simulation.results.txt', quote = FALSE, row.names = FALSE,
col.names = !file.exists('simulation.results.txt'), append = TRUE)
}
}