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
|
# ==================================================================================
# = 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)
# 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) {
exoFree <- NULL
if (slopes) {
exoFree <- matrix(TRUE, length(manifests), 2)
exoFree[5,] <- FALSE
exoFree[4,1] <- FALSE
}
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, exoFree=exoFree),
# 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)
omxCheckCloseEnough(jointRAM1$data$observedStats$z1.vcov,
vcov(lm(z1 ~ z1c + z2c, jointdata)))
# glm binomial is set up a little different
# omxCheckCloseEnough(jointRAM1$data$observedStats$z2.vcov,
# vcov(glm(z2 ~ z1c + z2c, jointdata, family = binomial(link="probit"))),
# .01)
omxCheckCloseEnough(jointRAM1$data$observedStats$z3.vcov,
vcov(lm(z3 ~ z1c + z2c, jointdata)))
perm <- c(2,3,4,1)
omxCheckCloseEnough(jointRAM1$data$observedStats$z4.vcov,
vcov(MASS::polr(z4 ~ z2c, jointdata, method="probit"))[perm,perm], 1e-7)
omxCheckCloseEnough(jointRAM1$data$observedStats$z5.vcov,
vcov(MASS::polr(z5 ~ 1, jointdata, method="probit")), 1e-8)
summary(jointRAM1)
# plot(jointRAM1) # (if using umx)
# Where do these come from?
#cat(deparse(round(unname(coef(jointRAM1)),4)))
param <- c(0.5929, 0.5773, 0.661, 0.6089, 0.1722, 0.0851, 0.0274, 0.5475,
0.4579, 7.9062, 2.0767, 0.0563, -0.3714, 0.1242, 0.7867, -0.6529,
-0.2965)
omxCheckCloseEnough(coef(jointRAM1), param, 1e-3)
#cat(deparse(round(unname(c(jointRAM1$output$standardErrors)),4)))
param.se <- c(0.0616, 0.1021, 0.0682, 0.094, 0.0657, 0.0566, 0.0682, 0.057,
0.065, 0.0754, 0.0601, 0.0741, 0.0783, 0.0728, 0.0925, 0.0653,
0.0586)
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, .05)
# ===============================================================
# = 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)
}
}
|
