(Hopefully) fixed the error with the numerous bivariate saturated models. I believe...
[openmx:openmx.git] / demo / BivariateSaturated_MatrixRawCholesky.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 # -----------------------------------------------------------------------------
18 # Program: BivariateSaturated_MatrixRawCholesky.R  
19 # Author: Hermine Maes
20 # Date: 2009.08.01 
21 #
22 # ModelType: Saturated
23 # DataType: Continuous
24 # Field: None
25 #
26 # Purpose: 
27 #      Bivariate Saturated model to estimate means and (co)variances 
28 #      using Cholesky Decomposition
29 #      Matrix style model input - Raw data input
30 #
31 # RevisionHistory:
32 #      Hermine Maes -- 2009.10.08 updated & reformatted
33 #      Ross Gore -- 2011.06.15 added Model, Data & Field metadata     
34 # -----------------------------------------------------------------------------
35
36 require(OpenMx)
37 require(MASS)
38 # Load Libraries
39 # -----------------------------------------------------------------------------
40
41 set.seed(200)
42 rs=.5
43 xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
44 testData <- xy
45 testData <- testData[, order(apply(testData, 2, var))[2:1]] #put the data columns in order from largest to smallest variance
46 selVars <- c("X","Y")
47 dimnames(testData) <- list(NULL, selVars)
48 summary(testData)
49 cov(testData)
50 # Simulate Data
51 # -----------------------------------------------------------------------------
52
53 bivSatModel6 <- mxModel("bivSat6",
54     mxMatrix(
55                 type="Full", 
56                 nrow=2, 
57                 ncol=2, 
58                 free=c(TRUE,TRUE,FALSE,TRUE), 
59                 values=c(1,.2,0,1), 
60                 name="Chol"
61         ),
62         mxMatrix(
63                 type="Full", 
64                 nrow=1, 
65                 ncol=2, 
66                 free=TRUE, 
67                 values=c(0,0), 
68                 name="expMean"
69         ),
70     mxAlgebra(
71                 expression=Chol %*% t(Chol), 
72                 name="expCov"
73         ),
74         mxData(
75                 observed=testData, 
76                 type="raw"
77         ),
78         mxFitFunctionML(),mxExpectationNormal(
79                 covariance="expCov",
80                 means="expMean",
81                 dimnames=selVars
82         )
83 )
84
85 bivSatFit6 <- mxRun(bivSatModel6)
86 EM <- mxEval(expMean, bivSatFit6)
87 EC <- mxEval(expCov, bivSatFit6)
88 LL <- mxEval(objective,bivSatFit6)
89 # example 6: Saturated Model with Raw Data and Matrix-Style Input
90 # -----------------------------------------------------------------------------
91
92 Mx.EM <- matrix(c(0.03211188, -0.004889211),1,2)
93 Mx.EC <- matrix(c(1.0092891, 0.4813504, 0.4813504, 0.9935366),2,2)
94 Mx.LL <- 5415.772
95 # Mx answers hard-coded
96 # -------------------------------------
97 # example Mx..2: Saturated Model with Raw Data
98 # -----------------------------------------------------------------------------
99
100 omxCheckCloseEnough(LL,Mx.LL,.001)
101 omxCheckCloseEnough(EC,Mx.EC,.001)
102 omxCheckCloseEnough(EM,Mx.EM,.001)
103 # 6: RawMat Cholesky
104 # -------------------------------------
105 # Compare OpenMx results to Mx results 
106 # (LL: likelihood; EC: expected covariance, EM: expected means)
107 # -----------------------------------------------------------------------------