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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
# Maes, H. H., Neale, M. C., Kendler, K. S., Hewitt, J. K., Silberg, J. L., Foley, D. L., ... & Eaves, L. J. (1998). Assortative mating for major psychiatric diagnoses in two population-based samples. Psychological medicine, 28(6), 1389-1401.
#
# https://vipbg.vcu.edu/vipbg/Articles/psychological-assortative-1998.pdf
library(OpenMx)
library(testthat)
diagno <- toupper(c('alc', 'gad', 'mdd', 'pan', 'pho'))
parent <- c('wife', 'husband')
c1 <- apply(expand.grid(diagno, parent)[,c(2,1)], 1, paste, collapse='')
wife <- c1[1:5]
husband <- c1[6:10]
c2 <- matrix(0, length(c1), length(c1),
dimnames=list(c1,c1))
# Table 2 ABD upper triangle
c2['wifeALC', 'wifeGAD'] <- .399
c2['wifeALC', 'wifeMDD'] <- .341
c2['wifeALC', 'wifePAN'] <- .276
c2['wifeALC', 'wifePHO'] <- .284
c2['wifeGAD', 'wifeMDD'] <- .582
c2['wifeGAD', 'wifePAN'] <- .372
c2['wifeGAD', 'wifePHO'] <- .394
c2['wifeMDD', 'wifePAN'] <- .31
c2['wifeMDD', 'wifePHO'] <- .304
c2['wifePAN', 'wifePHO'] <- .371
# Table 2 ABD lower triangle
c2['husbandGAD', 'husbandALC'] <- .295
c2['husbandMDD', 'husbandALC'] <- .341
c2['husbandMDD', 'husbandGAD'] <- .576
c2['husbandPAN', 'husbandALC'] <- .283
c2['husbandPAN', 'husbandGAD'] <- .314
c2['husbandPAN', 'husbandMDD'] <- .576
c2['husbandPHO', 'husbandALC'] <- .249
c2['husbandPHO', 'husbandGAD'] <- .220
c2['husbandPHO', 'husbandMDD'] <- .181
c2['husbandPHO', 'husbandPAN'] <- .276
# Table 2 ABD
c2[wife,wife] + c2[husband, husband]
# Table 3 ABD
c2['wifeALC', husband] <- c(.119, .064, .19, .029, .005)
c2['wifeGAD', husband] <- c(.209, .208, .207, .119, .133)
c2['wifeMDD', husband] <- c(.132, .088, .162, .112, .105)
c2['wifePAN', husband] <- c(.064, .092, .151, .219, .102)
c2['wifePHO', husband] <- c(-.079, .13, .221, -.016, .064)
# Repair symmetry
c2[husband,wife] <- t(c2[wife, husband])
c2[wife, wife][lower.tri(diag(length(diagno)))] <-
t(c2[wife, wife])[lower.tri(diag(length(diagno)))]
c2[husband, husband][upper.tri(diag(length(diagno)))] <-
t(c2[husband, husband])[upper.tri(diag(length(diagno)))]
diag(c2) <- 1
c3 <- matrix(0, length(c1), length(c1),
dimnames=list(c1,c1))
# Table 2 AFT upper triangle
c3['wifeALC', 'wifeGAD'] <- .378
c3['wifeALC', 'wifeMDD'] <- .313
c3['wifeALC', 'wifePAN'] <- .076
c3['wifeALC', 'wifePHO'] <- .095
c3['wifeGAD', 'wifeMDD'] <- .709
c3['wifeGAD', 'wifePAN'] <- .471
c3['wifeGAD', 'wifePHO'] <- .228
c3['wifeMDD', 'wifePAN'] <- .367
c3['wifeMDD', 'wifePHO'] <- .227
c3['wifePAN', 'wifePHO'] <- .371
# Table 2 AFT lower triangle
c3['husbandGAD', 'husbandALC'] <- .125
c3['husbandMDD', 'husbandALC'] <- .185
c3['husbandMDD', 'husbandGAD'] <- .64
c3['husbandPAN', 'husbandALC'] <- .267
c3['husbandPAN', 'husbandGAD'] <- .483
c3['husbandPAN', 'husbandMDD'] <- .467
c3['husbandPHO', 'husbandALC'] <- .194
c3['husbandPHO', 'husbandGAD'] <- .312
c3['husbandPHO', 'husbandMDD'] <- .149
c3['husbandPHO', 'husbandPAN'] <- .344
# Table 2 AFT
c3[wife,wife] + c3[husband, husband]
# Table 3 AFT
c3['wifeALC', husband] <- c(.068, .213, .267, -.047, -.078)
c3['wifeGAD', husband] <- c(.107, -.014, .169, .116, -.006)
c3['wifeMDD', husband] <- c(.251, .014, .121, -.196, .074)
c3['wifePAN', husband] <- c(.002, .114, .015, .138, .062)
c3['wifePHO', husband] <- c(.003, .14, .033, .028, .071)
# Repair symmetry
c3[husband,wife] <- t(c3[wife, husband])
c3[wife, wife][lower.tri(diag(length(diagno)))] <-
t(c3[wife, wife])[lower.tri(diag(length(diagno)))]
c3[husband, husband][upper.tri(diag(length(diagno)))] <-
t(c3[husband, husband])[upper.tri(diag(length(diagno)))]
diag(c3) <- 1
diag2d <- matrix(apply(expand.grid(diagno, diagno)[,c(2,1)],
1, paste, collapse=""),
length(diagno), length(diagno))
paths <- c(
mxPath(wife, arrows=2, connect="unique.bivariate",
lbound=-.99, ubound=.99,
labels=paste0('w', diag2d[lower.tri(diag2d)])),
mxPath(husband, arrows=2, connect="unique.bivariate",
lbound=-.99, ubound=.99,
labels=paste0('h', diag2d[lower.tri(diag2d)])),
mxPath(wife, husband, arrows=0, connect="unique.pairs",
lbound=-.99, ubound=.99,
labels=paste0('d', c(diag2d))))
abd <- mxModel(
'abd', type="RAM",
manifestVars = c1,
mxData(c2, type = "cov", numObs = 854),
mxPath(c1, arrows=2, values=1, free=FALSE), paths)
aft <- mxModel(
'aft', type="RAM",
manifestVars = c1,
mxData(c3, type = "cov", numObs = 568),
mxPath(c1, arrows=2, values=1, free=FALSE), paths)
fig6 <- mxModel("fig6", abd, aft,
mxFitFunctionMultigroup(c("abd", "aft")))
if (0) {
fig6 <- mxOption(fig6,"Always Checkpoint","Yes")
fig6 <- mxOption(fig6,"Checkpoint Units","evaluations")
fig6 <- mxOption(fig6,"Checkpoint Count",1)
fig6 <- mxOption(fig6, "Checkpoint Fullpath", "/dev/fd/2")
}
fig6 <- mxRun(fig6)
fig6r <- mxRefModels(fig6, run=TRUE)
fig6s <- summary(fig6, refModels=fig6r)
expect_equal(fig6s$ChiDoF, 65)
expect_equal(fig6s$Chi, 544.0544, 1e-3)
expect_equal(fig6s$CFI, .8756, 1e-4)
expect_equal(length(mxGetExpected(fig6, "covariance")), 2)
|
