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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
|
#
# Copyright 2007-2018 by the individuals mentioned in the source code history
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#------------------------------------------------------------------------------
# Author: Michael D. Hunter
# Date: 2011.03.24
# Filename: LISRELFactorRegressionWithMeans_MatrixRaw.R
# Purpose: Create a test for the mxExpectationLISREL function using FIML. This
# test includes a model for the means and uses raw data. The known
# parameters are taken from model 3 (AKA threeLatentMultipleReg1 and
# threeLatentMultipleReg1Out) in
# models/passing/IntroSEM-ThreeLatentMultipleRegTest1.R. This is the same test as
# models/passing/LISRELFactorRegressionWithMeans_Matrix.R only it uses the raw
# data.
#------------------------------------------------------------------------------
# Revision History:
# Sat Mar 24 01:27:41 Central Daylight Time 2012 -- Michael Hunter created complete file from
# models/passing/LISRELFactorRegressionWithMeans_Matrix.R
# Sat Mar 24 01:51:53 Central Daylight Time 2012 -- Michael Hunter corrected the expected parameters
# Thu Jul 19 14:21:05 Central Daylight Time 2012 -- Michael Hunter corrected typo. Note: from March 24 until Jul 19, I thought I had a broken FIML re-ordering of data. In fact, I had a broken test of the re-ordering and FIML re-ordering worked all along.
#
#--------------------------------------------------------------------
# Load OpenMx
require(OpenMx)
#--------------------------------------------------------------------
# Read in and set up the data
IndManExo <- 1:8
IndManEnd <- 9:12
# The data
data(latentMultipleRegExample1)
# Rearange Columns to separate exogenous and endogenous variables
#rawlisdat <- latentMultipleRegExample1[, c(IndManEnd, IndManExo)]
rawlisdat <- latentMultipleRegExample1
# Take covariance and means
covlisdat <- cov(rawlisdat)
mealisdat <- colMeans(rawlisdat)
# Number of manifest and latent exogenous and endogenous variables
numLatExo <- 2
numLatEnd <- 1
numManExo <- 8
numManEnd <- 4
# Dimnames
LatExo <- paste('xi', 1:numLatExo, sep='')
LatEnd <- paste('eta', 1:numLatEnd, sep='')
ManExo <- names(rawlisdat)[IndManExo]
ManEnd <- names(rawlisdat)[IndManEnd]
#--------------------------------------------------------------------
# Specify the 13 extended LISREL matrices
lx <- mxMatrix("Full", numManExo, numLatExo,
free=c(F,T,T,T,F,F,F,F,F,F,F,F,F,T,T,T),
values=c(1, .2, .2, .2, 0, 0, 0, 0, 0, 0, 0, 0, 1, .2, .2, .2),
labels=c( paste('l', 1, 1:4, sep=''), rep(NA, 8), paste('l', 2, 5:8, sep='')),
name='LX',
dimnames=list(ManExo, LatExo)
) #DONE
ly <- mxMatrix("Full", numManEnd, numLatEnd,
free=c(F,T,T,T),
values=c(1, .2, .2, .2),
labels= paste('l', 3, 9:12, sep=''),
name='LY',
dimnames=list(ManEnd, LatEnd)
) #DONE
be <- mxMatrix("Zero", numLatEnd, numLatEnd, name='BE', dimnames=list(LatEnd, LatEnd)) #DONE
ga <- mxMatrix("Full", numLatEnd, numLatExo,
free=T,
values=.2,
labels=c('b13', 'b23'),
name='GA',
dimnames=list(LatEnd, LatExo)
) #DONE
ph <- mxMatrix("Symm", numLatExo, numLatExo,
free=c(T,T,T),
values=c(.8, .3, .8),
labels=c('varF1', 'covF1F2', 'varF2'),
name='PH',
dimnames=list(LatExo, LatExo)
) #DONE
ps <- mxMatrix("Symm", numLatEnd, numLatEnd,
free=T,
values=.8,
labels='varF3',
name='PS',
dimnames=list(LatEnd, LatEnd)
) #DONE
td <- mxMatrix("Diag", numManExo, numManExo,
free=T,
values=.8,
labels=paste('d', 1:8, sep=''),
name='TD',
dimnames=list(ManExo, ManExo)
) #DONE
te <- mxMatrix("Diag", numManEnd, numManEnd,
free=T,
values=.8,
labels=paste('e', 9:12, sep=''),
name='TE',
dimnames=list(ManEnd, ManEnd)
) #DONE
th <- mxMatrix("Zero", numManExo, numManEnd, name='TH', dimnames=list(ManExo, ManEnd)) #DONE
tx <- mxMatrix("Full", numManExo, 1,
free=T,
values=.1,
labels=paste('m', 1:8, sep=''),
name='TX',
dimnames=list(ManExo, "TXMeans")
) #DONE
ty <- mxMatrix("Full", numManEnd, 1,
free=T,
values=.1,
labels=paste('m', 9:12, sep=''),
name='TY',
dimnames=list(ManEnd, "TYMeans")
) #DONE
ka <- mxMatrix("Zero", numLatExo, 1, name='KA', dimnames=list(LatExo, "KAMeans")) #DONE
al <- mxMatrix("Zero", numLatEnd, 1, name='AL', dimnames=list(LatEnd, "ALMeans")) #DONE
#--------------------------------------------------------------------
# Define the model
lmod <- mxModel(
name='LISREL Factor Regression Model with Means',
mxData(observed=rawlisdat, type='raw', means=mealisdat),
lx, ly, be, ga, ph, ps, td, te, th, tx, ty, ka, al,
mxExpectationLISREL(LX=lx$name, LY=ly$name, BE=be$name,
GA=ga$name, PH=ph$name, PS=ps$name, TD=td$name,
TE=te$name, TH=th$name, TX=tx$name, TY=ty$name,
KA=ka$name, AL=al$name),
mxFitFunctionML()
)
#--------------------------------------------------------------------
# Run the model
# Uncomment the following lines when debugbing
#lmodRun <- mxRun(lmod, onlyFrontend=TRUE) # This runs fine.
#lmod <- mxOption(lmod, "Calculate Hessian", "No")
#lmod <- mxOption(lmod, "Standard Errors", "No")
#lmod <- mxOption(lmod, "Major iterations", 0)
lmodRun <- mxRun(lmod)
summary(lmodRun)
#--------------------------------------------------------------------
# Compare the estimate parameters to known values
# Recal that The known parameters are taken from model 3 (AKA
# threeLatentMultipleReg1 and threeLatentMultipleReg1Out) in
# models/passing/IntroSEM-ThreeLatentMultipleRegTest1.R.
expectedParam <- c(
0.809021, 1.148334, 1.305356, 0.805289, 1.232038, 1.189698,
0.87486, 0.786604, 0.701819, 0.990778, 0.456683, 1.92904,
0.153399, 1.316878, 0.770605, 1.133502, 1.067425, 1.063866,
1.053663, 0.847609, 0.761789, 1.186354, 1.009936, 0.942481,
0.971017, 0.880311, 0.956065, 0.060813, 0.037374, -0.048679,
-0.013194, 0.193348, 0.220002, 0.256785, 0.171892, 0.166191,
0.23483, 0.173025, 0.157615
)
names(expectedParam) <- c("l12", "l13", "l14", "l26", "l27", "l28", "l310",
"l311", "l312", "b13", "b23", "varF1", "covF1F2", "varF2",
"varF3", "d1", "d2", "d3", "d4", "d5", "d6",
"d7", "d8", "e9", "e10", "e11", "e12", "m1",
"m2", "m3", "m4", "m5", "m6", "m7", "m8",
"m9", "m10", "m11", "m12")
match(names(expectedParam), names(lmodRun$output$estimate))
omxCheckCloseEnough(lmodRun$output$estimate, expectedParam, epsilon=0.001)
#--------------------------------------------------------------------
# End
|
