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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
# ==================================================================================
# = This is a 1-factor model of 2 continuous variables, and 3 factors =
# = Two definition variables are created (coded with a "c" suffix ?) =
# = In RAM models, you add a definition variable as a fake latent with no variance =
# = with its mean set to "data.defnVar" =
# ==================================================================================
# TODO: Change the file name: "exoPredWLS.R" is not discoverable
library(OpenMx)
library(testthat)
# Make the results repeatable (unclear if this is needed?)
set.seed(1)
data("jointdata", package ="OpenMx", verbose= TRUE)
# ==========================================
# = Here's what our input data looks like: =
# ==========================================
str(jointdata)
# Make z1c = z1 + some noise;
jointdata$z1c <- with(jointdata, z1 * .1 + rnorm(nrow(jointdata)))
# and z2c which is just a column of noise with a mean = the mean factor level of z2 (no clue why as yet...)
jointdata$z2c <- with(jointdata, rnorm(nrow(jointdata), mean=unclass(z2)*.2))
# =============================================
# = Build a function to create WLS RAM models =
# =============================================
# note: This makes the script more complex, but will allow some programatic generation below
buildModel <- function(manifests= paste0('z', 1:5), type = 'WLS', slopes= TRUE) {
jr <- mxModel("JointRAM", type= "RAM",
manifestVars = manifests,
# TODO: is this legit (to use a var name found in the data as a latent name?)
latentVars = c('G','z1c','z2c'),
mxPath('one', c('z1', 'z3'), free= TRUE, labels= c('z1','z3')),
mxPath(paste0('z', c(1,3)), arrows= 2, free= TRUE, values= .5),
# ordinal variables have fixed variance and no mean
# TODO: why var = .5?
mxPath(paste0('z', c(2,4,5)), arrows= 2, free= FALSE, values= .5),
# latent scaled to var of 1 (mean is fixed at zero by default)
mxPath('G', arrows= 2, values= 1, free= FALSE),
mxPath('G', to = manifests, free= TRUE, values= 1, lbound= 0),
mxThreshold('z2', 1, free=TRUE, labels="z2_thresh1"),
mxThreshold('z4', 3, free=TRUE, labels=paste0("z4_thresh",1:3)),
mxThreshold('z5', 2, free=TRUE, labels=paste0("z5_thresh",1:2)),
# Note: Data are raw, despite our using WLS
mxData(jointdata, type = "raw", verbose= 0L),
# Note: this call to mxFitFunctionWLS is all that's
# required to make a model into WLS!
mxFitFunctionWLS(type)
)
# TODO: "slopes" is means?
if(slopes){
jr <- mxModel(jr,
mxPath('one', to = 'z1c', free= FALSE, labels= "data.z1c"),
mxPath('one', to = 'z2c', free= FALSE, labels= "data.z2c"),
mxPath('z1c', to = 'z1', labels= "r1"),
mxPath('z2c', to = 'z2', labels= "r2"))
}
# TODO: Shouldn't we have a call to mxExpectationRAM or LISREL??? here???
# mxExpectationRAM(M = NA, dimnames = NA, thresholds = "T", threshnames = ???)
return(jr)
}
jointRAM1 <- buildModel()
jointRAM1 <- mxRun(jointRAM1)
summary(jointRAM1)
expect_equal(length(mxGetExpected(jointRAM1, 'standvector')),
summary(jointRAM1)$observedStatistics)
# plot(jointRAM1) # (if using umx)
# Where do these come from?
# cat(deparse(round(unname(coef(jointRAM1)), 4)))
param <- c(0.5809, 0.5898, 0.6528, 0.6061, 0.1744, 0.0871, 0.0504, 0.5475,
0.4638, 7.9037, 2.0759, 0.0624, -0.3663, 0.127, 0.7916, -0.6475,
-0.296)
omxCheckCloseEnough(coef(jointRAM1), param, 1e-3)
param.se <- c(0.0613, 0.1056, 0.0684, 0.0941, 0.0665, 0.0541, 0.0648,
0.0559, 0.0645, 0.0736, 0.0593, 0.0745, 0.0776, 0.0726, 0.0921,
0.0655, 0.0585)
omxCheckCloseEnough(c(jointRAM1$output$standardErrors), param.se, 1e-3)
# ===============================================================
# = Switch jointRAM1 from MxFitFunctionWLS to an ML fitfunction =
# = to compare the estimates of these two approaches =
# ===============================================================
jointRAM2 <- mxModel(jointRAM1, mxFitFunctionML())
jointRAM2 <- mxRun(jointRAM2)
summary(jointRAM2)
estCmp <- cbind(coef(jointRAM2), coef(jointRAM1))
omxCheckCloseEnough(cor(estCmp)[2,1], 1, 1e-4)
seCmp <- cbind(jointRAM2$output$standardErrors, jointRAM1$output$standardErrors)
omxCheckCloseEnough(cor(seCmp)[2,1], 1, .04)
# ===============================================================
# = Iterate over model types allowed by the buildModel function =
# ===============================================================
mani = paste0('z', 5:1)
for (slopes in c(TRUE,FALSE)) {
for (type in c('WLS','DWLS','ULS')) {
jm3 <- buildModel(type=type, slopes=slopes)
#jm3$data$verbose <- 1L
jm3 <- mxRun(jm3)
jm4 <- mxModel(buildModel(manifests = mani, type=type, slopes=slopes), jm3$data)
jm4 <- mxRun(jm4)
}
}
# =========================
# = Make a LISREL version =
# =========================
jointLISREL <- mxModel("JointLISREL", type="LISREL",
manifestVars = list(endogenous= paste0('z',1:5)),
latentVars = list(endogenous= c('G','z1c','z2c')),
mxPath('one', paste0('z', c(1,3)), free= TRUE, labels= c('z1','z3')),
mxPath(paste0('z', c(1,3)), arrows= 2, free=TRUE, values= .5),
mxPath(paste0('z', c(2,4,5)), arrows= 2, free=FALSE, values= .5),
mxPath('G', arrows= 2, values= 1, free= FALSE),
mxPath('G', paste0('z', 1:5), free= TRUE, values= 1, lbound= 0, ubound= 4),
mxPath('one', 'z1c', free= FALSE, labels= "data.z1c"),
mxPath('one', 'z2c', free= FALSE, labels= "data.z2c"),
mxPath('z1c', 'z1', labels= "r1"),
mxPath('z2c', 'z2', labels= "r2"),
mxThreshold('z2', 1, free=TRUE, labels="z2_thresh1"),
mxThreshold('z4', 3, free=TRUE, labels=paste0("z4_thresh",1:3)),
mxThreshold('z5', 2, free=TRUE, labels=paste0("z5_thresh",1:2)),
mxData(jointdata, "raw", verbose=0L),
mxFitFunctionWLS()
# Shouldn't we have a call to mxExpectationRAM or LISREL??? here???
# mxExpectationRAM(M = NA, dimnames = NA, thresholds = "T", threshnames = ???)
)
# =================================================================
# = TODO: How would a user ever discover this?: What does it do?? =
# =================================================================
# jointLISREL$expectation$verbose <- 1L
jointLISREL <- mxRun(jointLISREL)
expect_equal(length(mxGetExpected(jointLISREL, 'standvector')),
summary(jointLISREL)$observedStatistics)
# Compare parameter estimates from the RAM and LISREL models
print(max(abs(coef(jointLISREL) - coef(jointRAM1))))
omxCheckCloseEnough(coef(jointLISREL), coef(jointRAM1), 2e-5)
# tmp <- c(jointRAM1$output$standardErrors)
# names(tmp) <- c()
# cat(deparse(round(tmp,4)))
numPeople <- 100
personData <- data.frame(
snp=rnorm(numPeople),
isMale=as.numeric(rbinom(numPeople,1,.5)),
phenotype=rnorm(numPeople),
snpsex=0)
m1 <- mxModel(
"gwsem", type="RAM",
manifestVars = c('phenotype'),
latentVars = c('snp', 'sex', 'snpsex'),
# residual variances
mxPath(c('phenotype'), arrows=2, values=1),
mxPath('one', 'sex', free=FALSE, labels="data.isMale"),
mxPath('one', 'snp', free=FALSE, labels="data.snp"),
mxAlgebra(data.snp * data.isMale, name="snpsexAlg",
dimnames=list(NULL, 'snpsex')),
mxPath('one', 'snpsex', free=FALSE, labels="data.snpsex"),
mxPath('one', 'phenotype'),
mxPath(c('snp','sex', 'snpsex'), 'phenotype'),
mxData(personData, type="raw", algebra='snpsexAlg'),
mxFitFunctionWLS(allContinuousMethod = "marginals"))
m1 <- mxRun(m1)
personData$snpsex <- with(personData, snp * isMale)
c1 <- coef(lm(phenotype ~ isMale + snp + snpsex, personData))
omxCheckCloseEnough(c1[-1], coef(m1)[c('gwsem.A[1,3]', 'gwsem.A[1,2]', 'gwsem.A[1,4]')], 1e-6)
|
