Convert mxFIMLObjective, mxMLObjective, and mxRAMObjective
[openmx:openmx.git] / models / passing / univSatM.R
1 #
2 #   Copyright 2007-2013 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 require(OpenMx)
18
19 #Simulate Data
20 set.seed(100)
21 x <- rnorm (1000, 0, 1)
22 testData <- as.matrix(x)
23 selVars <- c("X")
24 dimnames(testData) <- list(NULL, selVars)
25 summary(testData)
26 colMeans(testData)
27 var(testData)
28
29 #example 1: Saturated Model with Cov Matrices and Paths input
30 univSatModel1 <- mxModel("univSat1",
31     manifestVars= selVars,
32     mxPath(
33         from=c("X"), 
34         arrows=2, 
35         free=T, 
36         values=1, 
37         lbound=.01, 
38         labels="vX"
39     ),
40     mxData(
41         observed=var(testData), 
42         type="cov", 
43         numObs=1000 
44     ),
45     type="RAM"
46     )
47 univSatFit1 <- mxRun(univSatModel1)
48 EC1 <- mxEval(S, univSatFit1)
49 LL1 <- mxEval(objective, univSatFit1)
50 SL1 <- univSatFit1@output$SaturatedLikelihood
51 Chi1 <- LL1-SL1
52
53 #example 1m: Saturated Model with Cov Matrices & Means and Paths input
54 univSatModel1m <- mxModel("univSat1m",
55     manifestVars= selVars,
56     mxPath(
57         from=c("X"), 
58         arrows=2, 
59         free=T, 
60         values=1, 
61         lbound=.01, 
62         labels="vX"
63     ),
64     mxPath(
65         from="one", 
66         to="X", 
67         arrows=1, 
68         free=T, 
69         values=0, 
70         labels="mX"
71     ),
72     mxData(
73         observed=var(testData), 
74         type="cov", 
75         numObs=1000, 
76         means=colMeans(testData)
77     ),
78     type="RAM"
79     )
80 univSatFit1m <- mxRun(univSatModel1m)
81 EM1m <- mxEval(M, univSatFit1m)
82 EC1m <- mxEval(S, univSatFit1m)
83 LL1m <- mxEval(objective,univSatFit1m);
84 SL1m <- univSatFit1m@output$SaturatedLikelihood
85 Chi1m <- LL1m-SL1m
86
87 #example 2: Saturated Model with Raw Data and Path input
88 univSatModel2 <- mxModel("univSat2",
89     manifestVars= selVars,
90     mxPath(
91         from=c("X"), 
92         arrows=2, 
93         free=T, 
94         values=1, 
95         lbound=.01, 
96         labels="vX"
97     ),
98     mxPath(
99         from="one", 
100         to="X", 
101         arrows=1, 
102         free=T, 
103         values=0, 
104         labels="mX"
105     ),    
106     mxData(
107         observed=testData, 
108         type="raw", 
109     ),
110     type="RAM"
111     )
112 univSatFit2 <- mxRun(univSatModel2)
113 EM2 <- mxEval(M, univSatFit2)
114 EC2 <- mxEval(S, univSatFit2)
115 LL2 <- mxEval(objective,univSatFit2);
116
117 #example 2s: Saturated Model with Raw Data and Path input built upon Cov/Means version
118 univSatModel2s <- mxModel(univSatModel1,
119     mxData(
120         observed=testData, 
121         type="raw"
122     ),
123     mxPath(
124         from="one", 
125         to="X", 
126         arrows=1, 
127         free=T, 
128         values=0, 
129         labels="mX"
130     ),        
131     name="univSat2s",
132     type="RAM"
133     )
134 univSatFit2s <- mxRun(univSatModel2s)
135 EM2s <- mxEval(M, univSatFit2s)
136 EC2s <- mxEval(S, univSatFit2s)
137 LL2s <- mxEval(objective,univSatFit2s);
138
139
140 #example 3: Saturated Model with Cov Matrices and Matrices input
141 univSatModel3 <- mxModel("univSat3",
142     mxMatrix(
143         type="Symm", 
144         nrow=1, 
145         ncol=1, 
146         free=T, 
147         values=1, 
148         dimnames=list(selVars,selVars), 
149         name="expCov"
150     ),
151     mxData(
152         observed=var(testData), 
153         type="cov", 
154         numObs=1000 
155     ),
156     mxFitFunctionML(),mxExpectationNormal(
157         covariance="expCov"
158     )
159     )
160 univSatFit3 <- mxRun(univSatModel3)
161 EC3 <- mxEval(expCov, univSatFit3)
162 LL3 <- mxEval(objective, univSatFit3)
163 SL3 <- univSatFit3@output$SaturatedLikelihood
164 Chi3 <- LL3-SL3
165
166 #example 3m: Saturated Model with Cov Matrices & Means and Matrices input
167 univSatModel3m <- mxModel("univSat3m",
168     mxMatrix(
169         type="Symm", 
170         nrow=1, 
171         ncol=1, 
172         free=T, 
173         values=1, 
174         dimnames=list(selVars,selVars), 
175         name="expCov"
176     ),
177     mxMatrix(
178         type="Full", 
179         nrow=1, 
180         ncol=1, 
181         free=T, 
182         values=0, 
183         dimnames=list(NULL,selVars), 
184         name="expMean"
185     ),
186     mxData(
187         observed=var(testData), 
188         type="cov", 
189         numObs=1000,
190         means=colMeans(testData) 
191     ),
192     mxFitFunctionML(),mxExpectationNormal(
193         covariance="expCov", 
194         means="expMean"
195     )
196     )
197 univSatFit3m <- mxRun(univSatModel3m)
198 EM3m <- mxEval(expMean, univSatFit3m)
199 EC3m <- mxEval(expCov, univSatFit3m)
200 LL3m <- mxEval(objective, univSatFit3m);
201 SL3m <- univSatFit3m@output$SaturatedLikelihood
202 Chi3m <- LL3m-SL3m
203
204 #examples 4: Saturated Model with Raw Data and Matrices input
205 univSatModel4 <- mxModel("univSat4",
206     mxMatrix(
207         type="Symm", 
208         nrow=1, 
209         ncol=1, 
210         free=T, 
211         values=1, 
212         dimnames=list(selVars,selVars), 
213         name="expCov"
214     ),
215     mxMatrix(
216         type="Full", 
217         nrow=1, 
218         ncol=1, 
219         free=T, 
220         values=0, 
221         dimnames=list(NULL,selVars), 
222         name="expMean"
223     ),
224     mxData(
225         observed=testData, 
226         type="raw", 
227     ),
228     mxFitFunctionML(),mxExpectationNormal(
229         covariance="expCov", 
230         means="expMean"
231     )
232     )
233 univSatFit4 <- mxRun(univSatModel4)
234 EM4 <- mxEval(expMean, univSatFit4)
235 EC4 <- mxEval(expCov, univSatFit4)
236 LL4 <- mxEval(objective, univSatFit4);
237
238
239 #Mx answers hard-coded
240 #example Mx..1: Saturated Model with Cov Matrices
241 Mx.EC1 <-  1.062112
242 Mx.LL1 <- -1.474434e-17
243
244 #example Mx..1m: Saturated Model with Cov Matrices & Means
245 Mx.EM1m <- 0.01680509
246 Mx.EC1m <- 1.062112
247 Mx.LL1m <- -1.108815e-13
248
249 Mx.EM2 <- 0.01680516
250 Mx.EC2 <- 1.061050
251 Mx.LL2 <- 2897.135
252
253
254 #OpenMx summary
255 cov <- rbind(cbind(EC1,EC1m,EC2),cbind(EC3,EC3m,EC4))
256 mean <- rbind(cbind(EM1m, EM2),cbind(EM3m,EM4))
257 like <- rbind(cbind(LL1,LL1m,LL2),cbind(LL3,LL3m,LL4))
258 cov; mean; like
259
260 #old Mx summary
261 Mx.cov <- cbind(Mx.EC1,Mx.EC1m,Mx.EC2)
262 Mx.mean <- cbind(Mx.EM1m,Mx.EM2)
263 Mx.like <- cbind(Mx.LL1,Mx.LL1m,Mx.LL2)
264 Mx.cov; Mx.mean; Mx.like
265
266
267 #Compare OpenMx results to Mx results (LL: likelihood; EC: expected covariance, EM: expected means)
268 #1:CovPat
269 omxCheckCloseEnough(Chi1,Mx.LL1,.001)
270 omxCheckCloseEnough(EC1,Mx.EC1,.001)
271 #1m:CovMPat 
272 omxCheckCloseEnough(Chi1m,Mx.LL1m,.001)
273 omxCheckCloseEnough(EC1m,Mx.EC1m,.001)
274 omxCheckCloseEnough(EM1m,Mx.EM1m,.001)
275 #2:RawPat 
276 omxCheckCloseEnough(LL2,Mx.LL2,.001)
277 omxCheckCloseEnough(EC2,Mx.EC2,.001)
278 omxCheckCloseEnough(EM2,Mx.EM2,.001)
279 #2:RawSPat
280 omxCheckCloseEnough(LL2s,Mx.LL2,.001)
281 omxCheckCloseEnough(EC2s,Mx.EC2,.001)
282 omxCheckCloseEnough(EM2s,Mx.EM2,.001)
283 #3:CovMat
284 omxCheckCloseEnough(Chi3,Mx.LL1,.001)
285 omxCheckCloseEnough(EC3,Mx.EC1,.001)
286 #3m:CovMPat 
287 omxCheckCloseEnough(Chi3m,Mx.LL1m,.001)
288 omxCheckCloseEnough(EC3m,Mx.EC1m,.001)
289 omxCheckCloseEnough(EM3m,Mx.EM1m,.001)
290 #4:RawMat
291 omxCheckCloseEnough(LL4,Mx.LL2,.001)
292 omxCheckCloseEnough(EC4,Mx.EC2,.001)
293 omxCheckCloseEnough(EM4,Mx.EM2,.001)
294