Updated copyright to 2013 for R/ demo/ models/passing and src/ folders, and also...
[openmx:openmx.git] / demo / TwoFactorModel_PathCov.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: TwoFactorModel_PathCov.R  
19 # Author: Ryne Estabrook
20 # Date: 2009.08.01 
21 #
22 # ModelType: Factor
23 # DataType: Continuous
24 # Field: None
25 #
26 # Purpose: 
27 #       Two Factor model to estimate factor loadings, residual variances and means
28 #       Path style model input - Covariance matrix data input
29 #
30 # RevisionHistory:
31 #       Hermine Maes -- 2009.10.08 updated & reformatted
32 #       Ross Gore -- 2011.06.06 added Model, Data & Field metadata
33 # -----------------------------------------------------------------------------
34
35 require(OpenMx)
36 # Load Library
37 # -----------------------------------------------------------------------------
38
39 myFADataCov <- matrix(
40       c(0.997, 0.642, 0.611, 0.672, 0.637, 0.677, 0.342, 0.299, 0.337,
41         0.642, 1.025, 0.608, 0.668, 0.643, 0.676, 0.273, 0.282, 0.287,
42         0.611, 0.608, 0.984, 0.633, 0.657, 0.626, 0.286, 0.287, 0.264,
43         0.672, 0.668, 0.633, 1.003, 0.676, 0.665, 0.330, 0.290, 0.274,
44         0.637, 0.643, 0.657, 0.676, 1.028, 0.654, 0.328, 0.317, 0.331,
45         0.677, 0.676, 0.626, 0.665, 0.654, 1.020, 0.323, 0.341, 0.349,
46         0.342, 0.273, 0.286, 0.330, 0.328, 0.323, 0.993, 0.472, 0.467,
47         0.299, 0.282, 0.287, 0.290, 0.317, 0.341, 0.472, 0.978, 0.507,
48         0.337, 0.287, 0.264, 0.274, 0.331, 0.349, 0.467, 0.507, 1.059),
49       nrow=9,
50       dimnames=list(
51           c("x1", "x2", "x3", "x4", "x5", "x6", "y1", "y2", "y3"),
52           c("x1", "x2", "x3", "x4", "x5", "x6", "y1", "y2", "y3")),
53 )
54
55 twoFactorCov <- myFADataCov[c("x1","x2","x3","y1","y2","y3"),c("x1","x2","x3","y1","y2","y3")]
56   
57 myFADataMeans <- c(2.988, 3.011, 2.986, 3.053, 3.016, 3.010, 2.955, 2.956, 2.967)
58 names(myFADataMeans) <- c("x1", "x2", "x3", "x4", "x5", "x6", "y1", "y2", "y3")
59   
60 twoFactorMeans <- myFADataMeans[c(1:3,7:9)]
61 # Prepare Data
62 # -----------------------------------------------------------------------------
63
64 twoFactorModel <- mxModel("Two Factor Model Path", type="RAM",
65     mxData(
66         observed=twoFactorCov, 
67         type="cov", 
68         numObs=500, 
69         means=twoFactorMeans
70     ),
71     manifestVars=c("x1", "x2", "x3", "y1", "y2", "y3"),
72     latentVars=c("F1","F2"),
73     mxPath(
74         from=c("x1", "x2", "x3", "y1", "y2", "y3"),
75         arrows=2, 
76         free=TRUE, 
77         values=1,
78         labels=c("e1","e2","e3","e4","e5","e6")
79     ),
80         # residual variances
81         # -------------------------------------
82     mxPath(
83         from=c("F1","F2"),
84         arrows=2,
85                 connect="unique.pairs",
86         free=TRUE,
87         values=c(1, .5, 1),
88         labels=c("varF1", "cov", "varF2")
89     ), 
90         # latent variances and covaraince
91         # -------------------------------------
92     mxPath(
93         from="F1",
94         to=c("x1","x2","x3"),
95         arrows=1,
96         free=c(FALSE,TRUE,TRUE),
97         values=c(1,1,1),
98         labels=c("l1","l2","l3")
99     ),
100         # factor loadings for x variables
101         # -------------------------------------
102     mxPath(
103         from="F2",
104         to=c("y1","y2","y3"),
105         arrows=1,
106         free=c(FALSE,TRUE,TRUE),
107         values=c(1,1,1),
108         labels=c("l4","l5","l6")
109     ),
110         # factor loadings for y variables
111         # -------------------------------------
112     mxPath(
113         from="one",
114         to=c("x1","x2","x3","y1","y2","y3","F1","F2"),
115         arrows=1,
116         free=c(T ,T, T, T, T, T, F, F),
117         values=c(1,1,1,1,1,1,0,0),
118         labels=c("meanx1","meanx2","meanx3","meany1","meany2","meany3",NA,NA)
119     )
120         # means
121         # -------------------------------------
122 ) # close model
123 # Create an MxModel object
124 # -----------------------------------------------------------------------------
125       
126 twoFactorFit <- mxRun(twoFactorModel)
127
128 summary(twoFactorFit)
129 twoFactorFit@output$estimate
130
131 omxCheckCloseEnough(twoFactorFit@output$estimate[["l2"]], 0.9720, 0.01)
132 omxCheckCloseEnough(twoFactorFit@output$estimate[["l3"]], 0.9310, 0.01)
133 omxCheckCloseEnough(twoFactorFit@output$estimate[["l5"]], 1.0498, 0.01)
134 omxCheckCloseEnough(twoFactorFit@output$estimate[["l6"]], 1.0533, 0.01)
135 omxCheckCloseEnough(twoFactorFit@output$estimate[["varF1"]], 0.6622, 0.01)
136 omxCheckCloseEnough(twoFactorFit@output$estimate[["varF2"]], 0.4510, 0.01)
137 omxCheckCloseEnough(twoFactorFit@output$estimate[["cov"]], 0.2958, 0.01)
138 omxCheckCloseEnough(twoFactorFit@output$estimate[["e1"]], 0.3348, 0.01)
139 omxCheckCloseEnough(twoFactorFit@output$estimate[["e2"]], 0.3994, 0.01)
140 omxCheckCloseEnough(twoFactorFit@output$estimate[["e3"]], 0.4101, 0.01)
141 omxCheckCloseEnough(twoFactorFit@output$estimate[["e4"]], 0.5420, 0.01)
142 omxCheckCloseEnough(twoFactorFit@output$estimate[["e5"]], 0.4809, 0.01)
143 omxCheckCloseEnough(twoFactorFit@output$estimate[["e6"]], 0.5586, 0.01)
144 omxCheckCloseEnough(twoFactorFit@output$estimate[["meanx1"]], 2.988, 0.01)
145 omxCheckCloseEnough(twoFactorFit@output$estimate[["meanx2"]], 3.011, 0.01)
146 omxCheckCloseEnough(twoFactorFit@output$estimate[["meanx3"]], 2.986, 0.01)
147 omxCheckCloseEnough(twoFactorFit@output$estimate[["meany1"]], 2.955, 0.01)
148 omxCheckCloseEnough(twoFactorFit@output$estimate[["meany2"]], 2.956, 0.01)
149 omxCheckCloseEnough(twoFactorFit@output$estimate[["meany3"]], 2.967, 0.01)
150 # Compare OpenMx results to Mx results 
151 # -----------------------------------------------------------------------------