Update copyright.
[openmx:openmx.git] / models / passing / LISRELExoEndoOnly.R
1 #
2 #   Copyright 2007-2014 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16
17 #------------------------------------------------------------------------------
18 # Author: Michael D. Hunter
19 # Date: 2011.04.09
20 # Filename: LISRELExoEndoOnly.R
21 # Purpose: Create a test for the mxExpectationLISREL function using only
22 #  exogenous or only endogenous variables.  This test was created based on
23 #  models/passing/LISRELFactorRegressionWithMeans_Matrix*.R.
24 #------------------------------------------------------------------------------
25
26 # Revision History:
27 # Mon Apr 09 19:16:11 Central Daylight Time 2012 -- Michael Hunter created test from LISRELFactorRegressionWithMeans_Matrix*.R
28
29
30
31 #--------------------------------------------------------------------
32 # Load OpenMx
33
34 require(OpenMx)
35
36
37 #--------------------------------------------------------------------
38 # Read in and set up the data
39
40 IndManExo <- 1:8
41 IndManEnd <- 9:12
42
43 # The data
44 data(latentMultipleRegExample1)
45
46 # Rearange Columns to separate exogenous and endogenous variables
47 rawlisdat <- latentMultipleRegExample1[, c(IndManEnd, IndManExo)]
48 rawlisy <- latentMultipleRegExample1[, IndManEnd]
49 rawlisx <- latentMultipleRegExample1[, IndManExo]
50
51
52 # Take covariance and means
53 covlisdat <- cov(rawlisdat)
54 mealisdat <- colMeans(rawlisdat)
55
56
57 # Number of manifest and latent exogenous and endogenous variables
58 numLatExo <- 2
59 numLatEnd <- 1
60 numManExo <- 8
61 numManEnd <- 4
62
63
64 # Dimnames
65 LatExo <- paste('xi', 1:numLatExo, sep='')
66 LatEnd <- paste('eta', 1:numLatEnd, sep='')
67 ManExo <- names(rawlisdat)[(numManEnd+1):(numManEnd+numManExo)]
68 ManEnd <- names(rawlisdat)[1:numManEnd]
69
70 #--------------------------------------------------------------------
71 # Specify the 13 extended LISREL matrices
72
73
74 lx <- mxMatrix("Full", numManExo, numLatExo,
75         free=c(F,T,T,T,F,F,F,F,F,F,F,F,F,T,T,T),
76         values=c(1, .2, .2, .2, 0, 0, 0, 0, 0, 0, 0, 0, 1, .2, .2, .2),
77         labels=c( paste('l', 1, 1:4, sep=''), rep(NA, 8),  paste('l', 2, 5:8, sep='')),
78         name='LX',
79         dimnames=list(ManExo, LatExo)
80 ) #DONE
81
82 ly <- mxMatrix("Full", numManEnd, numLatEnd,
83         free=c(F,T,T,T),
84         values=c(1, .2, .2, .2),
85         labels= paste('l', 3, 9:12, sep=''),
86         name='LY',
87         dimnames=list(ManEnd, LatEnd)
88 ) #DONE
89
90 be <- mxMatrix("Zero", numLatEnd, numLatEnd, name='BE', dimnames=list(LatEnd, LatEnd)) #DONE
91
92 ga <- mxMatrix("Full", numLatEnd, numLatExo,
93         free=T,
94         values=.2,
95         labels=c('b13', 'b23'),
96         name='GA',
97         dimnames=list(LatEnd, LatExo)
98 ) #DONE
99
100 ph <- mxMatrix("Symm", numLatExo, numLatExo,
101         free=c(T,T,T),
102         values=c(.8, .3, .8),
103         labels=c('varF1', 'covF1F2', 'varF2'),
104         name='PH',
105         dimnames=list(LatExo, LatExo)
106 ) #DONE
107
108 ps <- mxMatrix("Symm", numLatEnd, numLatEnd,
109         free=T,
110         values=.8,
111         labels='varF3',
112         name='PS',
113         dimnames=list(LatEnd, LatEnd)
114 ) #DONE
115
116 td <- mxMatrix("Diag", numManExo, numManExo,
117         free=T,
118         values=.8,
119         labels=paste('d', 1:8, sep=''),
120         name='TD',
121         dimnames=list(ManExo, ManExo)
122 ) #DONE
123
124 te <- mxMatrix("Diag", numManEnd, numManEnd,
125         free=T,
126         values=.8,
127         labels=paste('e', 9:12, sep=''),
128         name='TE',
129         dimnames=list(ManEnd, ManEnd)
130 ) #DONE
131
132 th <- mxMatrix("Zero", numManExo, numManEnd, name='TH', dimnames=list(ManExo, ManEnd)) #DONE
133
134 tx <- mxMatrix("Full", numManExo, 1,
135         free=T,
136         values=.1,
137         labels=paste('m', 1:8, sep=''),
138         name='TX',
139         dimnames=list(ManExo, "TXMeans")
140 ) #DONE
141
142 ty <- mxMatrix("Full", numManEnd, 1,
143         free=T,
144         values=.1,
145         labels=paste('m', 9:12, sep=''),
146         name='TY',
147         dimnames=list(ManEnd, "TYMeans")
148 ) #DONE
149
150 ka <- mxMatrix("Zero", numLatExo, 1, name='KA', dimnames=list(LatExo, "KAMeans")) #DONE
151
152 al <- mxMatrix("Zero", numLatEnd, 1, name='AL', dimnames=list(LatEnd, "ALMeans")) #DONE
153
154
155
156 #--------------------------------------------------------------------
157 # Define Endogenous-only the model
158
159 ymod <- mxModel(
160         name='LISREL Endogenous Model with Means',
161         mxData(observed=rawlisy, type='raw'),
162         ly, be, ps, te, ty, al,
163         mxExpectationLISREL(
164                 LY=ly@name,
165                 BE=be@name,
166                 PS=ps@name,
167                 TE=te@name,
168                 TY=ty@name,
169                 AL=al@name
170         ),
171         mxFitFunctionML()
172 )
173
174
175 #--------------------------------------------------------------------
176 # Run the Endogenous-only model
177
178
179 # Uncomment the following lines when debugbing
180 #ymodRun <- mxRun(ymod, onlyFrontend=TRUE) # This runs fine.
181 #ymod <- mxOption(ymod, "Calculate Hessian", "No")
182 #ymod <- mxOption(ymod, "Standard Errors", "No")
183 #ymod <- mxOption(ymod, "Major iterations", 1)
184
185
186 ymodRun <- mxRun(ymod)
187 summary(ymodRun)
188
189
190
191 #--------------------------------------------------------------------
192 #--------------------------------------------------------------------
193 # Define the model
194
195 xmod <- mxModel(
196         name='LISREL Exogenous Model with Means',
197         mxData(observed=rawlisx, type='raw'),
198         lx, ph, td, tx, ka,
199         mxExpectationLISREL(
200                 LX=lx@name,
201                 PH=ph@name,
202                 TD=td@name,
203                 TX=tx@name,
204                 KA=ka@name
205         ),
206         mxFitFunctionML()
207 )
208
209
210 #--------------------------------------------------------------------
211 # Run the Exogenous-only model
212
213
214 # Uncomment the following lines when debugbing
215 #xmodRun <- mxRun(xmod, onlyFrontend=TRUE) # This runs fine.
216 #xmod <- mxOption(xmod, "Calculate Hessian", "No")
217 #xmod <- mxOption(xmod, "Standard Errors", "No")
218 #xmod <- mxOption(xmod, "Major iterations", 1)
219
220
221 xmodRun <- mxRun(xmod)
222 summary(xmodRun)
223
224
225
226
227
228 #--------------------------------------------------------------------
229 # Create RAM models that mirror the LISREL ones
230
231 xman <- names(rawlisx)
232
233 xrmod <- mxModel(
234         mxData(observed=rawlisx, type='raw'),
235         name="RAM Exogenous",
236         type="RAM",
237         manifestVars=xman,
238         latentVars=c("f1", "f2"),
239         mxPath("f1", xman[1:4], free=c(F, T, T, T), values=c(1, .2, .2, .2)),
240         mxPath("f2", xman[5:8], free=c(F, T, T, T), values=c(1, .2, .2, .2)),
241         mxPath(c("f1", "f2"), connect="unique.pairs", free=T, arrows=2, values=c(.8, .3, .8)),
242         mxPath(xman, arrows=2, free=T, values=.8),
243         mxPath("one", c(xman, "f1", "f2"), free=c(rep(T, 8), F, F), values=c(rep(.1, 8), 0, 0))
244 )
245
246 xrmodRun <- mxRun(xrmod)
247
248
249
250
251
252 yman <- names(rawlisy)
253
254 yrmod <- mxModel(
255         mxData(observed=rawlisy, type='raw'),
256         name="RAM Endogenous",
257         type="RAM",
258         manifestVars=yman,
259         latentVars="f1",
260         mxPath("f1", yman[1:4], free=c(F, T, T, T), values=c(1, .2, .2, .2)),
261         mxPath("f1", free=T, arrows=2, values=.8),
262         mxPath(yman, arrows=2, free=T, values=.8),
263         mxPath("one", c(yman, "f1"), free=c(rep(T, 4), F), values=c(rep(.1, 4), 0))
264 )
265
266 yrmodRun <- mxRun(yrmod)
267 summary(yrmodRun)
268
269
270
271
272 #--------------------------------------------------------------------
273 # Compare the estimate parameters in LISREL and RAM
274
275 # Check the exogenous only model
276 omxCheckCloseEnough(xmodRun@output$estimate, xrmodRun@output$estimate[c(1:6, 15:17, 7:14, 18:25)], epsilon=0.001)
277
278
279 # Check the endogenoug only model
280 omxCheckCloseEnough(ymodRun@output$estimate, yrmodRun@output$estimate[c(1:3, 8, 4:7, 9:12)], epsilon=0.001)
281
282
283
284 #--------------------------------------------------------------------
285 #require(rbenchmark)
286
287
288 #benchmark(mxRun(xmod), mxRun(xrmod), replications=10)
289 # LISREL and RAM models here take about the same amount of time
290 # LISREL is acutally .9% faster
291
292 #benchmark(mxRun(ymod), mxRun(yrmod), replications=10)
293 # LISREL takes about 9% longer than RAM models here
294
295 # For the combined model, LISREL is about 10% slower.
296
297
298
299 #--------------------------------------------------------------------
300 # End
301