Renamed Row Objective FIML demo and added to Users Guide
[openmx:openmx.git] / demo / RowObjectiveFIMLBivariateSaturated.R
1 #
2 #   Copyright 2007-2010 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: FIMLRowObjectiveBivariateCorrelation.R
19 # Author: Hermine Maes 
20 # Date: 2009.08.01 
21 #
22 # Purpose: 
23 #      Optimization Example in OpenMx: Testing significance of correlation
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 mxRowObjective
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 testVars <- c('X','Y')
50 names(testData) <- testVars
51 summary(testData)
52 cov(testData)
53 # Simulate Data: two standardized variables X & Y with correlation of .5
54
55
56 # -----------------------------------------------------------------------------
57
58 bivCorModelSpec <- mxModel(
59     name="FIML Saturated Bivariate",
60     mxData(
61         observed=testData, 
62         type="raw",
63     ),
64     mxMatrix(
65         type="Full", 
66         nrow=1, 
67         ncol=2, 
68         free=TRUE, 
69         values=c(0,0), 
70         name="expMean"
71     ), 
72     mxMatrix(
73         type="Symm",
74         nrow=2, 
75         ncol=2,
76         values=c(.21, .2, .2, .21),
77         free=TRUE,
78         name='expCov'
79     )
80 )
81
82 bivCorFiltering <- mxModel(
83     model=bivCorModelSpec,
84     mxAlgebra(
85         expression=omxSelectRowsAndCols(expCov, existenceVector),
86         name="filteredExpCov",
87     ),
88     mxAlgebra(
89         expression=omxSelectCols(expMean, existenceVector),
90         name="filteredExpMean",
91     ),
92     mxAlgebra(
93         expression=sum(existenceVector),
94         name="numVar_i")
95 )
96
97 bivCorCalc <- mxModel(
98     model=bivCorFiltering,
99     mxAlgebra(
100         expression = log(2*pi),
101         name = "log2pi"
102     ),
103     mxAlgebra(
104         expression=log2pi %*% numVar_i + log(det(filteredExpCov)),
105         name ="firstHalfCalc",
106     ),
107     mxAlgebra(
108         expression=(filteredDataRow - filteredExpMean) %&% solve(filteredExpCov),
109         name = "secondHalfCalc",
110     )
111 )
112
113 bivCorRowObj <- mxModel(
114     model=bivCorCalc,
115     mxAlgebra(
116         expression=(firstHalfCalc + secondHalfCalc),
117         name="rowAlgebra",
118     ),
119     mxAlgebra(
120         expression=sum(rowResults),
121         name = "reduceAlgebra",
122     ),
123     mxRowObjective(
124         rowAlgebra='rowAlgebra',
125         reduceAlgebra='reduceAlgebra',
126         dimnames=c('X','Y'),
127     )
128 )
129
130 bivCorTotal <- bivCorRowObj
131
132 # Fit Saturated Model with Raw Data and Matrix-style Input.
133 # Estimate with Full Information Maximum Likelihood (FIML).
134 # FIML is implemented here as an mxRowObjective function for pedagogical reasons only.
135 # FIML for one row of data is
136 #   2*log(2*pi) + log(det(Cov)) + (Row - Mean) %*% solve(Cov) %*% t(Row - Mean)
137 #  where Cov is the filtered expected covariance matrix
138 #        Row is the filtered data row
139 #        Mean is the filtered expected means row vector
140 #        solve(*) is the inverse of *
141 #        t(*) is the transpose of *
142 #        det(*) is the determinant of *
143 # FIML for a whole data set is the sum of all the FIML rows.
144 # -----------------------------------------------------------------------------
145
146 bivCorFit <- mxRun(bivCorTotal)
147 EM <- mxEval(expMean, bivCorFit)
148 EC <- mxEval(expCov, bivCorFit)
149 LL <- mxEval(objective, bivCorFit)
150 # Run Model and Generate Output
151 # -----------------------------------------------------------------------------
152
153 Mx.EM <- matrix(c(0.03211656, -0.004883885), nrow=1, ncol=2)
154 Mx.EC <- matrix(c(1.0092853, 0.4813504, 0.4813504, 0.9935390), nrow=2, ncol=2)
155 Mx.LL <- 5415.772
156 # Mx Answers of Saturated Model Hard-coded
157 # -----------------------------------------------------------------------------
158
159 omxCheckCloseEnough(LL, Mx.LL, .001)
160 omxCheckCloseEnough(EC, Mx.EC, .001)
161 omxCheckCloseEnough(EM, Mx.EM, .001)
162 # Compare OpenMx Results to Mx Results 
163 # LL: likelihood; EC: expected covariance, EM: expected means
164 # -----------------------------------------------------------------------------