Compute condition number of information matrix
[openmx:openmx.git] / models / nightly / OneFactorOrdinal01_MatrixRaw.R
1 #
2 #   Copyright 2007-2012 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 # Program: OneFactorOrdinal01_MatrixRaw.R  
18 # Author: Michael Neale
19 # Date: 2010.08.14 
20 #
21 # ModelType: Factor
22 # DataType: Ordinal
23 # Field: None
24 #
25 # Purpose: 
26 #      One Factor model to estimate factor loadings, residual variances and means
27 #      Matrix style model input - Raw data input - Ordinal data
28 #
29 # RevisionHistory:
30 #      Hermine Maes -- 2010.09.07 updated & reformatted & respecified
31 #      Ross Gore -- 2011.06.06  added Model, Data & Field metadata
32 # -----------------------------------------------------------------------------
33
34 require(OpenMx)
35 require(MASS)
36 # Load Libraries
37 # -----------------------------------------------------------------------------
38
39 nVariables<-5
40 nFactors<-1
41 nThresholds<-3
42 nSubjects<-500
43 isIdentified<-function(nVariables,nFactors) as.logical(1+sign((nVariables*(nVariables-1)/2) - nVariables*nFactors + nFactors*(nFactors-1)/2))
44
45 isIdentified(nVariables,nFactors) # if this function returns FALSE then model is not identified, otherwise it is.
46 # Set up simulation parameters: 
47 # nVariables>=3, nThresholds>=1, 
48 # nSubjects>=nVariables*nThresholds 
49 # and model should be identified
50 # -------------------------------------
51
52 loadings <- matrix(.7,nrow=nVariables,ncol=nFactors)
53 residuals <- 1 - (loadings * loadings)
54 sigma <- loadings %*% t(loadings) + vec2diag(residuals)
55 mu <- matrix(0,nrow=nVariables,ncol=1)
56
57 set.seed(1234)
58 continuousData <- mvrnorm(n=nSubjects,mu,sigma)
59 # Simulate multivariate normal data
60 # -------------------------------------
61
62 quants<-quantile(continuousData[,1],  probs = c((1:nThresholds)/(nThresholds+1)))
63 ordinalData<-matrix(0,nrow=nSubjects,ncol=nVariables)
64 for(i in 1:nVariables)
65 {
66 ordinalData[,i] <- cut(as.vector(continuousData[,i]),c(-Inf,quants,Inf))
67 }
68 # Chop continuous variables into 
69 # ordinal data with nThresholds+1 
70 # approximately equal categories, 
71 # based on 1st variable
72 # -------------------------------------
73
74 ordinalData <- mxFactor(as.data.frame(ordinalData),levels=c(1:(nThresholds+1)))
75 # Make the ordinal variables into 
76 # R factors
77 # -------------------------------------
78
79 fruitynames<-paste("banana",1:nVariables,sep="")
80 names(ordinalData)<-fruitynames
81 # Name the variables
82 # -------------------------------------
83
84 # Simulate Data
85 # -----------------------------------------------------------------------------
86
87 oneFactorThresholdModel01 <- mxModel("oneFactorThresholdModel01",
88     mxMatrix(
89         type="Full", 
90         nrow=nVariables, 
91         ncol=nFactors, 
92         free=TRUE, 
93         values=0.2, 
94         lbound=-.99, 
95         ubound=2, 
96         name="facLoadings"
97     ),
98     mxMatrix(
99         type="Diag", 
100         nrow=nVariables, 
101         ncol=nVariables,
102         free=TRUE,
103         values=0.9,
104         name="resVariances"
105     ),
106     mxAlgebra(
107         expression=facLoadings %*% t(facLoadings) + resVariances, 
108         name="expCovariances"
109     ),
110     mxMatrix(
111         type="Full", 
112         nrow=1, 
113         ncol=nVariables,
114         free=T,
115         name="expMeans"
116     ),
117     mxMatrix(
118         type="Full", 
119         nrow=nThresholds, 
120         ncol=nVariables,
121         free=rep( c(F,F,rep(T,(nThresholds-2))), nVariables), 
122         values=rep( c(0,1,rep(.2,(nThresholds-2))), nVariables),
123         lbound=rep( c(-Inf,rep(.01,(nThresholds-1))), nVariables),
124         dimnames=list(c(), fruitynames),
125         name="thresholdDeviations"
126     ),
127     mxMatrix(
128         type="Lower",
129         nrow=nThresholds,
130         ncol=nThresholds,
131         free=FALSE,
132         values=1,
133         name="unitLower"
134     ),
135     mxAlgebra(
136         expression=unitLower %*% thresholdDeviations, 
137         name="expThresholds"
138     ),
139         mxMatrix(
140          type="Unit",
141          nrow=nThresholds,
142          ncol=1,
143          name="columnofOnes"
144     ),
145     mxAlgebra(
146          expression=expMeans %x% columnofOnes,
147          name="meansMatrix"
148     ),
149     mxAlgebra(
150          expression=sqrt(t(diag2vec(expCovariances))) %x% columnofOnes,
151          name="variancesMatrix"
152     ),
153     mxAlgebra(
154          expression=(expThresholds - meansMatrix) / variancesMatrix,
155          name="thresholdMatrix"
156     ),
157     mxMatrix( 
158         type="Iden", 
159         nrow=nVariables, 
160         ncol=nVariables, 
161         name="Identity"
162     ),
163     mxAlgebra(
164          expression=solve(sqrt(Identity * expCovariances)) %*% facLoadings,
165          name="facLoadingsMatrix"
166     ),
167     mxData(
168         observed=ordinalData, 
169         type='raw'
170     ),
171     mxFitFunctionML(),mxExpectationNormal(
172         covariance="expCovariances", 
173         means="expMeans", 
174         dimnames=fruitynames, 
175         thresholds="expThresholds"
176     )
177 )
178 # Create Factor Model with Raw Ordinal Data and Matrices Input
179 # -----------------------------------------------------------------------------
180
181 oneFactorThresholdFit01 <- mxRun(oneFactorThresholdModel01, suppressWarnings=TRUE)
182 # Fit the model with mxRun
183 # -----------------------------------------------------------------------------
184
185 summary(oneFactorThresholdFit01)
186 # Print a summary of the results
187 # -----------------------------------------------------------------------------