Update Users Guide pages A through F for the expectation-fit change.
[openmx:openmx.git] / docs / source / FactorAnalysisOrdinal_Matrix.rst
1 .. _factoranalysisordinal-matrix-specification:
2
3 Factor Analysis Ordinal, Matrix Specification
4 =============================================
5
6 This example builds on the 'Factor Analysis, Matrix Specification' example, and extends it to the ordinal case.  We will discuss a common factor model with several indicators.  This example can be found in the following files:
7
8 * http://openmx.psyc.virginia.edu/svn/trunk/demo/OneFactorModelOrdinal_MatrixRaw.R
9 * http://openmx.psyc.virginia.edu/svn/trunk/demo/OneFactorModelOrdinal01_MatrixRaw.R
10
11 Common Factor Model
12 -------------------
13
14 The common factor model is a method for modeling the relationships between observed variables believed to measure or indicate the same latent variable. While there are a number of exploratory approaches to extracting latent factor(s), this example uses structural modeling to fit a confirmatory factor model. The model for any person and path diagram of the common factor model for a set of variables :math:`x_{1}` - :math:`x_{5}` are given below.
15
16 .. math::
17    :nowrap:
18    
19    \begin{eqnarray*} 
20    x_{ij} = \mu_{j} + \lambda_{j} * \eta_{i} + \epsilon_{ij}
21    \end{eqnarray*}
22
23 .. image:: graph/OneFactorModel5vt.png
24     :height: 2in
25
26 The path diagram above displays 16 parameters (represented in the arrows: 5 manifest variances, five manifest means, five factor loadings and one factor variance).  However, given we are dealing with ordinal data in this example, we are estimating thresholds rather than means, with nThresholds being one less the number of categories in the variables, here 3. Furthermore, we must constrain either the factor variance or one factor loading to a constant to identify the model and scale the latent variable.  In this instance, we chose to constrain the variance of the factor.  We also need to constrain the total variances of the manifest variables, as ordinal variables do not have a scale of measurement.  As such, this model contains 20 free parameters and is not fully saturated. 
27
28 Data
29 ^^^^
30
31 Our first step to running this model is to include the data to be analyzed. The data for this example were simulated in R.  Given the focus of this documentation is on OpenMx, we will not discuss the details of the simulation here, but we do provide the code so that the user can simulate data in a similar way.
32
33 .. code-block:: r
34
35     # Step 1: set up simulation parameters 
36     # Note: nVariables>=5, nThresholds>=1, nSubjects>=nVariables x nThresholds
37     # (maybe more) and model should be identified
38     
39     nVariables<-5
40     nFactors<-1
41     nThresholds<-3
42     nSubjects<-500
43     isIdentified <- function(nVariables,nFactors) 
44         as.logical(1+sign((nVariables*(nVariables-1)/2) 
45         - nVariables*nFactors + nFactors*(nFactors-1)/2))
46     # if this function returns FALSE then model is not identified, otherwise it is.
47     isIdentified(nVariables,nFactors)
48
49     loadings <- matrix(.7,nrow=nVariables,ncol=nFactors)
50     residuals <- 1 - (loadings * loadings)
51     sigma <- loadings %*% t(loadings) + vec2diag(residuals)
52     mu <- matrix(0,nrow=nVariables,ncol=1)
53     
54     # Step 2: simulate multivariate normal data
55
56     set.seed(1234)
57     continuousData <- mvrnorm(n=nSubjects,mu,sigma)
58
59     # Step 3: chop continuous variables into ordinal data 
60     # with nThresholds+1 approximately equal categories, based on 1st variable
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
69     # Step 4: make the ordinal variables into R factors
70
71     ordinalData <- mxFactor(as.data.frame(ordinalData),levels=c(1:(nThresholds+1)))
72
73     # Step 5: name the variables
74
75     bananaNames <- paste("banana",1:nVariables,sep="")
76     names(ordinalData) <- bananaNames
77     
78
79 Model Specification
80 ^^^^^^^^^^^^^^^^^^^
81
82 The following code contains all of the components of our model. Before running a model, the OpenMx library must be loaded into R using either the ``require()`` or ``library()`` function. All objects required for estimation (data, matrices, an expectation function, and a fit function) are included in their functions. This code uses the ``mxModel`` function to create an ``MxModel`` object, which we will then run.  We pre-specify a number of 'variables', namely the number of variables analyzed ``nVariables``, in this case 5, the number of factors ``nFactors``, here one, and the number of thresholds ``nthresholds``, here 3 or one less than the number of categories in the simulated ordinal variable.
83
84 .. code-block:: r
85
86     oneFactorThresholdModel <- mxModel("oneFactorThresholdModel",
87         mxMatrix(
88             type="Full", 
89             nrow=nVariables, 
90             ncol=nFactors, 
91             free=TRUE, 
92             values=0.2, 
93             lbound=-.99, 
94             ubound=.99, 
95             name="facLoadings"
96         ),
97         mxMatrix(
98             type="Unit", 
99             nrow=nVariables, 
100             ncol=1, 
101             name="vectorofOnes"
102         ),
103         mxAlgebra(
104             expression=vectorofOnes - (diag2vec(facLoadings %*% t(facLoadings))) , 
105             name="resVariances"
106         ),
107         mxAlgebra(
108             expression=facLoadings %*% t(facLoadings) + vec2diag(resVariances), 
109             name="expCovariances"
110         ),
111         mxMatrix(
112             type="Zero", 
113             nrow=1, 
114             ncol=nVariables, 
115             name="expMeans"
116         ),
117         mxMatrix(
118             type="Full", 
119             nrow=nThresholds, 
120             ncol=nVariables,
121             free=TRUE, 
122             values=.2,
123             lbound=rep( c(-Inf,rep(.01,(nThresholds-1))) , nVariables),
124             dimnames=list(c(), bananaNames),
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         mxData(
140             observed=ordinalData, 
141             type='raw'
142         ),
143         mxExpectationNormal(
144             covariance="expCovariances", 
145             means="expMeans", 
146             dimnames=bananaNames, 
147             thresholds="expThresholds"
148         ),
149         mxFitFunctionML()
150     )
151
152
153 This ``mxModel`` function can be split into several parts. First, we give the model a name "Common Factor ThresholdModel Matrix Specification".
154
155 The second component of our code creates an ``MxData`` object. The example above, reproduced here, first references the object where our data is, then uses the ``type`` argument to specify that this is raw data.
156
157 .. code-block:: r
158
159     mxData(
160         observed=ordinalData, 
161         type="raw"
162     )
163
164 The first ``mxMatrix`` statement declares a ``Full`` **nVariables x nFactors** matrix of factor loadings to be estimated, called "facLoadings", where the rows represent the dependent variables and the column(s) represent the independent variable(s).  The common factor model requires that one parameter (typically either a factor loading or factor variance) be constrained to a constant value. In our model, we will constrain the factor variance to 1 for identification, and let all the factor loadings be freely estimated.  Even though we specify just one start value of 0.2, it is recycled for each of the elements in the matrix.  Given the factor variance is fixed to one, and the variances of the observed variables are fixed to one (see below), the factor loadings are standarized, and thus must lie between -.99 and .99 as indicated by the ``lbound`` and ``ubound`` values.
165
166 .. code-block:: r
167
168     # factor loadings
169     mxMatrix(
170         type="Full", 
171         nrow=nVariables, 
172         ncol=nFactors, 
173         free=TRUE, 
174         values=0.2, 
175         lbound=-.99, 
176         ubound=.99, 
177         name="facLoadings"
178     )
179
180 Note that if ``nFactors>1``, we could add  a ``standardized`` ``mxMatrix`` to estimate the correlation between the factors.  Such a matrix automatically has 1's on the diagonal, fixing the factor variances to one and thus allowing all the factor loadings to be estimated.  In the current example, all the factor loadings are estimated which implies that the factor variance is fixed to 1.  Alternatively, we could add a ``symmetric`` **1x1** ``mxMatrix`` to estimates the variance of the factor, when one of the factor loadings is fixed.
181
182 As our data are ordinal, we further need to constrain the variances of the observed variables to unity.  These variances are made up of the contributions of the latent common factor and the residual variances.  The amount of variance explained by the common factor is obtained by squaring the factor loadings.  We subtract the squared factor loadings from 1 to get the amount explained by the residual variance, thereby implicitly fixing the variances of the observed variables to 1.  To do this for all variables simultaneously, we use matrix algebra functions.  We first specify a vector of One's by declaring a ``Unit`` **nVariables x 1** matrix called ``vectorofOnes``.  We need to subtract the squared factor loadings which are on the diagonal of the matrix multiplication of the factor loading matrix ``facLoadings`` and its transpose.  To extract those into squared factor loadings into a vector, we use the ``diag2vec`` function.  This new vector is subtracted from the ``vectorofOnes`` using an ``mxAlgebra`` statement to generate the residual variances, and named ``resVariances``.
183
184 .. code-block:: r
185
186     mxMatrix(
187         type="Unit", 
188         nrow=nVariables, 
189         ncol=1, 
190         name="vectorofOnes"
191     )
192     # residuals
193     mxAlgebra(
194         expression=vectorofOnes - (diag2vec(facLoadings %*% t(facLoadings))) , 
195         name="resVariances"
196     )
197
198 We then use the reverse function ``vec2diag`` to put the residual variances on the diagonal and add the contributions through the common factor from the matrix multipication of the factor loadings matrix and its transpose to obtain the formula for the expected covariances, aptly named ``expCovariances``.
199
200 .. code-block:: r
201
202     mxAlgebra(
203         expression=facLoadings %*% t(facLoadings) + vec2diag(resVariances), 
204         name="expCovariances"
205     )
206     
207 When fitting to ordinal rather than continuous data, we estimate thresholds rather than means.  The matrix of thresholds is of size **nThresholds x nVariables** where ``nThresholds`` is one less than the number of categories for the ordinal variable(s).  We still specify a matrix of means, however, it is fixed to zero.  An alternative approach is to fix the first two thresholds (to zero and one, see below), which allows us to estimate means and variances in a similar way to fitting to continuous data.  Let's first specify the model with zero means and free thresholds.
208
209 The means are specified as a ``Zero`` **1 x nVariables** matrix, called "expMeans".  A means matrix always contains a single row, and one column for every manifest variable in the model.
210
211 .. code-block:: r
212
213     # expected means
214     mxMatrix(
215         type="Zero", 
216         nrow=1, 
217         ncol=nVariables, 
218         name="expMeans"
219     )
220     
221 The mean of the factor(s) is also fixed to 1, which is implied by not including a matrix for it.  Alternatively, we could explicitly add a ``Full`` **1 x nFactors** ``mxMatrix`` with a fixed value of zero for the factor mean(s), named "facMeans".  
222
223 We estimate the ``Full`` **nThresholds x nVariables** matrix.  To make sure that the thresholds systematically increase from the lowest to the highest, we estimate the first threshold and the increments compared to the previous threshold by constraining the increments to be positive.  This is accomplished through some R algebra, concatenating `minus infinity` and (nThreshold-1) times .01 as the lower bound for the remaining estimates.  This matrix of ``thresholdDeviations`` is then pre-multiplied by a ``lower`` triangular matrix of ones of size **nThresholds x nThresholds**  to obtain the expected thresholds in increasing order in the ``thresholdMatrix``.
224
225 .. code-block:: r
226
227     mxMatrix(
228          type="Full", 
229          nrow=nThresholds, 
230          ncol=nVariables,
231          free=TRUE, 
232          values=.2,
233          lbound=rep( c(-Inf,rep(.01,(nThresholds-1))) , nVariables),
234          dimnames=list(c(), bananaNames),
235          name="thresholdDeviations"
236      )
237      mxMatrix(
238          type="Lower",
239          nrow=nThresholds,
240          ncol=nThresholds,
241          free=FALSE,
242          values=1,
243          name="unitLower"
244      )
245      # expected thresholds
246      mxAlgebra(
247          expression=unitLower %*% thresholdDeviations, 
248          name="expThresholds"
249      )
250
251 The final parts of this model are the expectation function and the fit function.  The choice of expectation function determines the required arguments.  Here we fit to raw ordinal data, thus we specify the matrices for the expected covariance matrix of the data, as well as the expected means and thresholds previously specified.  We use ``dimnames`` to map the model for means, thresholds and covariances onto the observed variables.
252
253 .. code-block:: r
254
255     mxExpectationNormal(
256         covariance="expCovariances", 
257         means="expMeans", 
258         dimnames=bananaNames, 
259         thresholds="expThresholds"
260     ),
261     mxFitFunctionML()
262
263 The free parameters in the model can then be estimated using full information maximum likelihood (FIML) for covariances, means and thresholds.  FIML is specified by using raw data with the ``mxFitFunctionML``.  To estimate free parameters, the model is run using the ``mxRun`` function, and the output of the model can be accessed from the ``@output`` slot of the resulting model.  A summary of the output can be reached using ``summary()``.
264
265 .. code-block:: r
266
267     oneFactorFit <- mxRun(oneFactorThresholdModel)
268
269     oneFactorFit@output
270
271     summary(oneFactorFit)
272     
273 As indicate above, the model can be re-parameterized such that means and variances of the observed variables are estimated similar to the continuous case, by fixing the first two thresholds.  This basically rescales the parameters of the model.  Below is the full script:
274
275 .. code-block:: r
276
277     oneFactorThresholdModel01 <- mxModel("oneFactorThresholdModel01",
278         mxMatrix(
279             type="Full", 
280             nrow=nVariables, 
281             ncol=nFactors, 
282             free=TRUE, 
283             values=0.2, 
284             lbound=-.99, 
285             ubound=2, 
286             name="facLoadings"
287         ),
288         mxMatrix(
289             type="Diag", 
290             nrow=nVariables, 
291             ncol=nVariables,
292             free=TRUE,
293             values=0.9,
294             name="resVariances"
295         ),
296         mxAlgebra(
297             expression=facLoadings %*% t(facLoadings) + resVariances, 
298             name="expCovariances"
299         ),
300         mxMatrix(
301             type="Full", 
302             nrow=1, 
303             ncol=nVariables,
304             free=TRUE,
305             name="expMeans"
306         ),
307         mxMatrix(
308             type="Full", 
309             nrow=nThresholds, 
310             ncol=nVariables,
311             free=rep( c(F,F,rep(T,(nThresholds-2))), nVariables), 
312             values=rep( c(0,1,rep(.2,(nThresholds-2))), nVariables),
313             lbound=rep( c(-Inf,rep(.01,(nThresholds-1))), nVariables),
314             dimnames=list(c(), bananaNames),
315             name="thresholdDeviations"
316         ),
317         mxMatrix(
318             type="Lower",
319             nrow=nThresholds,
320             ncol=nThresholds,
321             free=FALSE,
322             values=1,
323             name="unitLower"
324         ),
325         mxAlgebra(
326             expression=unitLower %*% thresholdDeviations, 
327             name="expThresholds"
328         ),
329         mxMatrix(
330             type="Unit",
331             nrow=nThresholds,
332             ncol=1,
333             name="columnofOnes"
334         ),
335         mxAlgebra(
336             expression=expMeans %x% columnofOnes,
337             name="meansMatrix"
338         ),
339         mxAlgebra(
340             expression=sqrt(t(diag2vec(expCovariances))) %x% columnofOnes,
341             name="variancesMatrix"
342         ),
343         mxAlgebra(
344             expression=(expThresholds - meansMatrix) / variancesMatrix,
345             name="thresholdMatrix"
346         ),
347         mxMatrix( 
348             type="Iden", 
349             nrow=nVariables, 
350             ncol=nVariables, 
351             name="Identity"
352         ),
353         mxAlgebra(
354             expression=solve(sqrt(Identity * expCovariances)) %*% facLoadings,
355             name="standFacLoadings"
356         ),
357         mxData(
358             observed=ordinalData, 
359             type='raw'
360         ),
361         mxExpectationNormal(
362             covariance="expCovariances", 
363             means="expMeans", 
364             dimnames=bananaNames, 
365             thresholds="expThresholds"
366         ),
367         mxFitFunctionML()
368     )
369
370 We will only highlight the changes from the previous model specification.  By fixing the first and second threshold to 0 and 1 respectively for each variable, we are now able to estimate a mean and a variance for each variable instead.  If we are estimating the variances of the observed variables, the factor loadings are no longer standardized, thus we relax the upper boundary on the factor loading matrix ``facLoadings`` to be 2.  The residual variances are now directly estimated as a ``Diagonal`` matrix of size ``nVariables x nVariables``, and given a start value higher than that for the factor loadings.  As the residual variances are already on the diagonal of the ``resVariances`` matrix, we no longer need to add the ``vec2diag`` function to obtain the ``expCovariances`` matrix.
371
372 .. code-block:: r
373
374     mxMatrix(
375         type="Full", 
376         nrow=nVariables, 
377         ncol=nFactors, 
378         free=TRUE, 
379         values=0.2, 
380         lbound=-.99, 
381         ubound=2, 
382         name="facLoadings"
383     )
384     mxMatrix(
385         type="Diag", 
386         nrow=nVariables, 
387         ncol=nVariables,
388         free=TRUE,
389         values=0.9,
390         name="resVariances"
391     )
392     mxAlgebra(
393         expression=facLoadings %*% t(facLoadings) + resVariances, 
394         name="expCovariances"
395     )
396     
397 Next, we now estimate the means for the observed variables and thus change the ``expMeans`` matrix to a ``Full`` matrix, and set it free.  The most complicated change happens to the matrix of ``thresholdDeviations``.  Its type and dimensions stay the same.  However, we now fix the first two thresholds, but allow the remainder of the thresholds (in this case, just one) to be estimated.  We use the R ``rep`` function to make this happen.  The ``values`` statement now has the fixed value of 0 for the first threshold, the fixed value of 1 for the second threshold, and the start value of .2 for the remaining threshold(s).  Finally, no change is required for the ``lbound`` matrix, which is still necessary to keep the estimated increments (third threshold and possible more) positive.
398
399 .. code-block:: r
400
401     mxMatrix(
402         type="Full", 
403         nrow=1, 
404         ncol=nVariables,
405         free=TRUE,
406         name="expMeans"
407     )
408     mxMatrix(
409         type="Full", 
410         nrow=nThresholds, 
411         ncol=nVariables,
412         free=rep( c(F,F,rep(T,(nThresholds-2))), nVariables), 
413         values=rep( c(0,1,rep(.2,(nThresholds-2))), nVariables),
414         lbound=rep( c(-Inf,rep(.01,(nThresholds-1))), nVariables),
415         dimnames=list(c(), bananaNames),
416         name="thresholdDeviations"
417     )
418
419 These are all the changes required to fit the alternative specification, which should give the same likelihood and goodness-of-fit statistics as the original one.  We have added some matrices and algebra to calculate the 'standardized' thresholds and factor loadings which should be equal to those obtained with the original specification.  To standardize the thresholds, the respective mean is subtracted from the thresholds, by expanding the means matrix to the same size as the threshold matrix.  The result is divided by the corresponding standard deviation.  To standardize the factor loadings, they are pre-multiplied by the inverse of the standard deviations.
420  
421 .. code-block:: r
422     
423     mxMatrix(
424         type="Unit",
425         nrow=nThresholds,
426         ncol=1,
427         name="columnofOnes"
428     )
429     mxAlgebra(
430         expression=expMeans %x% columnofOnes,
431         name="meansMatrix"
432     )
433     mxAlgebra(
434         expression=sqrt(t(diag2vec(expCovariances))) %x% columnofOnes,
435         name="variancesMatrix"
436     )
437     mxAlgebra(
438         expression=(expThresholds - meansMatrix) / variancesMatrix,
439         name="thresholdMatrix"
440     )
441     mxMatrix( 
442         type="Iden", 
443         nrow=nVariables, 
444         ncol=nVariables, 
445         name="Identity"
446     )
447     mxAlgebra(
448         expression=solve(sqrt(Identity * expCovariances)) %*% facLoadings,
449         name="facLoadingsMatrix"
450     )
451