(Hopefully) fixed the error with the numerous bivariate saturated models. I believe...
[openmx:openmx.git] / demo / BivariateHeterogeneity_MatrixRaw.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: BivariateHeterogeneity_MatrixRaw.R  
19 # Author: Hermine Maes
20 # Date: 2009.08.01 
21 #
22 # ModelType: Heterogeneity
23 # DataType: Continuous
24 # Field: None
25 #
26 # Purpose: 
27 #      Bivariate Heterogeneity model to test differences in means and 
28 #      variances across multiple groups
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
38 require(MASS)
39 # Load Libraries
40 # -----------------------------------------------------------------------------
41
42 set.seed(200)
43 rs=.5
44 xy1 <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
45 xy1 <- xy1[, order(apply(xy1, 2, var))[2:1]] #put the data columns in order from largest to smallest variance
46 set.seed(200)
47 # group 1
48 # --------------------------------------
49
50 rs=.4
51 xy2 <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
52 # group 2
53 # --------------------------------------
54 # Simulate Data
55 # -----------------------------------------------------------------------------
56
57
58 selVars <- c("X","Y")
59 summary(xy1)
60 cov(xy1)
61 dimnames(xy1) <- list(NULL, selVars)
62 summary(xy2)
63 cov(xy2)
64 dimnames(xy2) <- list(NULL, selVars)
65 # Print Descriptive Statistics
66 # -----------------------------------------------------------------------------
67
68 bivHetModel <- mxModel("bivariate Heterogeneity Matrix Specification",
69     mxModel("group1",
70         mxMatrix(
71             type="Lower", 
72             nrow=2, 
73             ncol=2, 
74             free=T, 
75             values=.5,
76             labels=c("vX1", "cXY1", "vY1"),
77             name="Chol1"
78         ), 
79         mxAlgebra(
80             expression=Chol1 %*% t(Chol1), 
81             name="EC1"
82         ), 
83         mxMatrix(
84             type="Full", 
85             nrow=1, 
86             ncol=2, 
87             free=T, 
88             values=c(0,0), 
89             labels=c("mX1", "mY1"), 
90             name="EM1"
91         ), 
92         mxData(
93             xy1, 
94             type="raw"
95         ), 
96         mxExpectationNormal(
97             "EC1", 
98             "EM1",
99             selVars),
100         mxFitFunctionML()
101         ),
102     mxModel("group2",
103         mxMatrix(
104             type="Lower", 
105             nrow=2, 
106             ncol=2, 
107             free=T, 
108             values=.5,
109             labels=c("vX2", "cXY2", "vY2"),
110             name="Chol2"
111         ), 
112         mxAlgebra(
113             expression=Chol2 %*% t(Chol2), 
114             name="EC2" 
115         ), 
116         mxMatrix(
117             type="Full", 
118             nrow=1, 
119             ncol=2, 
120             free=T, 
121             values=c(0,0), 
122             labels=c("mX2", "mY2"), 
123             name="EM2"
124         ), 
125         mxData(
126             xy2, 
127             type="raw"
128         ), 
129         mxExpectationNormal(
130             "EC2", 
131             "EM2",
132             selVars),
133         mxFitFunctionML()
134         ),
135     mxAlgebra(
136         group1.fitfunction + group2.fitfunction, 
137         name="h12"
138     ),
139     mxFitFunctionAlgebra("h12")
140 )
141
142 bivHetFit <- mxRun(bivHetModel)
143     EM1Het <- mxEval(group1.EM1, bivHetFit)
144     EM2Het <- mxEval(group2.EM2, bivHetFit)
145     EC1Het <- mxEval(group1.EC1, bivHetFit)
146     EC2Het <- mxEval(group2.EC2, bivHetFit)
147     LLHet <- mxEval(fitfunction, bivHetFit)
148 # Fit Heterogeneity Model
149 # -----------------------------------------------------------------------------
150
151 bivHomModel <- bivHetModel
152     bivHomModel[['group2.Chol2']]@labels <- bivHomModel[['group1.Chol1']]@labels
153     bivHomModel[['group2.EM2']]@labels <- bivHomModel[['group1.EM1']]@labels
154 bivHomFit <- mxRun(bivHomModel)
155     EM1Hom <- mxEval(group1.EM1, bivHomFit)
156     EM2Hom <- mxEval(group2.EM2, bivHomFit)
157     EC1Hom <- mxEval(group1.EC1, bivHomFit)
158     EC2Hom <- mxEval(group2.EC2, bivHomFit)
159     LLHom <- mxEval(fitfunction, bivHomFit)
160
161     Chi= LLHom-LLHet
162     LRT= rbind(LLHet,LLHom,Chi)
163     LRT
164 # Fit Homnogeneity Model
165 # -----------------------------------------------------------------------------
166
167
168
169 Mx.EM1Het <- matrix(c(0.03211284, -0.004889846),1,2)
170 Mx.EC1Het <- matrix(c(1.0092856, 0.4813512, 0.4813512, 0.9935414),2,2)
171 Mx.EM2Het <- matrix(c(0.03341992, -0.007112054),1,2)
172 Mx.EC2Het <- matrix(c(1.012324, 0.3799160, 0.379916, 0.9956605),2,2)
173 Mx.LLHet <- 10944.873
174 # 1: Heterogeneity Model
175 # -------------------------------------
176
177 Mx.EM1Hom <- matrix(c(0.03276872, -0.0059975),1,2)
178 Mx.EC1Hom <- matrix(c(1.0108055, 0.4306339, 0.4306339, 0.9946009),2,2)
179 Mx.EM2Hom <- matrix(c(0.03276872, -0.0059975),1,2)
180 Mx.EC2Hom <- matrix(c(1.0108055, 0.4306339, 0.4306339, 0.9946009),2,2)
181 Mx.LLHom <- 10954.368
182 # 2: Homogeneity Model
183 # -------------------------------------
184 # Mx answers hard-coded
185 # -----------------------------------------------------------------------------
186
187
188 cov <- rbind(cbind(EC1Het,EC2Het),cbind(EC1Hom,EC2Hom))
189 mean <- rbind(cbind(EM1Het, EM2Het),cbind(EM1Hom,EM2Hom))
190 like <- rbind(cbind(LLHet),cbind(LLHom))
191 cov; mean; like
192 # OpenMx summary
193 # -----------------------------------------------------------------------------
194
195 Mx.cov <- rbind(cbind(Mx.EC1Het,Mx.EC2Het),cbind(Mx.EC1Hom,Mx.EC2Hom))
196 Mx.mean <- rbind(cbind(Mx.EM1Het, Mx.EM2Het),cbind(Mx.EM1Hom,Mx.EM2Hom))
197 Mx.like <- rbind(cbind(Mx.LLHet),cbind(Mx.LLHom))
198 Mx.cov; Mx.mean; Mx.like
199 # old Mx summary
200 # -----------------------------------------------------------------------------
201
202
203
204 omxCheckCloseEnough(LLHet,Mx.LLHet,.001)
205 omxCheckCloseEnough(EC1Het,Mx.EC1Het,.001)
206 omxCheckCloseEnough(EM1Het,Mx.EM1Het,.001)
207 omxCheckCloseEnough(EC2Het,Mx.EC2Het,.001)
208 omxCheckCloseEnough(EM2Het,Mx.EM2Het,.001)
209
210 omxCheckCloseEnough(LLHom,Mx.LLHom,.001)
211 omxCheckCloseEnough(EC1Hom,Mx.EC1Hom,.001)
212 omxCheckCloseEnough(EM1Hom,Mx.EM1Hom,.001)
213 omxCheckCloseEnough(EC2Hom,Mx.EC2Hom,.001)
214 omxCheckCloseEnough(EM2Hom,Mx.EM2Hom,.001)
215 # Compare OpenMx results to Mx results 
216 # (LL: likelihood; EC: expected covariance, EM: expected means)
217 # -----------------------------------------------------------------------------