(Hopefully) fixed the error with the numerous bivariate saturated models. I believe...
[openmx:openmx.git] / demo / BivariateCorrelation.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: BivariateCorrelation.R  
19 # Author: Hermine Maes
20 # Date: 2009.08.01 
21 #
22 # ModelType: Saturated
23 # DataType: Continuous
24 # Field: None
25 #
26 # Purpose: 
27 #      Optimization Example in OpenMx: Testing significance of correlation
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 Library
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: two standardized variables X & Y with correlation of .5
49 # -----------------------------------------------------------------------------
50
51
52 bivCorModel <- mxModel("bivCor",
53     mxMatrix(
54         type="Full", 
55         nrow=1, 
56         ncol=2, 
57         free=TRUE, 
58         values=c(0,0), 
59         name="expMean"
60     ), 
61     mxMatrix(
62         type="Lower", 
63         nrow=2, 
64         ncol=2, 
65         free=TRUE,
66         values=.5, 
67         name="Chol"
68     ), 
69     mxAlgebra(
70         expression=Chol %*% t(Chol), 
71         name="expCov", 
72     ), 
73     mxData(
74         observed=testData, 
75         type="raw"
76     ), 
77     mxExpectationNormal(
78         covariance="expCov", 
79         means="expMean",
80         dimnames=selVars),
81     mxFitFunctionML()
82     )
83 # Fit Saturated Model with Raw Data and Matrix-style Input
84 # -----------------------------------------------------------------------------
85
86
87 bivCorFit <- mxRun(bivCorModel)
88 EM <- mxEval(expMean, bivCorFit)
89 EC <- mxEval(expCov, bivCorFit)
90 LL <- mxEval(fitfunction, bivCorFit)
91 # Run Model and Generate Output
92 # -----------------------------------------------------------------------------
93
94
95 bivCorModelSub <-mxModel(bivCorModel,
96     mxMatrix(
97         type="Diag", 
98         nrow=2, 
99         ncol=2, 
100         free=TRUE, 
101         name="Chol"
102     )
103 )
104 # Specify SubModel testing Covariance=Zero
105 # -----------------------------------------------------------------------------
106
107
108 bivCorFitSub <- mxRun(bivCorModelSub)
109 EMs <- mxEval(expMean, bivCorFitSub)
110 ECs <- mxEval(expCov, bivCorFitSub)
111 LLs <- mxEval(fitfunction, bivCorFitSub)
112 Chi= LLs-LL;
113 LRT= rbind(LL,LLs,Chi); LRT
114 # Run Model and Generate Output
115 # -----------------------------------------------------------------------------
116
117
118 Mx.EM <- matrix(c(0.03211656, -0.004883885), 1, 2)
119 Mx.EC <- matrix(c(1.0092853, 0.4813504, 0.4813504, 0.9935390), 2, 2)
120 Mx.LL <- 5415.772
121 # Mx Answers of Saturated Model Hard-coded
122 # -----------------------------------------------------------------------------
123
124 omxCheckCloseEnough(LL,Mx.LL,.001)
125 omxCheckCloseEnough(EC,Mx.EC,.001)
126 omxCheckCloseEnough(EM,Mx.EM,.001)
127 # Compare OpenMx Results to Mx Results 
128 # LL: likelihood; EC: expected covariance, EM: expected means
129 # -----------------------------------------------------------------------------