Allow multiple 'references' in single call to omxInterval function.
[openmx:openmx.git] / models / failing / SimpleConfidenceIntervals.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 require(OpenMx)
17
18 #Prepare Data
19 # -----------------------------------------------------------------------
20 data(twinData)
21 summary(twinData)
22 selVars <- c('bmi1','bmi2')
23 mzfData <- as.matrix(subset(twinData, zyg==1, c(bmi1,bmi2)))
24 dzfData <- as.matrix(subset(twinData, zyg==3, c(bmi1,bmi2)))
25 colMeans(mzfData,na.rm=TRUE)
26 colMeans(dzfData,na.rm=TRUE)
27 cov(mzfData,use="complete")
28 cov(dzfData,use="complete")
29
30 #Fit ACE Model with RawData and Matrices Input
31 # -----------------------------------------------------------------------
32 twinACEModel <- mxModel("twinACE",
33         mxModel("common",
34                 # Matrices X, Y, and Z to store a, c, and e path coefficients
35                 mxMatrix(
36                         type="Full", 
37                         nrow=1, 
38                         ncol=1, 
39                         free=TRUE,  
40                         values=.6,  
41                         label="a", 
42                         name="X"
43                 ), 
44                 mxMatrix(
45                         type="Full", 
46                         nrow=1, 
47                         ncol=1, 
48                         free=TRUE,  
49                         values=.6,  
50                         label="c", 
51                         name="Y"
52                 ),
53                 mxMatrix(
54                         type="Full", 
55                         nrow=1, 
56                         ncol=1, 
57                         free=TRUE,  
58                         values=.6,  
59                         label="e", 
60                         name="Z"
61                 ),
62                 # Matrices A, C, and E compute variance components
63                 mxAlgebra(
64                         expression=X %*% t(X), 
65                         name="A"
66                 ), 
67                 mxAlgebra(
68                         expression=Y %*% t(Y), 
69                         name="C"
70                 ), 
71                 mxAlgebra(
72                         expression=Z %*% t(Z), 
73                         name="E"
74                 ),
75                 mxMatrix(
76                         type="Full", 
77                         nrow=1, 
78                         ncol=2, 
79                         free=TRUE, 
80                         values= 20,
81                         label="mean", 
82                         name="expMean"
83                 ),
84             # Algebra for expected variance/covariance matrix in MZ
85             mxAlgebra(
86                         expression= rbind  (cbind(A+C+E , A+C),
87                                                                 cbind(A+C   , A+C+E)), 
88                         name="expCovMZ"
89                 ),
90             # Algebra for expected variance/covariance matrix in DZ
91             # note use of 0.5, converted to 1*1 matrix
92             mxAlgebra(
93                         expression= rbind  (cbind(A+C+E     , 0.5%x%A+C),
94                                                                 cbind(0.5%x%A+C , A+C+E)), 
95                         name="expCovDZ"
96                 )
97         ),
98         mxModel("MZ",
99             mxData(
100                         observed=mzfData, 
101                         type="raw"
102                 ), 
103             mxFIMLObjective(
104                         covariance="common.expCovMZ", 
105                         means="common.expMean", 
106                         dimnames=selVars
107                 )
108         ),
109         mxModel("DZ",
110             mxData(
111                         observed=dzfData, 
112                         type="raw"
113                 ), 
114             mxFIMLObjective(
115                         covariance="common.expCovDZ", 
116                         means="common.expMean", 
117                         dimnames=selVars
118                 )
119         ),
120         mxAlgebra(
121                 expression=MZ.objective + DZ.objective, 
122                 name="twin"
123         ), 
124         mxAlgebraObjective("twin"),
125     omxInterval(c("common.A[1,1]", "common.C[1,1]", "common.E[1,1]"), 3.84, 3.84)
126 )       
127 #Run ACE model
128 # -----------------------------------------------------------------------
129 twinACEFit <- mxRun(twinACEModel, intervals=TRUE)
130
131 CIaupper <- mxModel(twinACEFit, name = 'A_CIupper',
132                 mxMatrix("Full", values=mxEval(objective, twinACEFit), name="oldfit"), 
133                 mxAlgebra(((oldfit + 3.84) - (MZ.objective + DZ.objective))^2 - common.A,name="upperCIa"), 
134                 mxAlgebraObjective("upperCIa"))
135
136 CIalower <- mxModel(twinACEFit, name = 'A_CIlower',
137                 mxMatrix("Full", values=mxEval(objective, twinACEFit), name="oldfit"),
138                 mxAlgebra(((oldfit + 3.84) - (MZ.objective + DZ.objective))^2 + common.A,name="lowerCIa"),
139                 mxAlgebraObjective("lowerCIa"))
140
141 runCIalower <- suppressWarnings(mxRun(CIalower, intervals=FALSE))
142 runCIaupper <- suppressWarnings(mxRun(CIaupper, intervals=FALSE))
143
144 CIcupper <- mxModel(twinACEFit, name = 'C_CIupper',
145                 mxMatrix("Full", values=mxEval(objective, twinACEFit), name="oldfit"), 
146                 mxAlgebra(((oldfit + 3.84) - (MZ.objective + DZ.objective))^2 - common.C,name="upperCIc"), 
147                 mxAlgebraObjective("upperCIc"))
148
149 CIclower <- mxModel(twinACEFit, name = 'C_CIlower',
150                 mxMatrix("Full", values=mxEval(objective, twinACEFit), name="oldfit"),
151                 mxAlgebra(((oldfit + 3.84) - (MZ.objective + DZ.objective))^2 + common.C,name="lowerCIc"),
152                 mxAlgebraObjective("lowerCIc"))
153
154 runCIclower <- suppressWarnings(mxRun(CIclower, intervals=FALSE))
155 runCIcupper <- suppressWarnings(mxRun(CIcupper, intervals=FALSE))
156
157 CIeupper <- mxModel(twinACEFit, name = 'E_CIupper',
158                 mxMatrix("Full", values=mxEval(objective, twinACEFit), name="oldfit"), 
159                 mxAlgebra(((oldfit + 3.84) - (MZ.objective + DZ.objective))^2 - common.E,name="upperCIe"), 
160                 mxAlgebraObjective("upperCIe"))
161
162 CIelower <- mxModel(twinACEFit, name = 'E_CIlower',
163                 mxMatrix("Full", values=mxEval(objective, twinACEFit), name="oldfit"),
164                 mxAlgebra(((oldfit + 3.84) - (MZ.objective + DZ.objective))^2 + common.E,name="lowerCIe"),
165                 mxAlgebraObjective("lowerCIe"))
166
167 runCIelower <- suppressWarnings(mxRun(CIelower, intervals=FALSE))
168 runCIeupper <- suppressWarnings(mxRun(CIeupper, intervals=FALSE))
169
170
171 omxCheckCloseEnough(twinACEFit@output$confidenceIntervalCodes[1,1], runCIalower@output$status[1], .1)
172 omxCheckCloseEnough(twinACEFit@output$confidenceIntervalCodes[1,2], runCIaupper@output$status[1], .1)
173 omxCheckCloseEnough(twinACEFit@output$confidenceIntervals[1, 1], mxEval(common.A, runCIalower), .001)
174 omxCheckCloseEnough(twinACEFit@output$confidenceIntervals[1, 2], mxEval(common.A, runCIaupper), .001)
175
176 omxCheckCloseEnough(twinACEFit@output$confidenceIntervals[2, 1], mxEval(common.C, runCIclower), .001)
177 omxCheckCloseEnough(twinACEFit@output$confidenceIntervals[2, 2], mxEval(common.C, runCIcupper), .001)
178 omxCheckCloseEnough(twinACEFit@output$confidenceIntervalCodes[2,1], runCIclower@output$status[1], .1)
179 omxCheckCloseEnough(twinACEFit@output$confidenceIntervalCodes[2,2], runCIcupper@output$status[1], .1)
180
181
182 omxCheckCloseEnough(twinACEFit@output$confidenceIntervalCodes[3,1], runCIelower@output$status[1], .1)
183 omxCheckCloseEnough(twinACEFit@output$confidenceIntervalCodes[3,2], runCIeupper@output$status[1], .1)
184 omxCheckCloseEnough(twinACEFit@output$confidenceIntervals[3, 1], mxEval(common.E, runCIelower), .001)
185 omxCheckCloseEnough(twinACEFit@output$confidenceIntervals[3, 2], mxEval(common.E, runCIeupper), .001)
186
187 twinACEFit@output$confidenceIntervalCodes
188
189