Add checks for AIC, BIC
[openmx:openmx.git] / models / passing / SummaryCheck.R
1 #
2 #   Copyright 2007-2014 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 # Author: Michael D. Hunter
18 # Date: 2014.06.10
19 # Filename: SummaryCheck.R
20 # Purpose: Actually test that pieces of MxSummary() do what they're supposed 
21 #  to do.
22 # Procedure: If you're going to muck around with some part of summary,
23 #  first make sure this file tests that part.  If this file does not test the
24 #  part of summary that you want to change, then you must add to this test file
25 #  before modifying summary.  Once the test has been added, muck around with
26 #  summary to your heart's content.
27 # TODO add multigroup case to make sure this works
28 # TODO add model with constraint to make sure it adds back DoF
29 # TODO add RMSEA confidence interval
30 #------------------------------------------------------------------------------
31
32
33 #------------------------------------------------------------------------------
34 # Load package and data
35
36 require(OpenMx)
37
38 data(demoTwoFactor)
39
40 #------------------------------------------------------------------------------
41 # Specify raw data model
42
43 vpf <- 5
44 varnames <- names(demoTwoFactor)
45 latnames <- c("F1", "F2")
46
47 fl <- mxPath(from=rep(latnames, each=vpf), to=varnames, values=.8, labels=paste("load_", varnames, sep=""), free=TRUE)
48 fc <- mxPath(latnames, connect="unique.pairs", arrows=2, free=c(FALSE, TRUE, FALSE), values=c(1, .5, 1), labels=c("varF1", "covF1F2", "varF2"))
49 mm <- mxPath(from="one", to=varnames, free=TRUE, values=colMeans(demoTwoFactor), labels=paste("mean_", varnames, sep=""))
50 rv <- mxPath(from=varnames, arrows=2, free=TRUE, values=.2, labels=paste("resid_", varnames, sep=""))
51
52 raw <- mxModel("Raw Test Model to Check MxSummary",
53         type="RAM", manifestVars=varnames, latentVars=latnames,
54         fl, fc, mm, rv,
55         mxData(demoTwoFactor, "raw")
56 )
57
58 raw.fit <- mxRun(raw)
59
60 #------------------------------------------------------------------------------
61 # Specify Saturated raw data model
62
63 sat.fit <- omxSaturatedModel(raw, run=TRUE)
64
65
66 #------------------------------------------------------------------------------
67 # Specify cov data models
68 cov <- mxModel("Covariance Test Model to Check MxSummary",
69         type="RAM", manifestVars=varnames, latentVars=latnames,
70         fl, fc, rv,
71         mxData(cov(demoTwoFactor), "cov", numObs=nrow(demoTwoFactor))
72 )
73
74 cov.fit <- mxRun(cov)
75
76 covm <- mxModel("Covariance and Means Test Model to Check MxSummary",
77         type="RAM", manifestVars=varnames, latentVars=latnames,
78         fl, fc, rv, mm,
79         mxData(cov(demoTwoFactor), "cov", numObs=nrow(demoTwoFactor), means=colMeans(demoTwoFactor))
80 )
81
82 covm.fit <- mxRun(covm)
83
84
85 #------------------------------------------------------------------------------
86 # Get summaries of models
87
88 raw.sum <- summary(raw.fit)
89 sat.sum <- summary(sat.fit)
90 raws.sum <- summary(raw.fit, SaturatedLikelihood=sat.fit)
91 cov.sum <- summary(cov.fit)
92 covs.sum <- summary(covm.fit)
93
94
95 #------------------------------------------------------------------------------
96 # Compare Fit statistics to what they should be.
97
98 # Things to check for raw w/o sat, raw w/ sat, cov, cov & means, multigroup
99 #       number of observed statistics
100 #       TODO add multigroup case to make sure this works
101 #       degrees of freedom
102 #       -2LL, sat -2LL, numObs
103 #       chi-square, chi dof, chi p
104 #       AIC df AIV, param, BIC df, BIC param, BIC sample size
105 #       CFI, TLI, RMSEA
106 #       TODO add RMSEA confidence interval
107
108 #       number of observed statistics
109 omxCheckEquals(raw.sum$observedStatistics, 5000)
110 omxCheckEquals(raws.sum$observedStatistics, 5000)
111 omxCheckEquals(cov.sum$observedStatistics, 55)
112 omxCheckEquals(covs.sum$observedStatistics, 65)
113
114 #       TODO add multigroup case to make sure this works
115
116
117 #       degrees of freedom
118 # TODO add model with constraint to make sure it adds back DoF
119 ep <- 31
120 omxCheckEquals(raw.sum$degreesOfFreedom, 5000-ep)
121 omxCheckEquals(raws.sum$degreesOfFreedom, 5000-ep)
122 omxCheckEquals(sat.sum$degreesOfFreedom, 5000-sat.sum$estimatedParameters)
123 omxCheckEquals(raws.sum$satDoF, 5000-sat.sum$estimatedParameters)
124 omxCheckEquals(cov.sum$degreesOfFreedom, 55-(ep-10)) #no means estimated
125 omxCheckEquals(covs.sum$degreesOfFreedom, 65-ep)
126
127
128 #       -2LL, sat -2LL, numObs
129 omxCheckWithinPercentError(raw.sum$Minus2LogLikelihood, 9236.675, percent=1e-4)
130 omxCheckTrue(is.na(raw.sum$SaturatedLikelihood))
131 omxCheckWithinPercentError(raws.sum$SaturatedLikelihood, 9186.911, percent=1e-4)
132 omxCheckEquals(raw.sum$numObs, 500)
133 omxCheckEquals(raws.sum$numObs, 500)
134 omxCheckEquals(cov.sum$numObs, 500)
135 omxCheckEquals(covs.sum$numObs, 500)
136
137
138 #       chi-square, chi dof, chi p
139 omxCheckTrue(is.na(raw.sum$Chi))
140 omxCheckTrue(!is.na(raws.sum$Chi))
141 #cov.sum$Chi
142 #covs.sum$Chi
143
144 # ChiDoF = df - sat.df, i.e. df = obsStat-ep, sat.df = obsStat-sat.ep, so ChiDoF = sat.ep - ep
145 chi.df <- 34
146 #omxCheckTrue(is.na(raw.sum$ChiDoF)) # is still populated, but it reported as NA when satLik is NA
147 omxCheckEquals(raws.sum$ChiDoF, chi.df)
148 omxCheckEquals(cov.sum$ChiDoF, chi.df)
149 omxCheckEquals(covs.sum$ChiDoF, chi.df)
150
151 #raw.sum$p
152 #raws.sum$p
153 #cov.sum$p
154 #covs.sum$p
155
156
157 #       AIC df, AIC param, BIC df, BIC param, BIC sample size
158 #       CFI, TLI, RMSEA
159 #       TODO add RMSEA confidence interval
160
161 omxCheckTrue(is.na(raw.sum$RMSEA))
162 omxCheckTrue(!is.na(raws.sum$RMSEA))
163
164 omxCheckCloseEnough(raw.sum$informationCriteria['AIC:','par'], 9298.67, .01)
165 omxCheckCloseEnough(raw.sum$informationCriteria['BIC:','par'], 9429.32, .01)
166
167 omxCheckCloseEnough(sat.sum$informationCriteria['AIC:','par'], 9316.91, .01)
168 omxCheckCloseEnough(sat.sum$informationCriteria['BIC:','par'], 9590.86, .01)
169
170 omxCheckCloseEnough(raws.sum$informationCriteria['AIC:','par'], 9298.67, .01)
171 omxCheckCloseEnough(raws.sum$informationCriteria['BIC:','par'], 9429.32, .01)
172
173 omxCheckCloseEnough(cov.sum$informationCriteria['AIC:','par'], 91.66, .01)
174 omxCheckCloseEnough(cov.sum$informationCriteria['BIC:','par'], 180.17, .01)
175
176 omxCheckCloseEnough(covs.sum$informationCriteria['AIC:','par'], 111.66, .01)
177 omxCheckCloseEnough(covs.sum$informationCriteria['BIC:','par'], 242.31, .01)