-
Notifications
You must be signed in to change notification settings - Fork 0
/
ACE_M_separate_moderator_values.R
132 lines (106 loc) · 4.71 KB
/
ACE_M_separate_moderator_values.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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#y_dz = Item responses of DZ twins (matrix)
#y_mz = Item responses of MZ twins (matrix)
#n_mz = Number of MZ twin pairs
#n_dz = Number of DZ twin pairs
#n_items = Number of phenotypic items administered
#b = Vector with item difficulty parameters, assumed known here
#x_MZ = moderator variable values for every MZ twin (matrix)
#x_DZ = moderator variable values for every DZ twin (matrix)
#Required structure of the y_mz/y_dz data matrix:
#y_dz[i,k] = kth datapoint from the ith DZ twin pair
#y_mz[i,k] = kth datapoint from the ith MZ twin pair
#This results in a matrix of n_mz (or, in case of y_dz, n_dz)
#rows and 2*n_items columns. e.g. y_mz[1,22] is the response
#of MZ twin 1 from family 1 to item 22 if n_items = 22
#Required structure of the x_MZ/x_DZ data matrix:
#x_MZ[i,1:2] = Vector with moderator variable values
#for the first (x_MZ[i,1]) and second (x_MZ[i,2])
#twin of every ith MZ twin pair.
#x_DZ[i,1:2] = Vector with moderator variable values
#for the first (x_MZ[i,1]) and second (x_MZ[i,2])
#twin of every ith MZ twin pair.
#JAGS uses precision parameters for the variance parameters.
#Therefore, after running the script, these precision parameters
#have to be inverted. For example:
#var_a <- 1/outputAnalysis$tau_a[,,1] with the rjags package
model{
##MZ twins
for (fam in 1:n_mz){
c_mz[fam] ~ dnorm(0, 1)
a_mz[fam] ~ dnorm(0, 1)
tau_c_mz[fam,1] <- exp(beta_0c + beta_1c*x_MZ[fam,1])
tau_c_mz[fam,2] <- exp(beta_0c + beta_1c*x_MZ[fam,2])
tau_a_mz[fam,1] <- exp(beta_0a + beta_1a*x_MZ[fam,1])
tau_a_mz[fam,2] <- exp(beta_0a + beta_1a*x_MZ[fam,2])
tau_e_mz[fam,1] <- 1/(exp(beta_0e + beta_1e*x_MZ[fam,1] ))
tau_e_mz[fam,2] <- 1/(exp(beta_0e + beta_1e*x_MZ[fam,2] ))
c_mz_twin1[fam] <- c_mz[fam] * sqrt(tau_c_mz[fam,1])
c_mz_twin2[fam] <- c_mz[fam] * sqrt(tau_c_mz[fam,2])
a_mz_twin1[fam] <- a_mz[fam] * sqrt(tau_a_mz[fam,1])
a_mz_twin2[fam] <- a_mz[fam] * sqrt(tau_a_mz[fam,2])
pheno_mz[fam,1] ~ dnorm(mu_mz + beta_1m_mz * x_MZ[fam,1] +
beta_1mc_mz * x_MZ[fam,2] +
c_mz_twin1[fam] + a_mz_twin1[fam],
tau_e_mz[fam,1])
pheno_mz[fam,2] ~ dnorm(mu_mz + beta_1m_mz * x_MZ[fam,2] +
beta_1mc_mz * x_MZ[fam,1] +
c_mz_twin2[fam] + a_mz_twin2[fam],
tau_e_mz[fam,2])
#### Measurement model for phenotypic item data:
for (k in 1:n_items){
logit(p[fam,k]) <- pheno_mz[fam,1] - b[k]
y_mz[fam,k] ~ dbern(p[fam,k])
}
for (k in (n_items+1):(2*n_items)){
logit(p[fam,k]) <- pheno_mz[fam,2] - b[k-n_items]
y_mz[fam,k] ~ dbern(p[fam,k])
}
} #end MZ twins
##DZ twins
for (fam in 1:n_dz){
c_dz[fam] ~ dnorm(0,1)
a1_dz[fam] ~ dnorm(0,2)
a2_dz[fam,1] ~ dnorm(a1_dz[fam], 2)
a2_dz[fam,2] ~ dnorm(a1_dz[fam], 2)
tau_c_dz[fam,1] <- exp(beta_0c + beta_1c*x_DZ[fam,1])
tau_c_dz[fam,2] <- exp(beta_0c + beta_1c*x_DZ[fam,2])
tau_a_dz[fam,1] <- exp(beta_0a + beta_1a*x_DZ[fam,1])
tau_a_dz[fam,2] <- exp(beta_0a + beta_1a*x_DZ[fam,2])
tau_e_dz[fam,1] <- 1/(exp(beta_0e + beta_1e*x_DZ[fam,1] ))
tau_e_dz[fam,2] <- 1/(exp(beta_0e + beta_1e*x_DZ[fam,2] ))
c_dz_twin1[fam] <- c_dz[fam] * sqrt(tau_c_dz[fam,1])
c_dz_twin2[fam] <- c_dz[fam] * sqrt(tau_c_dz[fam,2])
a_dz_twin1[fam] <- a2_dz[fam,1] * sqrt(tau_a_dz[fam,1])
a_dz_twin2[fam] <- a2_dz[fam,2] * sqrt(tau_a_dz[fam,2])
pheno_dz[fam,1] ~ dnorm(mu_dz + beta_1m_dz * x_DZ[fam,1] +
beta_1mc_dz * x_DZ[fam,2] +
c_dz_twin1[fam] + a_dz_twin1[fam],
tau_e_dz[fam,1])
pheno_dz[fam,2] ~ dnorm(mu_dz + beta_1m_dz * x_DZ[fam,2] +
beta_1mc_dz * x_DZ[fam,1] +
c_dz_twin2[fam] + a_dz_twin2[fam],
tau_e_dz[fam,2])
#### Measurement model for phenotypic item data:
for (k in 1:n_items){
logit(p_dz[fam,k]) <- pheno_dz[fam,1] - b[k]
y_dz[fam,k] ~ dbern(p_dz[fam,k])
}
for (k in (n_items+1):(2*n_items)){
logit(p_dz[fam,k]) <- pheno_dz[fam,2] - b[k-n_items]
y_dz[fam,k] ~ dbern(p_dz[fam,k])
}
} #end DZ twins
#Priors
mu_mz ~ dnorm(0, .1)
mu_dz ~ dnorm(0, .1)
beta_1m_mz ~ dnorm(0, .1)
beta_1m_dz ~ dnorm(0, .1)
beta_1mc_mz ~ dnorm(0, .1)
beta_1mc_dz ~ dnorm(0, .1)
beta_1a ~ dnorm(0, .1)
beta_1c ~ dnorm(0, .1)
beta_1e ~ dnorm(0, .1)
beta_0a ~ dnorm(-1, .5)
beta_0c ~ dnorm(-1, .5)
beta_0e ~ dnorm(-1, .5)
}