Corrected LL tolerances to be within .01% of MX values. Now passes.
[openmx:openmx.git] / demo / BivariateSaturated.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 # Program: BivariateSaturated.R  
18 # Author: Hermine Maes
19 # Date: 2009.08.01 
20 #
21 # ModelType: Saturated
22 # DataType: Continuous
23 # Field: None
24 #
25 # Purpose: 
26 #      Bivariate Saturated model to estimate means and (co)variances 
27 #      using Cholesky Decomposition
28 #      Two matrix styles - Two data styles
29 #
30 # RevisionHistory:
31 #      Hermine Maes -- 2009.10.08 updated & reformatted
32 #      Ross Gore -- 2011.06.15 added Model, Data & Field metadata
33 # -----------------------------------------------------------------------
34
35 require(OpenMx)
36 require(MASS)
37 # Load Libraries
38 # -----------------------------------------------------------------------------
39
40 set.seed(200)
41 rs=.5
42 xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
43 testData <- xy
44 testData <- testData[, order(apply(testData, 2, var))[2:1]] #put the data columns in order from largest to smallest variance
45 # Note: Users do NOT have to re-order their data columns.  This is only to make data generation the same on different operating systems: to fix an inconsistency with the mvrnorm function in the MASS package.
46 selVars <- c('X','Y')
47 dimnames(testData) <- list(NULL, selVars)
48 summary(testData)
49 cov(testData)
50 # Simulate Data
51 # -----------------------------------------------------------------------------
52
53 bivSatModel1 <- mxModel("bivSat1",
54     manifestVars= selVars,
55     mxPath(
56         from=c("X", "Y"), 
57         arrows=2, 
58         free=T, 
59         values=1, 
60         lbound=.01, 
61         labels=c("varX","varY")
62     ),
63     mxPath(
64         from="X", 
65         to="Y", 
66         arrows=2, 
67         free=T, 
68         values=.2, 
69         lbound=.01, 
70         labels="covXY"
71     ),
72     mxData(
73         observed=cov(testData), 
74         type="cov", 
75         numObs=1000 
76     ),
77     type="RAM"
78     )
79 bivSatFit1 <- mxRun(bivSatModel1)
80 EC1 <- mxEval(S, bivSatFit1)
81 LL1 <- mxEval(fitfunction, bivSatFit1)
82 SL1 <- summary(bivSatFit1)$SaturatedLikelihood
83 Chi1 <- LL1-SL1
84 # example 1: Saturated Model with Cov Matrices and Path-Style Input
85 # -----------------------------------------------------------------------------
86
87 bivSatModel1m <- mxModel("bivSat1m",
88     manifestVars= selVars,
89     mxPath(
90         from=c("X", "Y"), 
91         arrows=2, 
92         free=T, 
93         values=1, 
94         lbound=.01, 
95         labels=c("varX","varY")
96     ),
97     mxPath(
98         from="X", 
99         to="Y", 
100         arrows=2, 
101         free=T, 
102         values=.2, 
103         lbound=.01, 
104         labels="covXY"
105     ),
106     mxPath(
107         from="one", 
108         to=c("X", "Y"), 
109         arrows=1, 
110         free=T, 
111         values=.01, 
112         labels=c("meanX","meanY")
113     ),
114     mxData(
115         observed=cov(testData), 
116         type="cov", 
117         numObs=1000, 
118         means=colMeans(testData)
119     ),
120     type="RAM"
121     )
122 bivSatFit1m <- mxRun(bivSatModel1m)
123 EM1m <- mxEval(M, bivSatFit1m)
124 EC1m <- mxEval(S, bivSatFit1m)
125 LL1m <- mxEval(fitfunction, bivSatFit1m)
126 SL1m <- summary(bivSatFit1m)$SaturatedLikelihood
127 Chi1m <- LL1m-SL1m
128 # example 1m: Saturated Model with Cov Matrices & Means and Path-Style Input
129 # -----------------------------------------------------------------------------
130
131 bivSatModel2 <- mxModel("bivSat2",
132     manifestVars= selVars,
133     mxPath(
134         from=c("X", "Y"), 
135         arrows=2, 
136         free=T, 
137         values=1, 
138         lbound=.01, 
139         labels=c("varX","varY")
140     ),
141     mxPath(
142         from="X", 
143         to="Y", 
144         arrows=2, 
145         free=T, 
146         values=.2, 
147         lbound=.01, 
148         labels="covXY"
149     ),
150     mxPath(from="one", 
151         to=c("X", "Y"),
152         arrows=1, 
153         free=T), 
154     mxData(
155         observed=testData, 
156         type="raw", 
157     ),
158     type="RAM"
159     )
160 bivSatFit2 <- mxRun(bivSatModel2)
161 EM2 <- mxEval(M, bivSatFit2)
162 EC2 <- mxEval(S, bivSatFit2)
163 LL2 <- mxEval(fitfunction, bivSatFit2)
164 SL2 <- summary(bivSatFit1)$SaturatedLikelihood
165 Chi2 <- LL2-SL2
166 # example 2: Saturated Model with Raw Data and Path input
167 # -----------------------------------------------------------------------------
168
169 bivSatModel2s <- mxModel(bivSatModel1,
170     mxData(
171         observed=testData, 
172         type="raw", 
173     ),
174     mxPath(from="one", 
175         to=c("X", "Y"),
176         arrows=1, 
177         free=T),     
178     name = "bivSat2s"
179     )
180 bivSatFit2s <- mxRun(bivSatModel2s)
181 EM2s <- mxEval(M, bivSatFit2s)
182 EC2s <- mxEval(S, bivSatFit2s)
183 LL2s <- mxEval(fitfunction, bivSatFit2s)
184 # example 2s: Saturated Model with Raw Data and Path input built upon Cov/Means version
185 # -----------------------------------------------------------------------------
186
187 bivSatModel3 <- mxModel("bivSat3",
188     mxMatrix(
189         type="Symm", 
190         nrow=2, 
191         ncol=2, 
192         free=T, 
193         values=c(1,.5,1), 
194         dimnames=list(selVars,selVars), 
195         name="expCov"
196     ),
197     mxData(
198         observed=cov(testData), 
199         type="cov", 
200         numObs=1000 
201     ),
202     mxExpectationNormal(
203         covariance="expCov"
204     ),
205     mxFitFunctionML()
206     )
207 bivSatFit3 <- mxRun(bivSatModel3)
208 EC3 <- mxEval(expCov, bivSatFit3)
209 LL3 <- mxEval(fitfunction,bivSatFit3)
210 SL3 <- summary(bivSatFit3)$SaturatedLikelihood
211 Chi3 <- LL3-SL3
212 # example 3: Saturated Model with Cov Matrices and Matrix-Style Input
213 # -----------------------------------------------------------------------------
214
215 bivSatModel3m <- mxModel("bivSat3m",
216     mxMatrix(
217         type="Symm", 
218         nrow=2, 
219         ncol=2, 
220         free=T, 
221         values=c(1,.5,1), 
222         name="expCov"
223     ),
224     mxMatrix(
225         type="Full", 
226         nrow=1, 
227         ncol=2, 
228         free=T, 
229         values=c(0,0), 
230         name="expMean"
231     ),
232     mxData(
233         observed=cov(testData), 
234         type="cov", 
235         numObs=1000, 
236         means=colMeans(testData) 
237     ),
238     mxExpectationNormal(
239         covariance="expCov",
240         means="expMean",
241         dimnames=selVars
242     ),
243     mxFitFunctionML()
244     )
245 bivSatFit3m <- mxRun(bivSatModel3m)
246 EM3m <- mxEval(expMean, bivSatFit3m)
247 EC3m <- mxEval(expCov, bivSatFit3m)
248 LL3m <- mxEval(fitfunction,bivSatFit3m)
249 SL3m <- summary(bivSatFit3m)$SaturatedLikelihood
250 Chi3m <- LL3m-SL3m
251 # example 3m: Saturated Model with Cov Matrices & Means and Matrix-Style Input
252 # -----------------------------------------------------------------------------
253
254 bivSatModel4 <- mxModel("bivSat4",
255 mxMatrix(
256     type="Symm", 
257     nrow=2, 
258     ncol=2, 
259     free=T, 
260     values=c(1,.5,1), 
261     name="expCov"
262 ),
263 mxMatrix(
264     type="Full", 
265     nrow=1, 
266     ncol=2, 
267     free=T, 
268     values=c(0,0), 
269     name="expMean"
270 ),
271 mxData(
272     observed=testData, 
273     type="raw", 
274 ),
275 mxExpectationNormal(
276     covariance="expCov",
277     means="expMean",
278     dimnames=selVars
279 ),
280 mxFitFunctionML()
281 )
282 bivSatFit4 <- mxRun(bivSatModel4)
283 EM4 <- mxEval(expMean, bivSatFit4)
284 EC4 <- mxEval(expCov, bivSatFit4)
285 LL4 <- mxEval(fitfunction,bivSatFit4)
286 # examples 4: Saturated Model with Raw Data and Matrix-Style Input
287 # -----------------------------------------------------------------------------
288
289 bivSatModel5 <- mxModel("bivSat5",
290     mxMatrix(
291         type="Full", 
292         nrow=2, 
293         ncol=2, 
294         free=c(T,T,F,T), 
295         values=c(1,.2,0,1), 
296         name="Chol"
297     ),
298     mxAlgebra(
299         Chol %*% t(Chol), 
300         name="expCov", 
301         dimnames=list(selVars,selVars)
302     ),
303     mxMatrix(
304         type="Full", 
305         nrow=1, 
306         ncol=2, 
307         free=T, 
308         values=c(0,0), 
309         dimnames=list(NULL,selVars), 
310         name="expMean"
311     ),
312     mxData(
313         observed=cov(testData), 
314         type="cov", 
315         numObs=1000 
316     ),
317     mxExpectationNormal(
318         covariance="expCov"
319     ),
320     mxFitFunctionML()
321     )
322 bivSatFit5 <- mxRun(bivSatModel5)
323 EC5 <- mxEval(expCov, bivSatFit5)
324 LL5 <- mxEval(fitfunction,bivSatFit5)
325 SL5 <- summary(bivSatFit5)$SaturatedLikelihood
326 Chi5 <- LL5-SL5
327 # example 5: Saturated Model with Cov Matrices and Matrix-Style Input
328 # -----------------------------------------------------------------------------
329
330 bivSatModel5m <- mxModel("bivSat5m",
331     mxMatrix(
332         type="Full", 
333         nrow=2, 
334         ncol=2, 
335         free=c(T,T,F,T), 
336         values=c(1,.2,0,1), 
337         name="Chol"
338     ),
339     mxAlgebra(
340         Chol %*% t(Chol), 
341         name="expCov", 
342         dimnames=list(selVars,selVars)
343     ),
344     mxMatrix(
345         type="Full", 
346         nrow=1, 
347         ncol=2, 
348         free=T, 
349         values=c(0,0), 
350         dimnames=list(NULL,selVars), 
351         name="expMean"
352     ),
353     mxData(
354         observed=cov(testData), 
355         type="cov", 
356         numObs=1000, 
357         means=colMeans(testData) 
358     ),
359     mxExpectationNormal(
360         covariance="expCov",
361         means="expMean"
362     ),
363     mxFitFunctionML()
364     )
365 bivSatFit5m <- mxRun(bivSatModel5m)
366 EM5m <- mxEval(expMean, bivSatFit5m)
367 EC5m <- mxEval(expCov, bivSatFit5m)
368 LL5m <- mxEval(fitfunction,bivSatFit5m);
369 SL5m <- summary(bivSatFit5m)$SaturatedLikelihood
370 Chi5m <- LL5m-SL5m
371 # example 5m: Saturated Model with Cov Matrices & Means and Matrix-Style Input
372 # -----------------------------------------------------------------------------
373
374 bivSatModel6 <- mxModel("bivSat6",
375     mxMatrix(
376         type="Full", 
377         nrow=2, 
378         ncol=2, 
379         free=c(T,T,F,T), 
380         values=c(1,.2,0,1), 
381         name="Chol"
382     ),
383     mxAlgebra(
384         Chol %*% t(Chol), 
385         name="expCov", 
386         dimnames=list(selVars,selVars)
387     ),
388     mxMatrix(
389         type="Full", 
390         nrow=1, 
391         ncol=2, 
392         free=T, 
393         values=c(0,0), 
394         dimnames=list(NULL,selVars), 
395         name="expMean"
396     ),
397     mxData(
398         observed=testData, 
399         type="raw", 
400     ),
401     mxExpectationNormal(
402         covariance="expCov",
403         means="expMean"
404     ),
405     mxFitFunctionML()
406     )
407 bivSatFit6 <- mxRun(bivSatModel6)
408 EM6 <- mxEval(expMean, bivSatFit6)
409 EC6 <- mxEval(expCov, bivSatFit6)
410 LL6 <- mxEval(fitfunction,bivSatFit6)
411 # examples 6: Saturated Model with RawData and Matrix-Style Input
412 # -----------------------------------------------------------------------------
413
414
415
416 Mx.EC1 <- matrix(c(1.0102951, 0.4818317, 0.4818317, 0.9945329),2,2)
417 Mx.LL1 <- -2.258885e-13
418 # example Mx..1: Saturated Model 
419 # with Cov Matrices
420 # -------------------------------------
421
422 Mx.EM1m <- matrix(c(0.03211648, -0.004883811),1,2)
423 Mx.EC1m <- matrix(c(1.0102951, 0.4818317, 0.4818317, 0.9945329),2,2)
424 Mx.LL1m <- -5.828112e-14
425 # example Mx..1m: Saturated Model 
426 # with Cov Matrices & Means
427 # -------------------------------------
428
429 Mx.EM2 <- matrix(c(0.03211188, -0.004889211),1,2)
430 Mx.EC2 <- matrix(c(1.0092891, 0.4813504, 0.4813504, 0.9935366),2,2)
431 Mx.LL2 <- 5415.772
432 # example Mx..2: Saturated Model with
433 # Raw Data
434 # -------------------------------------
435 # Mx answers hard-coded
436 # -----------------------------------------------------------------------------
437
438 cov <- rbind(cbind(EC1,EC1m,EC2),cbind(EC3,EC3m,EC4))
439 mean <- rbind(cbind(EM1m, EM2),cbind(EM3m,EM4))
440 like <- rbind(cbind(LL1,LL1m,LL2),cbind(LL3,LL3m,LL4))
441 cov; mean; like
442 # OpenMx summary
443 # -----------------------------------------------------------------------
444
445
446 Mx.cov <- cbind(Mx.EC1,Mx.EC1m,Mx.EC2)
447 Mx.mean <- cbind(Mx.EM1m,Mx.EM2)
448 Mx.like <- cbind(Mx.LL1,Mx.LL1m,Mx.LL2)
449 Mx.cov; Mx.mean; Mx.like
450 # old Mx summary
451 # -----------------------------------------------------------------------
452
453
454 omxCheckCloseEnough(Chi1,Mx.LL1,.01)
455 omxCheckCloseEnough(EC1,Mx.EC1,.001)
456 #1:CovPat
457 # -------------------------------------
458
459 omxCheckCloseEnough(Chi1m,Mx.LL1m,.01)
460 omxCheckCloseEnough(EC1m,Mx.EC1m,.001)
461 omxCheckCloseEnough(EM1m,Mx.EM1m,.001)
462 #1m:CovMPat 
463 # -------------------------------------
464
465 omxCheckCloseEnough(LL2,Mx.LL2,.01)
466 omxCheckCloseEnough(EC2,Mx.EC2,.001)
467 omxCheckCloseEnough(EM2,Mx.EM2,.001)
468 #2:RawPat 
469 # -------------------------------------
470
471 omxCheckCloseEnough(LL2s,Mx.LL2,.01)
472 omxCheckCloseEnough(EC2s,Mx.EC2,.001)
473 omxCheckCloseEnough(EM2s,Mx.EM2,.001)
474 #2:RawSPat
475 # -------------------------------------
476
477 omxCheckCloseEnough(Chi3,Mx.LL1,.01)
478 omxCheckCloseEnough(EC3,Mx.EC1,.001)
479 #3:CovMat
480 # -------------------------------------
481
482 omxCheckCloseEnough(Chi3m,Mx.LL1m,.01)
483 omxCheckCloseEnough(EC3m,Mx.EC1m,.001)
484 omxCheckCloseEnough(EM3m,Mx.EM1m,.001)
485 #3m:CovMPat 
486 # -------------------------------------
487
488 omxCheckCloseEnough(LL4,Mx.LL2,.01)
489 omxCheckCloseEnough(EC4,Mx.EC2,.001)
490 omxCheckCloseEnough(EM4,Mx.EM2,.001)
491 #4:RawMat
492 # -------------------------------------
493
494 omxCheckCloseEnough(Chi5,Mx.LL1,.01)
495 omxCheckCloseEnough(EC5,Mx.EC1,.001)
496 #5:CovMat Cholesky
497 # -------------------------------------
498
499 omxCheckCloseEnough(Chi5m,Mx.LL1m,.01)
500 omxCheckCloseEnough(EC5m,Mx.EC1m,.001)
501 omxCheckCloseEnough(EM5m,Mx.EM1m,.001)
502 #5m:CovMPat Cholesky
503 # -------------------------------------
504
505
506 omxCheckCloseEnough(LL6,Mx.LL2,.01)
507 omxCheckCloseEnough(EC6,Mx.EC2,.001)
508 omxCheckCloseEnough(EM6,Mx.EM2,.001)
509 #6:RawMat Cholesky
510 # -------------------------------------
511 # Compare OpenMx results to Mx results 
512 # (LL: likelihood; EC: expected covariance, EM: expected means)
513 # -----------------------------------------------------------------------
514