Update demos directory to reflect new style header and new commenting style.
[openmx:openmx.git] / demo / OneFactorModel_PathCov.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: OneFactorModel_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 #      One Factor model to estimate factor loadings, 
28 #      residual variances and means
29 #      Path style model input - Covariance matrix data input
30 #
31 # RevisionHistory:
32 #      Hermine Maes -- 2009.10.08 updated & reformatted
33 #      Ross Gore -- 2011.06.06 added Model, Data & Field metadata
34 # -----------------------------------------------------------------------------
35
36 require(OpenMx)
37 # Load Library
38 # -----------------------------------------------------------------------------
39
40 myFADataCov<-matrix(
41         c(0.997, 0.642, 0.611, 0.672, 0.637, 0.677,
42           0.642, 1.025, 0.608, 0.668, 0.643, 0.676,
43           0.611, 0.608, 0.984, 0.633, 0.657, 0.626,
44           0.672, 0.668, 0.633, 1.003, 0.676, 0.665,
45           0.637, 0.643, 0.657, 0.676, 1.028, 0.654,
46           0.677, 0.676, 0.626, 0.665, 0.654, 1.020),
47         nrow=6,
48         dimnames=list(
49                 c("x1","x2","x3","x4","x5","x6"),
50                 c("x1","x2","x3","x4","x5","x6"))
51 )
52
53 myFADataMeans <- c(2.988, 3.011, 2.986, 3.053, 3.016, 3.010)
54 names(myFADataMeans) <- c("x1","x2","x3","x4","x5","x6")
55 # Prepare Data
56 # -----------------------------------------------------------------------------
57
58 oneFactorModel <- mxModel("Common Factor Model Path Specification", 
59         type="RAM",
60         mxData(
61                 observed=myFADataCov, 
62                 type="cov", 
63                 numObs=500,
64                 mean=myFADataMeans
65         ),
66         manifestVars=c("x1","x2","x3","x4","x5","x6"),
67         latentVars="F1",
68         mxPath(
69                 from=c("x1","x2","x3","x4","x5","x6"),
70                 arrows=2,
71                 free=TRUE,
72                 values=c(1,1,1,1,1,1),
73                 labels=c("e1","e2","e3","e4","e5","e6")
74         ),
75         # residual variances
76         # -------------------------------------
77         mxPath(from="F1",
78                 arrows=2,
79                 free=TRUE,
80                 values=1,
81                 labels ="varF1"
82         ),
83         # -------------------------------------
84         # latent variance
85         mxPath(from="F1",
86                 to=c("x1","x2","x3","x4","x5","x6"),
87                 arrows=1,
88                 free=c(FALSE,TRUE,TRUE,TRUE,TRUE,TRUE),
89                 values=c(1,1,1,1,1,1),
90                 labels =c("l1","l2","l3","l4","l5","l6")
91         ),
92         # factor loadings
93         # --------------------------------------
94         mxPath(from="one",
95                 to=c("x1","x2","x3","x4","x5","x6","F1"),
96                 arrows=1,
97                 free=c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),
98                 values=c(1,1,1,1,1,1,0),
99                 labels =c("meanx1","meanx2","meanx3","meanx4","meanx5","meanx6",NA)
100         )
101         # means
102         # ------------------------------------- 
103 ) # close model
104 # Create an MxModel object
105 # -----------------------------------------------------------------------------
106
107 oneFactorFit <- mxRun(oneFactorModel)
108
109 summary(oneFactorFit)
110 oneFactorFit@output$estimate
111
112 omxCheckCloseEnough(oneFactorFit@output$estimate[["l2"]], 0.999, 0.01)
113 omxCheckCloseEnough(oneFactorFit@output$estimate[["l3"]], 0.959, 0.01)
114 omxCheckCloseEnough(oneFactorFit@output$estimate[["l4"]], 1.028, 0.01)
115 omxCheckCloseEnough(oneFactorFit@output$estimate[["l5"]], 1.008, 0.01)
116 omxCheckCloseEnough(oneFactorFit@output$estimate[["l6"]], 1.021, 0.01)
117 omxCheckCloseEnough(oneFactorFit@output$estimate[["varF1"]], 0.645, 0.01)
118 omxCheckCloseEnough(oneFactorFit@output$estimate[["e1"]], 0.350, 0.01)
119 omxCheckCloseEnough(oneFactorFit@output$estimate[["e2"]], 0.379, 0.01)
120 omxCheckCloseEnough(oneFactorFit@output$estimate[["e3"]], 0.389, 0.01)
121 omxCheckCloseEnough(oneFactorFit@output$estimate[["e4"]], 0.320, 0.01)
122 omxCheckCloseEnough(oneFactorFit@output$estimate[["e5"]], 0.370, 0.01)
123 omxCheckCloseEnough(oneFactorFit@output$estimate[["e6"]], 0.346, 0.01)
124 omxCheckCloseEnough(oneFactorFit@output$estimate[["meanx1"]], 2.988, 0.01)
125 omxCheckCloseEnough(oneFactorFit@output$estimate[["meanx2"]], 3.011, 0.01)
126 omxCheckCloseEnough(oneFactorFit@output$estimate[["meanx3"]], 2.986, 0.01)
127 omxCheckCloseEnough(oneFactorFit@output$estimate[["meanx4"]], 3.053, 0.01)
128 omxCheckCloseEnough(oneFactorFit@output$estimate[["meanx5"]], 3.016, 0.01)
129 omxCheckCloseEnough(oneFactorFit@output$estimate[["meanx6"]], 3.010, 0.01)
130 # Compare OpenMx results to Mx results 
131 # -----------------------------------------------------------------------------