(Hopefully) fixed the error with the numerous bivariate saturated models. I believe...
[openmx:openmx.git] / demo / BivariateSaturated_PathRaw.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 # Program: BivariateSaturated_PathRaw.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 #      Path style model input - Raw data input
28 #
29 # RevisionHistory:
30 #      Hermine Maes -- 2009.10.08 updated & reformatted
31 #      Ross Gore -- 2011.06.15 added Model, Data & Field metadata
32 # -----------------------------------------------------------------------------
33
34 require(OpenMx)
35 require(MASS)
36 # Load Libraries
37 # -----------------------------------------------------------------------------
38
39 set.seed(200)
40 rs=.5
41 xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
42 testData <- xy
43 testData <- testData[, order(apply(testData, 2, var))[2:1]] #put the data columns in order from largest to smallest variance
44 selVars <- c("X","Y")
45 dimnames(testData) <- list(NULL, selVars)
46 summary(testData)
47 cov(testData)
48 # Simulate Data
49 # -----------------------------------------------------------------------------
50
51 bivSatModel2 <- mxModel("bivSat2",
52     manifestVars= selVars,
53     mxPath(
54         from=c("X", "Y"), 
55         arrows=2, 
56         free=T, 
57         values=1, 
58         lbound=.01, 
59         labels=c("varX","varY")
60     ),
61     mxPath(
62         from="X", 
63         to="Y", 
64         arrows=2, 
65         free=T, 
66         values=.2, 
67         lbound=.01, 
68         labels="covXY"
69     ),
70     mxPath(
71         from="one",
72         to=c("X", "Y"),
73         arrows=1,
74         free=T),
75     mxData(
76         observed=testData, 
77         type="raw", 
78     ),
79     type="RAM"
80 )
81
82 bivSatFit2 <- mxRun(bivSatModel2)
83 EM2 <- mxEval(M, bivSatFit2)
84 EC2 <- mxEval(S, bivSatFit2)
85 LL2 <- mxEval(objective, bivSatFit2)
86 SL2 <- summary(bivSatFit2)$SaturatedLikelihood
87 Chi2 <- LL2-SL2
88 # Example 2: Saturated Model with Raw Data and Path input
89 # -----------------------------------------------------------------------------
90
91 Mx.EM2 <- matrix(c(0.03211188, -0.004889211),1,2)
92 Mx.EC2 <- matrix(c(1.0092891, 0.4813504, 0.4813504, 0.9935366),2,2)
93 Mx.LL2 <- 5415.772
94 # example Mx..2: Saturated Model with Raw Data
95 # Mx answers hard-coded
96 # -----------------------------------------------------------------------------
97
98 omxCheckCloseEnough(LL2,Mx.LL2,.001)
99 omxCheckCloseEnough(EC2,Mx.EC2,.001)
100 omxCheckCloseEnough(EM2,Mx.EM2,.001)
101 # 2:RawPat 
102 # -------------------------------------
103 # Compare OpenMx results to Mx results 
104 # (LL: likelihood; EC: expected covariance, EM: expected means)
105 # -----------------------------------------------------------------------------