Made a note to users in the bivariate saturated files so they don't think they have...
[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 # 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.
47 selVars <- c("X","Y")
48 dimnames(testData) <- list(NULL, selVars)
49 summary(testData)
50 cov(testData)
51 # Simulate Data
52 # -----------------------------------------------------------------------------
53
54 bivSatModel6 <- mxModel("bivSat6",
55     mxMatrix(
56                 type="Full", 
57                 nrow=2, 
58                 ncol=2, 
59                 free=c(TRUE,TRUE,FALSE,TRUE), 
60                 values=c(1,.2,0,1), 
61                 name="Chol"
62         ),
63         mxMatrix(
64                 type="Full", 
65                 nrow=1, 
66                 ncol=2, 
67                 free=TRUE, 
68                 values=c(0,0), 
69                 name="expMean"
70         ),
71     mxAlgebra(
72                 expression=Chol %*% t(Chol), 
73                 name="expCov"
74         ),
75         mxData(
76                 observed=testData, 
77                 type="raw"
78         ),
79         mxFitFunctionML(),mxExpectationNormal(
80                 covariance="expCov",
81                 means="expMean",
82                 dimnames=selVars
83         )
84 )
85
86 bivSatFit6 <- mxRun(bivSatModel6)
87 EM <- mxEval(expMean, bivSatFit6)
88 EC <- mxEval(expCov, bivSatFit6)
89 LL <- mxEval(objective,bivSatFit6)
90 # example 6: Saturated Model with Raw Data and Matrix-Style Input
91 # -----------------------------------------------------------------------------
92
93 Mx.EM <- matrix(c(0.03211188, -0.004889211),1,2)
94 Mx.EC <- matrix(c(1.0092891, 0.4813504, 0.4813504, 0.9935366),2,2)
95 Mx.LL <- 5415.772
96 # Mx answers hard-coded
97 # -------------------------------------
98 # example Mx..2: Saturated Model with Raw Data
99 # -----------------------------------------------------------------------------
100
101 omxCheckCloseEnough(LL,Mx.LL,.001)
102 omxCheckCloseEnough(EC,Mx.EC,.001)
103 omxCheckCloseEnough(EM,Mx.EM,.001)
104 # 6: RawMat Cholesky
105 # -------------------------------------
106 # Compare OpenMx results to Mx results 
107 # (LL: likelihood; EC: expected covariance, EM: expected means)
108 # -----------------------------------------------------------------------------