(Hopefully) fixed the error with the numerous bivariate saturated models. I believe...
[openmx:openmx.git] / demo / RowObjectiveFIMLBivariateSaturated.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: RowObjectiveFIMLBivariateSaturated.R
19 # Author: Hermine Maes 
20 # Date: 2009.08.01 
21 #
22 # Purpose: 
23 #      Demonstrate an mxFitFunctionRow function implementation of FIML
24 #
25 # ModelType: Saturated
26 # DataType: Simulated Continuous
27 # Field: None 
28 #
29 # RevisionHistory:
30 #       HermineMaes -- 2009.10.08 updated & reformatted
31 #       RossGore -- 2011.04.10 modified to implement FIML via mxFitFunctionRow
32 #       MikeHunter -- 2011.04.11 debugged Gore implementation above
33 #       MikeHunter -- 2011.05.03 Renamed from BivariateCorrelation.R to 
34 #                                        FIMLRowObjectiveBivariateCorrelation.R
35 #       MikeHunter -- 2011.05.05 modified to use omxSelect* functions
36 #       MikeHunter -- 2011.05.26 Adjusted spacing & comments for readability.
37 #   MikeHunter -- 2011.06.30 Adjusted model structure to match docs file in User's Guide and corrected formula error.
38 # -----------------------------------------------------------------------------
39
40 require(OpenMx)
41 require(MASS)
42 # Load Library
43 # -----------------------------------------------------------------------------
44
45 set.seed(200)
46 rs <- .5
47 xy <- mvrnorm (1000, c(0,0), matrix(c(1, rs, rs, 1), nrow=2, ncol=2))
48 testData <- as.data.frame(xy)
49 testData <- testData[, order(apply(testData, 2, var))[2:1]] #put the data columns in order from largest to smallest variance
50 testVars <- c('X','Y')
51 names(testData) <- testVars
52 summary(testData)
53 cov(testData)
54 # Simulate Data: two standardized variables X & Y with correlation of .5
55
56
57 # -----------------------------------------------------------------------------
58
59 bivCorModelSpec <- mxModel(
60     name="FIML Saturated Bivariate",
61     mxData(
62         observed=testData, 
63         type="raw",
64     ),
65     mxMatrix(
66         type="Full", 
67         nrow=1, 
68         ncol=2, 
69         free=TRUE, 
70         values=c(0,0), 
71         name="expMean"
72     ), 
73     mxMatrix(
74         type="Symm",
75         nrow=2, 
76         ncol=2,
77         values=c(.21, .2, .2, .21),
78         free=TRUE,
79         name='expCov'
80     )
81 )
82
83 bivCorFiltering <- mxModel(
84     model=bivCorModelSpec,
85     mxAlgebra(
86         expression=omxSelectRowsAndCols(expCov, existenceVector),
87         name="filteredExpCov",
88     ),
89     mxAlgebra(
90         expression=omxSelectCols(expMean, existenceVector),
91         name="filteredExpMean",
92     ),
93     mxAlgebra(
94         expression=sum(existenceVector),
95         name="numVar_i")
96 )
97
98 bivCorCalc <- mxModel(
99     model=bivCorFiltering,
100     mxAlgebra(
101         expression = log(2*pi),
102         name = "log2pi"
103     ),
104     mxAlgebra(
105         expression=log2pi %*% numVar_i + log(det(filteredExpCov)),
106         name ="firstHalfCalc",
107     ),
108     mxAlgebra(
109         expression=(filteredDataRow - filteredExpMean) %&% solve(filteredExpCov),
110         name = "secondHalfCalc",
111     )
112 )
113
114 bivCorRowObj <- mxModel(
115     model=bivCorCalc,
116     mxAlgebra(
117         expression=(firstHalfCalc + secondHalfCalc),
118         name="rowAlgebra",
119     ),
120     mxAlgebra(
121         expression=sum(rowResults),
122         name = "reduceAlgebra",
123     ),
124     mxFitFunctionRow(
125         rowAlgebra='rowAlgebra',
126         reduceAlgebra='reduceAlgebra',
127         dimnames=c('X','Y'),
128     )
129 )
130
131 bivCorTotal <- bivCorRowObj
132
133 # Fit Saturated Model with Raw Data and Matrix-style Input.
134 # Estimate with Full Information Maximum Likelihood (FIML).
135 # FIML is implemented here as an mxFitFunctionRow function for pedagogical reasons only.
136 # FIML for one row of data is
137 #   2*log(2*pi) + log(det(Cov)) + (Row - Mean) %*% solve(Cov) %*% t(Row - Mean)
138 #  where Cov is the filtered expected covariance matrix
139 #        Row is the filtered data row
140 #        Mean is the filtered expected means row vector
141 #        solve(*) is the inverse of *
142 #        t(*) is the transpose of *
143 #        det(*) is the determinant of *
144 # FIML for a whole data set is the sum of all the FIML rows.
145 # -----------------------------------------------------------------------------
146
147 bivCorFit <- mxRun(bivCorTotal)
148 EM <- mxEval(expMean, bivCorFit)
149 EC <- mxEval(expCov, bivCorFit)
150 LL <- mxEval(objective, bivCorFit)
151 # Run Model and Generate Output
152 # -----------------------------------------------------------------------------
153
154 Mx.EM <- matrix(c(0.03211656, -0.004883885), nrow=1, ncol=2)
155 Mx.EC <- matrix(c(1.0092853, 0.4813504, 0.4813504, 0.9935390), nrow=2, ncol=2)
156 Mx.LL <- 5415.772
157 # Mx Answers of Saturated Model Hard-coded
158 # -----------------------------------------------------------------------------
159
160 omxCheckCloseEnough(LL, Mx.LL, .001)
161 omxCheckCloseEnough(EC, Mx.EC, .001)
162 omxCheckCloseEnough(EM, Mx.EM, .001)
163 # Compare OpenMx Results to Mx Results 
164 # LL: likelihood; EC: expected covariance, EM: expected means
165 # -----------------------------------------------------------------------------