changes to doc files
[openmx:openmx.git] / docs / source / AlternativeApproaches.rst
1 Alternative Approaches
2 ==================================
3
4 In the Beginner's Guide, we introduced the OpenMx framework for building and fitting models to data.  Here we will discuss some basic examples in more detail, providing a parallel treatment of the different formats and types of data and their associated functions to apply them to and the different methods and styles to build and fit models.  We will start with a single variable - univariate - example and generate data of different type and in different formats.  
5 Three types of data can be used: 
6
7 "(1)" continuous data, 
8 "(2)" ordinal data or 
9 "(3)" binary data,  
10
11 and the two supported data formats are 
12
13 "(i)" summary statistics, i.e. covariance matrices and possibly means, and 
14 "(ii)" raw data format.  
15
16 We will illustrate all of them, as arguments of functions may differ.  There are currently two supported methods for specifying models: 
17
18 "(a)" path specification and 
19 "(b)" matrix specification.  
20
21 Each example is presented using both approaches, so you can get a sense of their advantages/ disadvantages, and see which best fits your style.  In the 'path specification' model style you specify a model in terms of paths in a path diagram; the 'matrix specification' model style relies on matrices and matrix algebra to produce OpenMx code.  Furthermore, we will introduce the different styles available to specify models and data, namely 
22
23 "(A)" the piecewise style, 
24 "(B)" the stepwise style and
25 "(C)" the classic style.  
26
27 This will be done for a bivariate extension of the first example.
28
29 The code used to run the examples in this chapter is available here, and you may wish to access it while working through this manual. 
30
31 * http://openmx.psyc.virginia.edu/repoview/1/trunk/demo/AlternativeApproaches.R
32
33 * http://www.vipbg.vcu.edu/OpenMx/html/AlternativeApproaches.R
34
35
36 Introduction
37 ------------
38
39 Our first example fits a model to one variable - and is thus considered a univariate model.  The simplest model is one in which we estimate the variance of the variable.  Additionally we can estimate both the mean and the variance - a so called saturated model.  You might wonder why we use a modeling framework to estimate the mean and variance of a variable which we can do with anything from a calculator to using any basic spreadsheet or statistical program.  That is a fair question if you're just interested in the value of the mean and/or variance.  However, if you want to test, for example the mean of your variable is equal to a certain fixed value or to the mean of another variable or the same variable in another group, the modeling framework becomes relevant.  
40
41 Data Handling
42 -------------
43
44 You probably have your favorite data set ready to go, but before reading in data from an external file, we simulate a simple dataset directly in R, and use some of R's great capabilities.  As this is not an R manual, we just provide the code here with minimal explanation. There are several helpful sites for learning R, for instance http://www.statmethods.net/
45     
46 .. cssclass:: input
47 ..
48
49    OpenMx Code
50    
51    .. code-block:: r
52        
53         # Simulate Data
54         set.seed(100)
55         x <- rnorm (1000, 0, 1)
56         univData <- as.matrix(x)
57         dimnames(univData) <- list(NULL, "X")
58         summary(univData)
59         mean(univData)
60         var(univData)
61
62
63 The first line is a comment (starting with a #).  We set a seed for the simulation so that we generate the same data each time and get a reproducible answer.  We then create a variable **x** for 1000 subjects, with a mean of 0 and a variance of 1, using R's normal distribution function ``rnorm``.  We read the data in as a matrix into an object ``univData`` and give the variable a name *"X"* using the ``dimnames`` command.  We can easily produce some descriptive statistics in R using built-in functions ``summary``, ``mean`` and ``var``, just to make sure the data look like what we expect.  The output generated looks like this:   
64     
65     .. code-block:: r 
66
67          summary(univData)
68                 V1          
69          Min.   :-3.32078  
70          1st Qu.:-0.64970  
71          Median : 0.03690  
72          Mean   : 0.01681  
73          3rd Qu.: 0.70959  
74          Max.   : 3.30415  
75         > mean(univData)
76         [1] 0.01680509
77         > var(univData)
78               [,1]
79         [1,] 1.062112
80
81 For our second example, we will be fitting models to two variables which may be correlated.
82 The data used for the second example were generated using the multivariate normal function (``mvrnorm`` in the R package MASS).  The ``mvrnorm`` has three arguments: (i) sample size, (ii) vector of means, and (iii) covariance matrix.  We are simulating data on two variables named *X* and *Y* for 1000 individuals with means of zero, variances of one and a covariance of 0.5 using the following R code, and saved is as ``bivData``.  Note that we can now use the R function ``colMeans`` to generate the predicted means for the columns of our data frame and the ``cov`` to generate the observed covariance matrix.
83
84 .. cssclass:: input
85 ..
86
87    OpenMx Code
88    
89    .. code-block:: r
90
91          # Simulate Data
92          require(MASS)
93          set.seed(200)
94          rs=.5
95          xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
96          bivData <- xy
97          dimnames(bivData) <- list(NULL, c('X','Y'))
98          summary(bivData)
99          colMeans(bivData)
100          cov(bivData)
101
102 Notice that the simulated data are close to what we expected.
103    
104     ..  code-block:: r
105      
106          > summary(bivData)
107                 X                   Y            
108           Min.   :-3.296159   Min.   :-2.942561  
109           1st Qu.:-0.596177   1st Qu.:-0.633711  
110           Median :-0.010538   Median :-0.004139  
111           Mean   :-0.004884   Mean   : 0.032116  
112           3rd Qu.: 0.598326   3rd Qu.: 0.739236  
113           Max.   : 4.006771   Max.   : 4.173841  
114          >     colMeans(bivData)
115                     X            Y 
116          -0.004883811  0.032116480 
117          >     cov(bivData)
118                    X         Y
119          X 0.9945328 0.4818317
120          Y 0.4818317 1.0102951
121                   
122
123 Data Types
124 ^^^^^^^^^^
125
126 Continuous Data
127 +++++++++++++++
128
129 The data we simulated are continuous in nature and follow a normal distribution.  This can easily be verified by R's excellent graphical capabilities.  Here we show the R code and a basic histogram of the *univData* we generated.  
130
131 .. cssclass:: input
132 ..
133
134    OpenMx Code
135    
136    .. code-block:: r
137        
138          hist(univData)
139     
140 .. image:: graph/histogram_testData.png
141
142 This is the ideal type of data, as many of the models we fit to them assume that the data are normally distributed.  However, reality is often different and it might be necessary to apply a transformation to the original data to better approximate a normal distribution.  When there are 15 or more possible values for the variable of interest, it is appropriate to treat them as continuous.  Note that although the simulated data have many more than 15 different values, values are binned to simplify the graph.
143
144 Continuous data can be summarized by their mean and standard deviation.  Two or more variables are summarized by a vector of means and a covariance matrix which provides information on the variance of each of the variables as well as the covariances between the variables.
145
146 Categorical Data
147 ++++++++++++++++
148
149 A lot of variables, however, are not measured on a continuous scale, but using a limited number of categories.  If the categories are ordered in a logical way, we refer to them as *ordinal* variables and often assume that the underlying construct follows a normal distribution.  This assumption can actually be tested for any ordinal variable with a minimum of three categories, when more than one variable is available or the same variable is measured in related individuals or over time.
150
151 Categorical data contain less information than continuous data, and are summarized by thresholds which predict the proportion of individuals in a specific category.  As the sum of the proportions of each of the categories adds up to one, there is no information about the variance.  The relationship between two or more variables can be summarized in a correlation matrix.  Rather than estimating two (or more) thresholds and a correlation(s), one could fix the first threshold to zero and the second to one and estimate the means and covariance matrices instead, which can be interpreted in the same way as for continuous variables.  The estimated proportion in each of the categories can then be calculated by rescaling the statistics.
152
153 Often, unfortunately, variables are only measured with two categories (Yes/No, affected/unaffected, etc.) and are called *binary* variables.  The only statistic to be estimated in the univariate case is the threshold and no information is available about the variance.  With two or more variables, their relationship is also summarized in a correlation matrix.
154
155 The power of your study is directly related to the type of variable analyzed, and is typically higher for continuous variables compared to categorical variables, with ordinal variables providing more power than binary variables.  Whenever possible, use continuous variables or at least ordinal variables.
156
157 As a lot of real data are only available as categorical variables, we will generate both an ordinal and a binary variable from the simulated continuous variable in 'univData'.  The code below uses the ``cut`` and ``breaks`` commands to cut the continuous variable into 5 ordered categories.
158
159 .. cssclass:: input
160 ..
161
162    OpenMx Code
163    
164    .. code-block:: r
165        
166          univDataOrd <- data.frame(X=cut(univData[,1], breaks=5, ordered_result=T, 
167                                    labels=c(0,1,2,3,4)) )
168          table(univDataOrd)
169
170 A summary of the resulting data set can be generated as follows:
171
172     ..  code-block:: r
173     
174         > table(univDataOrd)
175         univDataOrd
176           0   1   2   3   4 
177          28 216 483 244  29
178        
179     
180 A similar approach could be used to create a binary variable.  However, here we show an alternative way to generate a binary variable using a specific cutoff using the ``ifelse`` command.  We will assign a value of 1 when the value of our original standardized continuous variable is above 0.5; otherwise a value of 0 will be assigned.
181
182 .. cssclass:: input
183 ..
184
185    OpenMx Code
186    
187    .. code-block:: r
188        
189          univDataBin <- data.frame(X=ifelse(univData[,1] >.5,1,0))
190          table(univDataBin)
191
192 The resulting data set table looks as follows: 
193
194     ..  code-block:: r
195     
196         > table(univDataBin)
197         univDataBin
198           0   1 
199         680 320
200         
201 We will go through the same steps to generate ordinal and binary data from the simulated bivariate data.  Given we need to repeat the same statement for the two variables, we employ a ``for`` statement.
202
203 .. cssclass:: input
204 ..
205
206    OpenMx Code
207    
208    .. code-block:: r
209        
210         bivDataOrd <- data.frame(bivData)
211         for (i in 1:2) { bivDataOrd[,i] <- cut(bivData[,i], breaks=5, ordered_result=T, 
212                                                labels=c(0,1,2,3,4)) }
213            table(bivDataOrd[,1],bivDataOrd[,2])
214         bivDataBin <- data.frame(bivData)
215         for (i in 1:2) { bivDataBin[,i] <- ifelse(bivData[,i] >.5,1,0) }
216            table(bivDataBin[,1],bivDataBin[,2])
217         
218
219 Data Formats
220 ^^^^^^^^^^^^
221
222 Raw Data
223 ++++++++
224
225 To make these data available for statistical modeling in OpenMx, we need to create an "MxData" object which is accomplished with the ``mxData`` function.  Remember to load the OpenMx package first.
226
227 .. cssclass:: input
228 ..
229
230    OpenMx Code
231    
232    .. code-block:: r
233        
234         require(OpenMx)
235         obsRawData <- mxData( observed=univData, type="raw" )
236         selVars <- "X"
237
238 First, we read the data matrix in with the ``observed`` argument.  Then, we tell OpenMx what format or type the data is in, in this case we're reading in the raw data.  We save this MxData object as *obsRawData*.  As later on, we need to be able to map our data onto the model, we typically create a vector with the variable labels of the variable(s) we are analyzing.  To make our scripts more readable, we use consistent names for objects - something you can decide to copy or change according to your preferences - and we use *selVars* for the variables we select for analysis.  In this example, it is a single variable *X*.
239
240     ..  code-block:: r
241     
242         > obsRawData
243         MxData 'data' 
244         type : 'raw' 
245         numObs : '1000' 
246         Data Frame or Matrix : 
247                             X
248            [1,] -5.021924e-01
249            [2,]  1.315312e-01
250            ....
251            [1000,] -2.141428e+00
252            Means : NA 
253            Acov : NA 
254            Thresholds : NA 
255         
256 A look at this newly created object shows that it was given the  ``name`` *data*, which is done by default.  It has the ``type`` that we specified, and ``numObs`` are automatically counted for us.  The actual data for the variable *X* are then listed; we only show the first two values.
257
258 In a similar manner we create a MxData object for the second example.  We read in the ``observed`` 'bivData', and indicate the ``type`` as raw.  We refer to this object as *obsBivData*.
259
260 .. cssclass:: input
261 ..
262
263    OpenMx Code
264    
265    .. code-block:: r
266        
267          obsBivData <- mxData( observed=bivData, type="raw" )
268
269 If we want to fit models to categorical data, we need to read in the ordinal or binary data.  However, when your data are ordinal or binary, OpenMx expects them to be 'ordered factors'.  To ensure that your data have the appropriate format, it is recommended/required to apply the ``mxFactor`` command to the categorical variables, where the ``x`` argument reads in a vector of data or a data.frame, and ``levels`` expects a vector of possible values for those data.  We save the resulting objects as *univDataOrdF* and *univDataBinF*, or *bivDataOrdF* and *bivDataBinF* for the corresponding data in the second example.
270
271 .. cssclass:: input
272 ..
273
274    OpenMx Code
275    
276    .. code-block:: r
277        
278           univDataOrdF <- mxFactor( x=univDataOrd, levels=c(0:4) )
279           univDataBinF <- mxFactor( x=univDataBin, levels=c(0,1) )
280           bivDataOrdF  <- mxFactor( x=bivDataOrd, levels=c(0:4) )
281           bivDataBinF  <- mxFactor( x=bivDataBin, levels=c(0,1) )
282
283
284 Next, we generate the corresponding MxData objects.
285
286 .. cssclass:: input
287 ..
288
289    OpenMx Code
290    
291    .. code-block:: r
292        
293          obsRawDataOrd <- mxData( observed=univDataOrdF, type="raw" )
294          obsRawDataBin <- mxData( observed=univDataBinF, type="raw" )
295          obsBivDataOrd <- mxData( observed=bivDataOrdF, type="raw" )
296          obsBivDataBin <- mxData( observed=bivDataBinF, type="raw" )
297          
298
299
300 Summary Stats
301 +++++++++++++
302
303 Covariances
304 ...........
305
306 While analyzing raw data is the standard in most statistical modeling these days, this was not the case in a previous generation of computers, which could only deal with summary statistics.  As fitting models to summary statistics still is much faster then using raw data (unless your data set is small), it is sometimes useful for didactic purposes.  Furthermore, sometimes one has access only to the summary statistics.  In the case where the dataset is complete, in other words there are no missing data, there is no advantage to using raw data.  For our example, we can easily create a covariance matrix based on our data set by using R's ``var()`` function, in the case of analyzing a single variable, or ``cov()`` function, when analyzing more than one variable.  This can be done prior to or directly when creating the MxData object.   Its first argument, ``observed``, reads in the data from an R matrix or data.frame, with the ``type`` given in the second argument, followed by the ``numObs`` argument which is necessary when reading in summary statistics.
307
308 .. cssclass:: input
309 ..
310
311    OpenMx Code
312    
313    .. code-block:: r
314        
315          univDataCov <- var(univData)
316          obsCovData  <- mxData( observed=univDataCov, type="cov", numObs=1000 )
317         
318 or 
319
320 .. cssclass:: input
321 ..
322
323    OpenMx Code
324    
325    .. code-block:: r
326        
327          obsCovData <- mxData( observed=var(univData), type="cov", numObs=1000 )
328
329 Given our first example has only one variable, we use the ``var()`` function (as there is no covariance for a single variable).  When summary statistics are used as input, the number of observations (``numObs``) needs to be supplied.  The resulting MxData object looks as follows:
330
331     ..  code-block:: r
332     
333         > obsCovData
334         MxData 'data' 
335         type : 'cov' 
336         numObs : '1000' 
337         Data Frame or Matrix : 
338                  X
339         X 1.062112
340         Means : NA
341         Acov : NA 
342         Thresholds : NA
343     
344 The differences with the previous data objects are that the type is now 'cov' and the actual data frame is now a single value, the variance of the 1000 data points.
345
346 Covariances and Means
347 .....................
348
349 In addition to the observed covariance matrix, a fourth argument ``means`` can be added for the vector of observed means from the data, calculated using the R ``colMeans`` command.
350
351 .. cssclass:: input
352 ..
353
354    OpenMx Code
355    
356    .. code-block:: r
357        
358          obsCovMeanData <- mxData( observed=var(univData),  type="cov", numObs=1000, 
359                                    means=colMeans(univData) )
360
361 You can verify that the new *obsCovMeanData* object now has a value for the observed means as well.
362
363 For the second, bivariate example the only change we'd have to make - besides reading in the *bivData* - is the use of ``cov`` instead of ``var`` to generate the object for the observed covariance matrix.
364
365 Correlations
366 ............
367
368 To analyze categorical data, we can also fit the models to summary statistics, in this case, correlation matrices, as indicated by using the ``cor()`` R command to generate them and by the ``type=`` =cor, which also requires the ``numObs`` argument to indicate how many observations (data records) are in the dataset.
369  
370 .. cssclass:: input
371 ..
372
373    OpenMx Code
374    
375    .. code-block:: r
376        
377          obsOrdData <- mxData( observed=cor(univDataOrdF), type="cor", numObs=1000 )
378
379
380 We will start by fitting a simple univariate model to the continuous data and then show which changes have to be made when dealing with ordinal or binary variables.  For the continuous data example, we will start with fitting the model to the summary statistics prior to fitting to raw data and show their equivalence (in the absence of missing data).
381
382
383 Model Handling
384 --------------
385
386 Path Method 
387 ^^^^^^^^^^^
388
389 Summary Stats 
390 +++++++++++++
391
392 If we have data on a single variable X summarized in its variance, the basic univariate model will simply estimate the variance of the variable X.  We call this model saturated because there is a free parameter corresponding to each and every observed statistic.  Here we have covariance matrix input only, so we can estimate one variance.  This model can be represented by the following path diagram:
393
394 .. image:: graph/UnivariateSaturatedModelNoMean.png
395
396 Model Building
397 ..............
398
399 When using the path specification, it is easiest to work from the path diagram.  Assuming you are familiar with path analysis (*for those who are not, there are several excellent introductions, see [LI1986]*), we have a box for the observed/manifest variable *X*, and one double headed arrow, labeled /sigma^2_x.  To indicate which variable we are analyzing, we use the ``manifestVars`` argument, which takes a vector of labels.  In this example, we are selecting one variable, which we pre-specified in the *selVars* object.
400
401 .. cssclass:: input
402 ..
403
404    OpenMx Code
405    
406    .. code-block:: r
407        
408          selVars   <- c("X")
409          manifestVars=selVars
410
411 We have already built the MxData object above, so here we will build the model by specifying the relevant paths.  Our first model only has one path which has two arrows and goes from the variable X to the variable X.  That path represents the variance of X which we aim to estimate.  Let's see how this translates into the ``mxPath`` object.
412
413 The ``mxPath`` command indicates where the path originates (``from``) and where it ends (``to``).  If the ``to`` argument is omitted, the path ends at the same variable where it started.  The ``arrows`` argument distinguishes one-headed arrows (if arrows=1) from two-headed arrows (if arrows=2).  The ``free`` command is used to specify which elements are free or fixed with a ``TRUE`` or ``FALSE`` option.  If the ``mxPath`` command creates more than one path, a single ``T`` implies that all paths created here are free.  If some of the paths are free and others fixed, a list is expected.  The same applies for the ``values`` command which is used to assign starting values or fixed final values, depending on the corresponding 'free' status.  Optionally, lower and upper bounds can be specified (using ``lbound`` and ``ubound``, again generally for all the paths or specifically for each path).  Labels can also be assigned using the ``labels`` command which expects as many labels (in quotes) as there are elements.  Thus for our example, we specify only a ``from`` argument, as the double-headed arrow (``arrows`` =2) goes back to *X*.  This path is estimated (``free`` =TRUE), and given a start value of 1 (``values`` =1) and has to be positive (``lbound`` =.01).  Finally we assign it a label (``labels`` ="vX").  The generated MxPath object is called *expVariance*.
414
415 .. cssclass:: input
416 ..
417
418    OpenMx Code
419    
420    .. code-block:: r
421        
422          expVariance <- mxPath(
423              from=c("X"), arrows=2, 
424              free=TRUE, 
425              values=1, 
426              lbound=.01, 
427              labels="vX"
428          ),
429
430 Note that all arguments could be listed on one line; in either case they are separated by comma's:
431
432 .. cssclass:: input
433 ..
434
435    OpenMx Code
436    
437    .. code-block:: r
438
439          expVariance <- mxPath( from=c("X"), arrows=2, 
440                                 free=TRUE, values=1, lbound=.01, labels="vX" )
441
442 The resulting MxPath object looks as follows:
443
444     ..  code-block:: r
445     
446        > expVariance
447         mxPath 
448         $from:  'X' 
449         $to:  'X' 
450         $arrows:  2 
451         $values:  1 
452         $free:  TRUE 
453         $labels:  vX 
454         $lbound:  0.01 
455         $ubound:  NA
456         $connect:  single   
457     
458 To evaluate the model that we have built, we need an expectation and a fit function that obtains the best solution for the model given the data.  When using the path specification, both are automatically generated by invoking the ``type=RAM`` argument in the model.  The 'RAM' objective function has a predefined structure.
459
460 .. cssclass:: input
461 ..
462
463    OpenMx Code
464    
465    .. code-block:: r
466        
467          type="RAM"
468
469 Internally, OpenMx translates the paths into RAM notation in the form of the matrices **A**, **S**, and **F** [see RAM1990].  Before we can 'run' the model through the optimizer, we need to put all the arguments into an MxModel using the ``mxModel`` command.  Its first argument is a name, and therefore is in quotes.  We then add all the arguments we have built so far, including the list of variables to be analyzed in ``manifestVars``, the MxData object, and the predicted model specified using paths.
470
471 .. cssclass:: input
472 ..
473
474    OpenMx Code
475    
476    .. code-block:: r
477        
478          univSatModel1 <- mxModel("univSat1", manifestVars=selVars, obsCovData, 
479                                   expVariance, type="RAM" )
480
481 We can inspect the MxModel object generated by this statement.
482
483     ..  code-block:: r
484     
485         > univSatModel1
486         MxModel 'univSat1' 
487         type : RAM 
488         $matrices : 'A', 'S', and 'F' 
489         $algebras :  
490         $constraints :  
491         $intervals :  
492         $latentVars : none
493         $manifestVars : 'X' 
494         $data : 1 x 1 
495         $data means : NA
496         $data type: 'cov' 
497         $submodels :  
498         $expectation : MxExpectationRAM 
499         $fitfunction : MxFitFunctionML 
500         $compute : NULL 
501         $independent : FALSE 
502         $options :  
503         $output : FALSE 
504     
505 Note that only the relevant arguments have been updated, and that the path information has been stored in the **A**, **S**, and **F** matrices.  The free parameter for the variance "vX" ends up in the **S** matrix which holds the symmetric (double-headed) paths.  Here we print the details for this **S** matrix:
506
507     ..  code-block:: r
508     
509        > univSatModel1$matrices$S
510         SymmMatrix 'S' 
511     
512         $labels
513           X   
514         X "vX"
515     
516         $values
517           X
518         X 1
519     
520         $free
521              X
522         X TRUE
523     
524         $lbound
525              X
526         X 0.01
527     
528         $ubound: No upper bounds assigned.
529     
530     
531 Model Fitting
532 .............
533
534 So far, we have specified the model, but nothing has been evaluated.  We have 'saved' the specification in the object ``univSatModel1``.  This object is evaluated when we invoke the ``mxRun`` command with the MxModel object as its argument.
535
536 .. cssclass:: input
537 ..
538
539    OpenMx Code
540    
541    .. code-block:: r
542        
543          univSatFit1 <- mxRun(univSatModel1)
544
545 You can verify that the arguments of the *univSatModel1* and *univSatFit1* look mostly identical.  What we expect to be updated with the estimated value of variance is the element of the **S** matrix, which we can output as follows:
546
547     ..  code-block:: r
548     
549         > univSatFit1$matrices$S$values
550                  X
551         X 1.062112
552         
553 An alternative form of extracting values from a matrix is:
554
555     ..  code-block:: r
556     
557         > univSatFit1[['S']]$values
558                  X
559         X 1.062112
560     
561 There are actually a variety of ways to generate output.  We will promote the use of the ``mxEval`` command, which takes two arguments: an ``expression`` and a ``model`` object.  The ``expression`` can be a matrix or algebra  defined in the model, new calculations using any of these matrices/algebras of the model, the objective function, etc.  Here we use ``mxEval`` to simply list the values of the **S** matrix, which formats the output slightly differently as a typical R matrix object, and call it *EC1*, short for the expected covariance:
562
563     ..  code-block:: r
564
565         EC1 <- mxEval(S, univSatFit1)
566         >        EC1
567                      X
568             X 1.062112
569     
570     
571 We can then use any regular R function in the ``mxEval`` command to generate derived fit statistics, some of which are built in as standard.  When fitting to covariance matrices, the saturated likelihood can be easily obtained and subtracted from the likelihood of the data to obtain a Chi-square goodness-of-fit.  The saturated likelihood, here named 'SL1' is obtained from the ``$output$Saturated`` argument of the fitted object *univSatFit1* which contains a range of statistics.  We get the likelihood of the data, here referred to as *LL1*, from the ``$output$fit`` argument of the fitted object *univSatFit1*.
572
573 .. cssclass:: input
574 ..
575
576    OpenMx Code
577    
578    .. code-block:: r
579        
580          SL1  <- univSatFit1$output$Saturated
581          LL1  <- univSatFit1$output$fit
582          Chi1 <- LL1-SL1
583
584 The output of these objects like as follows
585
586     ..  code-block:: r
587     
588         > SL1
589         [1] 1059.199
590         > LL1
591         [1] 1059.199
592         > Chi1
593         [1] 0
594         
595         
596 An alternative to requesting specific output is to generate the default summary of the model, which can be done with the ``summary`` function, and can also be saved in another R object, i.e. *univSatSumm1*.
597
598 .. cssclass:: input
599 ..
600
601    OpenMx Code
602    
603    .. code-block:: r
604        
605          summary(univSatFit1)
606          univSatSumm1 <- summary(univSatFit1)
607
608 This output includes a summary of the data (if available), a list of all the free parameters with their estimates (if the model contains free parameters), their confidence intervals (if requested), a list of goodness-of-fit statistics, and a list of job statistics (timestamps and OpenMx version).
609
610     ..  code-block:: r
611     
612        > univSatSumm1
613         data:
614         $univSat1.data
615         $univSat1.data$cov
616                  X
617         X 1.062112
618     
619     
620         free parameters:
621           name matrix row col Estimate  Std.Error Std.Estimate     Std.SE lbound ubound
622         1   vX      S   X   X 1.062112 0.04752282            1 0.04474372   0.01              
623     
624         observed statistics:  1 
625         estimated parameters:  1 
626         degrees of freedom:  0 
627         -2 log likelihood:  1059.199 
628         saturated -2 log likelihood:  1059.199 
629         number of observations:  1000 
630         chi-square:  0 
631         p:  1 
632         Information Criteria: 
633             df Penalty Parameters Penalty Sample-Size Adjusted
634         AIC          0           2.000000                   NA
635         BIC          0           6.907755             3.731699
636         CFI: NaN 
637         TLI: NaN 
638         RMSEA:  NA 
639         timestamp: 2014-04-02 18:41:35 
640         frontend time: 0.09399414 secs 
641         backend time: 0.007524967 secs 
642         independent submodels time: 5.602837e-05 secs 
643         wall clock time: 0.1015751 secs 
644         cpu time: 0.1015751 secs 
645         openmx version number: 999.0.0-3160 
646     
647
648 In addition to providing a covariance matrix as input data, we could add a means vector.  As this requires a few minor changes, lets highlight those.  The path diagram for this model, now including means (path from triangle of value 1) is as follows:
649
650 .. image:: graph/UnivariateSaturatedModel.png
651
652 We have to specify one additional ``mxPath`` command for the means.  In the path diagram, the means are specified by a triangle which has a fixed value of one, reflected in the ``from="one"`` argument, with the ``to`` argument referring to the variable whose mean is estimated.  Note that paths for means are always single headed.  We will save this path as the R object *expMean*.
653
654 .. cssclass:: input
655 ..
656
657    OpenMx Code
658    
659    .. code-block:: r
660        
661          expMean <- mxPath(from="one", to="X", arrows=1, free=TRUE, values=0, labels="mX")
662
663 This new path adds one additional parameter, called 'mX'.
664
665     ..  code-block:: r
666     
667         >  expMean
668         mxPath 
669         $from:  'one' 
670         $to:  'X' 
671         $arrows:  1 
672         $values:  0 
673         $free:  TRUE 
674         $labels:  mX 
675         $lbound:  NA 
676         $ubound:  NA
677         $connect:  single 
678     
679 The other required change is in the ``mxData`` command, which now takes a fourth argument ``means`` for the vector of observed means from the data, calculated using the R ``colMeans`` command.
680
681 .. cssclass:: input
682 ..
683
684    OpenMx Code
685    
686    .. code-block:: r
687        
688          obsCovMeanData <- mxData( observed=var(univData), type="cov", numObs=1000, 
689                                    means=colMeans(univData) )
690
691 As this new object will simply be added to the previous model, we can build onto our existing model.  Therefore, instead of using the first argument for the name, we use it in its other capacity, namely as the name of a previously defined MxModel object that is being modified.  In this case, we start with the previous model *univSatModel1*, which becomes the first argument of our new model *univSatModel1M*.  To change the name of the object, we add a ``name`` argument.  Note that the default order of arguments can be changed by adding the argument's syntax name.  We then add the new argument for the expected means, as well as the modified MxData object.
692
693 .. cssclass:: input
694 ..
695
696    OpenMx Code
697    
698    .. code-block:: r
699        
700          univSatModel1M <- mxModel(univSatModel1, name="univSat1M", expMean, obsCovMeanData )
701     
702 Note the following changes in the modified MxModel below.  First, the name is changed to 'univSat1M'.  Second, an additional matrix **M** was generated for the expected means vector.  Third, observed means were added, here referred to as '$data means'.
703
704     ..  code-block:: r
705     
706         >       univSatModel1M
707          MxModel 'univSat1M' 
708          type : RAM 
709          $matrices : 'A', 'S', 'F', and 'M' 
710          $algebras :  
711          $constraints :  
712          $intervals :  
713          $latentVars : none
714          $manifestVars : 'X' 
715          $data : 1 x 1 
716          $data means : 1 x 1 
717          $data type: 'cov' 
718          $submodels :  
719          $expectation : MxExpectationRAM 
720          $fitfunction : MxFitFunctionML 
721          $compute : NULL  
722          $independent : FALSE 
723          $options :  
724          $output : FALSE
725     
726 When a mean vector is supplied and a parameter added for the estimated mean, the RAM matrices **A**, **S** and **F** are augmented with an **M** matrix which can be extracted from the output in a similar way as the expected variance before, and is called 'EM1', short for expected mean.
727
728 .. cssclass:: input
729 ..
730
731    OpenMx Code
732    
733    .. code-block:: r
734        
735          univSatFit1M  <- mxRun(univSatModel1M)
736          EM1M          <- mxEval(M, univSatFit1M) 
737          univSatSumm1M <- summary(univSatFit1M)
738
739 The new summary object *univSatSumm1M* is different from the previous one in the following ways: the observed data means were added, an extra free parameter is listed and estimated, thus the fit statistics are updated.  Notice, however, that the likelihood of both models is the same.  (We have cut part of the summary that is not relevant here.)
740
741     ..  code-block:: r
742     
743         > univSatSumm1M
744         data:
745         $univSat1M.data
746         $univSat1M.data$cov
747                  X
748         X 1.062112
749     
750         $univSat1M.data$means
751                       X
752         [1,] 0.01680509
753     
754     
755         free parameters:
756           name matrix row col   Estimate  Std.Error Std.Estimate     Std.SE lbound ubound
757         1   vX      S   X   X 1.06211141 0.04752281            1 0.04474372   0.01       
758         2   mX      M   1   X 0.01680503 0.03259006           NA         NA                            
759     
760         observed statistics:  2 
761         estimated parameters:  2 
762         degrees of freedom:  0 
763         -2 log likelihood:  1059.199 
764         saturated -2 log likelihood:  1059.199 
765         number of observations:  1000 
766         chi-square:  8.867573e-12 
767         p:  0 
768         Information Criteria: 
769               df Penalty Parameters Penalty Sample-Size Adjusted
770         AIC 8.867573e-12            4.00000                   NA
771         BIC 8.867573e-12           13.81551             7.463399
772     
773
774 Raw Data 
775 ++++++++
776
777 Instead of fitting models to summary statistics, it is now popular to fit models directly to the raw data and using full information maximum likelihood (FIML).  Doing so requires specifying not only a model for the covariances, but also one for the means, just as in the case of fitting to covariance matrices and mean vectors described above. 
778
779 The only change required is in the MxData object, *obsRawData* defined above, which reads the raw data in directly from an R matrix or a data.frame into the ``observed`` first argument, and has ``type="raw"`` as its second argument.  A nice feature of OpenMx is that existing models can be easily modified.  Here we will start from the saturated model estimating covariances and means from summary statistics, namely *univSatModel1M*, as both expected means and covariances have to be modeled when fitting to raw data.
780
781 .. cssclass:: input
782 ..
783
784    OpenMx Code
785    
786    .. code-block:: r
787        
788          univSatModel2 <- mxModel(univSatModel1M, obsRawData )
789
790 The resulting model can be run as usual using ``mxRun``:
791
792 .. cssclass:: input
793 ..
794
795    OpenMx Code
796    
797    .. code-block:: r
798        
799          univSatFit2  <- mxRun(univSatModel2)
800          univSatSumm2 <- summary(univSatFit2)
801          EM2          <- mxEval(M, univSatFit2) 
802          EC2          <- mxEval(S, univSatFit2)
803          LL2          <- univSatFit2$output$fit
804
805 Note that the estimates for the expected means, as well as the expected covariance matrix are exactly the same as before, as we have no missing data.
806
807     ..  code-block:: r
808     
809         >        EM2
810                       X
811         [1,] 0.01680499
812         >        EC2
813                  X
814         X 1.061049
815         >        LL2
816         [1] 2897.135
817
818 The estimates for the predicted mean and covariance matrix are exactly the same as those obtained when fitting to summary statistics.  The likelihood, however, is different.??
819
820     ..  code-block:: r
821     
822        > univSatSumm2
823        data:
824        $univSat1M.data
825               X           
826         Min.   :-3.32078  
827         1st Qu.:-0.64970  
828         Median : 0.03690  
829         Mean   : 0.01681  
830         3rd Qu.: 0.70959  
831         Max.   : 3.30415  
832     
833        free parameters:
834          name matrix row col   Estimate  Std.Error Std.Estimate     Std.SE lbound ubound
835        1   vX      S   X   X 1.06104923 0.04745170            1 0.04472149   0.01       
836        2   mX      M   1   X 0.01680499 0.03257418           NA         NA                      
837     
838        observed statistics:  1000 
839        estimated parameters:  2 
840        degrees of freedom:  998 
841        -2 log likelihood:  2897.135 
842        saturated -2 log likelihood:  NA 
843        number of observations:  1000 
844        chi-square:  NA 
845        p:  NA 
846        Information Criteria: 
847            df Penalty Parameters Penalty Sample-Size Adjusted
848        AIC   901.1355           2901.135                   NA
849        BIC -3996.8043           2910.951             2904.599
850         
851     
852
853 Matrix Method
854 ^^^^^^^^^^^^^
855
856 The next example replicates these models using matrix-style coding.  In addition to the ``mxData``  and ``mxModel`` commands which were introduced before, the code to specify the model includes three new commands, (i) ``mxMatrix``, and (ii) ``mxExpectationNormal`` and ``mxFitFunctionML()``.
857
858 Summary Stats
859 +++++++++++++
860
861 Covariances 
862 ...........
863
864 Starting with the model fitted to the summary covariance matrix, the ``mxData`` is identical to that used in path style models, as is the case for all the corresponding models specified using paths or matrices. 
865
866 To specify the model, we now create a matrix for the expected covariance matrix using the ``mxMatrix`` command.  The first argument is its ``type``, symmetric for a covariance matrix.  The second and third arguments are the number of rows (``nrow``) and columns (``ncol``) – one each for a univariate model.  The ``free`` and ``values`` parameters work as in the path specification.  If only one element is given, it is applied to all elements of the matrix.  Alternatively, each element can be assigned its free/fixed status and starting value with a list command.  Note that in the current example, the matrix is a simple **1x1** matrix, but that will change rapidly in the later examples.
867
868 .. cssclass:: input
869 ..
870
871    OpenMx Code
872    
873    .. code-block:: r
874        
875          expCovMat <- mxMatrix( type="Symm", nrow=1, ncol=1, 
876                                 free=TRUE, values=1, name="expCov" )
877     
878 The resulting MxMatrix object *expCovMat* looks as follows.  Note that the starting value for the free parameter is 1 and that optionally labels can be assigned for the rows and columns of the matrix and lower and upper bounds can be assigned to limit the parameter space for the estimation:
879
880     ..  code-block:: r
881     
882        > expCovMat
883         SymmMatrix 'expCov' 
884     
885         $labels: No labels assigned.
886     
887         $values
888              [,1]
889         [1,]    1
890     
891         $free
892              [,1]
893         [1,] TRUE
894     
895         $lbound: No lower bounds assigned.
896     
897         $ubound: No upper bounds assigned.
898         
899 To link the model for the covariance matrix to the data, an ``mxExpectation`` needs to be specified which will then be evaluated with an ``mxFitFunction``.  The ``mxExpectationNormal`` command  takes two arguments, ``covariance`` to hold the expected covariance matrix (which we named "expCov" above using the ``mxMatrix`` command), and ``dimnames`` which allow the mapping of the observed data to the expected covariance matrix, i.e. the model.  ``mxFitFunctionML()`` will invoke the maximum likelihood ('ML'), to obtain the best estimates for the free parameters.
900
901 .. cssclass:: input
902 ..
903
904    OpenMx Code
905    
906    .. code-block:: r
907        
908          expectCov <- mxExpectationNormal( covariance="expCov", dimnames=selVars )
909          funML <- mxFitFunctionML()
910          
911 The internal name of an MxExpectationNormal object is by default *expectation* and that for an MxFitFunctionML object is by default *fitfunction*.  We can thus inspect these two objects by using the names of the resulting objects, here *expCovFun* and *ML* as shown below. The result of applying the fit function is not yet computed and thus reported as *<0 x 0 matrix>*; its arguments will change after running the model successfully.
912     
913     ..  code-block:: r
914     
915          > expectCov
916          MxExpectationNormal 'expectation' 
917          $covariance : 'expCov' 
918          $means : NA 
919          $dims : 'X' 
920          $thresholds : NA 
921          $threshnames : 'X'
922          
923          > funML
924          MxFitFunctionML 'fitfunction' 
925          $vector : FALSE 
926          <0 x 0 matrix>
927     
928 We can then simply combine the appropriate elements into a new model and fit it in the usual way to the data.  Please note that within the ``mxExpectationNormal`` function, we refer to the expected covariance matrix by its name within the ``mxMatrix`` function that created the matrix, namely *expCov*.  However when we combine the arguments into the ``mxModel`` function, we use the name of the MxMatrix and MxMLObjective objects, respectively *expCovMat*, *expCovFun* and *ML*, as shown below.   
929
930 .. cssclass:: input
931 ..
932
933    OpenMx Code
934    
935    .. code-block:: r
936        
937          univSatModel3 <- mxModel("univSat3", obsCovData, expCovMat, expectCov, funML)
938          univSatFit3   <- mxRun(univSatModel3)
939          univSatSumm3  <- summary(univSatFit3)
940
941 Note that the estimates for the free parameters and the goodness-of-fit statistics are exactly the same for the matrix method as they were for the path method.
942
943     ..  code-block:: r
944     
945         > univSatSumm3
946         data:
947         $univSat3.data
948         $univSat3.data$cov
949                  X
950         X 1.062112
951     
952     
953         free parameters:
954           name matrix row col Estimate  Std.Error lbound ubound
955         1 <NA> expCov   X   X 1.062112 0.04752287              
956     
957         observed statistics:  1 
958         estimated parameters:  1 
959         degrees of freedom:  0 
960         -2 log likelihood:  1059.199 
961         saturated -2 log likelihood:  1059.199 
962         number of observations:  1000 
963         chi-square:  0 
964         p:  1 
965         Information Criteria: 
966             df Penalty Parameters Penalty Sample-Size Adjusted
967         AIC          0           2.000000                   NA
968         BIC          0           6.907755             3.731699
969         
970 We can also obtain the values of the likelihood by accessing the fitted object with the default name for the fit function, here *univSatFit4$fitfunction*.  Note the the expectation part of the fitted object has not changed.
971         
972 ..  code-block:: r        
973         
974         > univSatFit3$expectation
975         MxExpectationNormal 'expectation' 
976         $covariance : 'expCov' 
977         $means : NA 
978         $dims : 'X' 
979         $thresholds : NA 
980         $threshnames : 'X' 
981         
982         > univSatFit3$fitfunction
983         MxFitFunctionML 'fitfunction' 
984         $vector : FALSE 
985                  [,1]
986         [1,] 1059.199
987         attr(,"expCov")
988                  [,1]
989         [1,] 1.062112
990         attr(,"expMean")
991         <0 x 0 matrix>
992         attr(,"gradients")
993         <0 x 0 matrix>
994         attr(,"SaturatedLikelihood")
995         [1] 1059.199
996         attr(,"IndependenceLikelihood")
997         [1] 1059.199
998         
999     
1000 Covariances + Means
1001 ...................
1002
1003 A means vector can also be added to the observed data as the fourth argument of the ``mxData`` command.  When means are requested to be modeled, a second ``mxMatrix`` command is required to specify the vector of expected means. In this case a matrix of ``type`` ='Full', with one row and one column, is assigned ``free`` =T with start value zero, and the name *expMean*.  The object is saved as *expMeanMat*.  
1004
1005 .. cssclass:: input
1006 ..
1007
1008    OpenMx Code
1009    
1010    .. code-block:: r
1011        
1012          expMeanMat <- mxMatrix( type="Full", nrow=1, ncol=1, 
1013                                  free=TRUE, values=0, name="expMean" )
1014     
1015 When we inspect this MxMatrix object, note that it looks rather similar to the *expCovMat* object, except for the name and type and start value.  Its estimate depends entirely on which argument of the expectation function it is supposed to represent.  As soon as we move to an example with more than one variable, the difference becomes more obvious as the expected means will be a vector while the expected covariance matrix will always be a symmetric matrix.
1016
1017     ..  code-block:: r
1018
1019        > exMeanMat
1020         SymmMatrix 'expMean' 
1021
1022         $labels: No labels assigned.
1023
1024         $values
1025              [,1]
1026         [1,]    0
1027
1028         $free
1029              [,1]
1030         [1,] TRUE
1031
1032         $lbound: No lower bounds assigned.
1033
1034         $ubound: No upper bounds assigned.
1035
1036 The second change is adding an additional argument ``means`` to the ``mxExpectationNormal`` function for the expected mean, here *expMean*.
1037
1038 .. cssclass:: input
1039 ..
1040
1041    OpenMx Code
1042    
1043    .. code-block:: r
1044        
1045          expextCovMean <- mxExpectationNormal( covariance="expCov", means="expMean", 
1046                                                dimnames=selVars )
1047
1048 We now create a new model based on the old one, give it a new name, read in the MxData object with covariance and mean, add the MxMatrix object for the means and change the expectation function to the one created above.
1049
1050 .. cssclass:: input
1051 ..
1052
1053    OpenMx Code
1054    
1055    .. code-block:: r
1056        
1057          univSatModel3M <- mxModel(univSatModel3, name="univSat3M", obsCovMeanData, 
1058                                    expMeanMat, expextCovMean, funML )
1059          univSatFit3M   <- mxRun(univSatModel3M)
1060          univSatSumm3M  <- summary(univSatFit3M)
1061
1062 You can verify that the only changes to the output are the addition of the means to the data and estimates, resulting in two observed statistics and two estimated parameters rather than one.  As a result the values AIC and BIC criteria have changed although the value for the likelihood is exactly the same as before.
1063
1064
1065 Raw Data 
1066 ++++++++
1067
1068 Finally, if we want to use the matrix specification with raw data, no changes are needed to the matrices for the means and covariances, or to the expectation which combines the two.  Instead of summary statistics, we now fit the model to the raw data, saved in the MxData object *obsRawData*.  The fit function is still the same ``mxFitFunctionML()`` but now uses FIML (Full Information Maximum Likelihood), appropriate for raw data to evaluate the likelihood of the data .
1069
1070 The MxModel object for the saturated model applied to raw data has a name *univSat4*, a MxData object *obsRawData*, a MxMatrix object for the expected covariance matrix *expCovMat*, a MxMatrix object for the expected means vector *expMeanMat*,  a mxExpectationNormal object *expCovMeanFun*, and a mxFitFunction object *ML*.
1071
1072 .. cssclass:: input
1073 ..
1074
1075    OpenMx Code
1076    
1077    .. code-block:: r
1078        
1079          univSatModel4 <- mxModel("univSat4", obsRawData, 
1080                                   expCovMat, expMeanMat, expectCovMean, funML )
1081          univSatFit4   <- mxRun(univSatModel4)
1082          univSatSumm4  <- summary(univSatFit4)
1083
1084     ..  code-block:: r
1085     
1086         > univSatSumm4
1087         data:
1088         $univSat4.data
1089                X           
1090          Min.   :-3.32078  
1091          1st Qu.:-0.64970  
1092          Median : 0.03690  
1093          Mean   : 0.01681  
1094          3rd Qu.: 0.70959  
1095          Max.   : 3.30415  
1096     
1097         free parameters:
1098           name  matrix row col   Estimate  Std.Error lbound ubound
1099         1 <NA>  expCov   X   X 1.06104925 0.04745032              
1100         2 <NA> expMean   1   X 0.01680499 0.03257294              
1101     
1102         observed statistics:  1000 
1103         estimated parameters:  2 
1104         degrees of freedom:  998 
1105         -2 log likelihood:  2897.135 
1106         saturated -2 log likelihood:  NA 
1107         number of observations:  1000 
1108         chi-square:  NA 
1109         p:  NA 
1110         Information Criteria: 
1111             df Penalty Parameters Penalty Sample-Size Adjusted
1112         AIC   901.1355           2901.135                   NA
1113         BIC -3996.8043           2910.951             2904.599
1114     
1115 Note that the output generated for the paths and matrices specification are again completely equivalent, regardless of whether the model was fitted to summary statistics or raw data.  In each of the four versions of the model fitted to the same data, the data objects were generated from the continuous data.  Similar models can be fit to categorical data, with one or more thresholds delineating the proportion of individual in each of the two or more categories, based on the assumption of an underlying (multi)normal probability density function.
1116
1117 Threshold Model
1118 +++++++++++++++
1119
1120 Binary Data
1121 ...........
1122
1123 We will show below - only for the version using the matrix method to build a model to be fitted to the raw data - which changes are required when the input data is categorical.  We'll start with a binary example, followed by an ordinal one.
1124
1125 First, we read in the binary data, *obsRawDataBin* created earlier.  Then we turn the symmetric predicted covariance matrix into a standardized matrix with the variance of categorical variables (on the diagonal) fixed to one.  To estimate the thresholds, we need to fix the mean to zero, by changing the ``type`` argument to 'Zero'.  The one new object that is required is a matrix for the thresholds which will be estimated.  For binary data, the threshold matrix is similar to the means matrix before.
1126
1127 .. cssclass:: input
1128 ..
1129
1130    OpenMx Code
1131    
1132    .. code-block:: r
1133        
1134         expCovMatBin  <- mxMatrix( type="Stand", nrow=1, ncol=1, 
1135                                    free=TRUE, values=.5, name="expCov" )
1136         expMeanMatBin <- mxMatrix( type="Zero", nrow=1, ncol=1, name="expMean" )
1137         expThreMatBin <- mxMatrix( type="Full", nrow=1, ncol=1, 
1138                                    free=TRUE, values=0, name="expThre" )
1139         
1140 Let's inspect the latter matrix.
1141
1142     ..  code-block:: r
1143
1144         > expThreMatBin
1145         FullMatrix 'expThre' 
1146
1147         $labels: No labels assigned.
1148
1149         $values
1150              [,1]
1151         [1,]    0
1152
1153         $free
1154              [,1]
1155         [1,] TRUE
1156
1157         $lbound: No lower bounds assigned.
1158
1159         $ubound: No upper bounds assigned.
1160         
1161 The final change is adding an additional ``threshold`` argument to the ``mxExpectationNormal`` function for the expected threshold, here *expThreMatBin*.
1162
1163 .. cssclass:: input
1164 ..
1165
1166    OpenMx Code
1167    
1168    .. code-block:: r
1169        
1170          expectBin <- mxExpectationNormal( covariance="expCov", means="expMean", 
1171                                            threshold="expThre", dimnames=selVars )
1172
1173 We then include all these objects into a model *univSat5* and fit it to the data.
1174
1175 .. cssclass:: input
1176 ..
1177
1178    OpenMx Code
1179    
1180    .. code-block:: r
1181        
1182          univSatModel5 <- mxModel("univSat5", obsRawDataBin, 
1183                                   expCovMatBin, expMeanMatBin, expThreMatBin, expectBin, funML )
1184          univSatFit5   <- mxRun(univSatModel5)
1185          univSatSumm5  <- summary(univSatFit5)
1186          
1187 The summary of the univariate model fitted to binary data includes a summary of the data.  Given binary data have no variance, it is fixed to one while the threshold is estimated.
1188
1189     ..  code-block:: r
1190
1191          > univSatSumm5
1192          data:
1193          $univSat5.data
1194           X      
1195           0:680  
1196           1:320  
1197
1198          free parameters:
1199            name  matrix row col  Estimate  Std.Error lbound ubound
1200          1 <NA> expThre   1   X 0.4676989 0.04124951              
1201
1202          observed statistics:  1000 
1203          estimated parameters:  1 
1204          degrees of freedom:  999 
1205          -2 log likelihood:  1253.739 
1206          saturated -2 log likelihood:  NA 
1207          number of observations:  1000 
1208          chi-square:  NA 
1209          p:  NA 
1210          Information Criteria: 
1211              df Penalty Parameters Penalty Sample-Size Adjusted
1212          AIC  -744.2611           1255.739                   NA
1213          BIC -5647.1086           1260.647             1257.471
1214          CFI: NA 
1215          TLI: NA 
1216          RMSEA:  NA 
1217          timestamp: 2012-02-24 00:32:39 
1218          frontend time: 0.1296248 secs 
1219          backend time: 0.007578135 secs 
1220          independent submodels time: 5.102158e-05 secs 
1221          wall clock time: 0.137254 secs 
1222          cpu time: 0.137254 secs 
1223          openmx version number: 999.0.0-1661 
1224         
1225
1226
1227 Ordinal Data
1228 ............
1229
1230 Next, we will show how to adapt the model to analyze an ordinal variable.  As the number of thresholds depends on the variable, we specify it first, by creating a number of thresholds *nth* object.  The matrices for the expected covariance matrices and expected means are the same as in the binary case.  The matrix for the thresholds, however, now has as many rows as there are thresholds.  Furthermore, start values should be increasing. Here, we estimate the thresholds directly though.
1231
1232 .. cssclass:: input
1233 ..
1234
1235    OpenMx Code
1236    
1237    .. code-block:: r
1238        
1239         nth <- 4
1240         expCovMatOrd  <- mxMatrix( type="Stand", nrow=1, ncol=1, 
1241                                    free=TRUE, values=.5, name="expCov" )
1242         expMeanMatOrd <- mxMatrix( type="Zero", nrow=1, ncol=1, name="expMean" )
1243         expThreMatOrd <- mxMatrix( type="Full", nrow=nth, ncol=1, 
1244                                    free=TRUE, values=c(-1.5,-.5,.5,1.5), name="expThre" )
1245     
1246 Here we print the matrix of thresholds:
1247
1248     ..  code-block:: r
1249
1250         > expThreMatOrd
1251         FullMatrix 'expThre' 
1252
1253         $labels: No labels assigned.
1254
1255         $values
1256              [,1]
1257         [1,] -1.5
1258         [2,] -0.5
1259         [3,]  0.5
1260         [4,]  1.5
1261
1262         $free
1263              [,1]
1264         [1,] TRUE
1265         [2,] TRUE
1266         [3,] TRUE
1267         [4,] TRUE
1268
1269         $lbound: No lower bounds assigned.
1270
1271         $ubound: No upper bounds assigned.
1272         
1273 The remainder of the model statements is almost identical to those of the binary model, except for replacing 'Bin' with 'Ord'.
1274
1275 .. cssclass:: input
1276 ..
1277
1278    OpenMx Code
1279    
1280    .. code-block:: r
1281        
1282         expFunOrd     <- mxExpectationNormal( covariance="expCov", means="expMean", 
1283                                               threshold="expThre", dimnames=selVars )
1284         univSatModel6 <- mxModel("univSat6", obsRawDataOrd, 
1285                                  expCovMatOrd, expMeanMatOrd, expThreMatOrd, expectOrd, funML )
1286         univSatFit6   <- mxRun(univSatModel6)
1287         univSatSumm6  <- summary(univSatFit6)
1288
1289 Thresholds
1290 ............
1291
1292 An alternative approach to ensure that the thresholds are increasing can be enforced through multiplying the threshold matrix with a lower triangular matrix of 'Ones' and bounding all threshold increments except the first to be positive. The first threshold will be estimated as before.  The remaining thresholds are estimated as increments from the previous thresholds.  To generalize this, we specify a start value for the lower threshold ('svLTh') and for the increments ('svITh'), and then create a vector of start values to match the number of thresholds ('svTh').  Similarly, a vector of lower bounds is defined with all thresholds, except the first bounded to be positive ('lbTh').  These start values and lower bounds are read in to a MxMatrix object, of size *nth x 1*, similar to the threshold matrix in the previous example.  Then, we create a lower triangular matrix of ones which will be pre-multiplied with the threshold matrix to generate the expected threshold matrix *expThreMatOrd*.  The rest of the model is not changed, except that all the intermediate matrices, named *threG* and *inc* also have to be included in the MxModel object *univSatModel6I*.
1293
1294 .. cssclass:: input
1295 ..
1296
1297    OpenMx Code
1298    
1299    .. code-block:: r
1300    
1301         svLTh     <- -1.5                              # start value for first threshold
1302         svITh     <- 1                                 # start value for increments
1303         svTh      <- (c(svLTh,(rep(svITh,nth-1))))     # start value for thresholds
1304         lbTh      <- c(-3,(rep(0.001,nth-1)))          # lower bounds for thresholds
1305         
1306         threG          <- mxMatrix( type="Full", nrow=nth, ncol=1, 
1307                                     free=TRUE, values=svTh, lbound=lbTh, name="Thre" )
1308         inc            <- mxMatrix( type="Lower", nrow=nth, ncol=nth, 
1309                                     free=FALSE, values=1, name="Inc" )        
1310         expThreMatOrd  <- mxAlgebra( expression= Inc %*% Thre, name="expThre" )
1311         expectOrd      <- mxExpectationNormal( covariance="expCov", means="expMean", 
1312                                                threshold="expThre", dimnames=selVars )
1313         univSatModel6I <- mxModel("univSat6", obsRawDataOrd, 
1314                                   expCovMatOrd, expMeanMatOrd, 
1315                                   Inc, Thre, expThreMatOrd, expectOrd, funML )
1316         univSatFit6I   <- mxRun(univSatModel6I, unsafe=T)
1317         univSatSumm6I  <- summary(univSatFit6I)
1318         
1319
1320 Approaches 
1321 ----------
1322
1323 Rarely will we analyze a single variable.  As soon as a second variable is added, not only can we estimate both means and  variances, but also a covariance between the two variables, as shown in the following path diagram:
1324
1325 .. image:: graph/BivariateSaturatedModel.png
1326     :height: 1.0in
1327   
1328 The path diagram for our bivariate example includes two boxes for the observed variables 'X' and 'Y', each with a two-headed arrow for the variance of each of the variables.  We also estimate a covariance between the two variables with the two-headed arrow connecting the two boxes.  The optional means are represented as single-headed arrows from a triangle to the two boxes.
1329
1330 As raw data are now standard for data analysis, we will focus this example on fitting directly to the raw data.  We will present the example in both the path and the matrix specification, and furthermore show not only the piecewise style but also the stepwise and the classic style of writing OpenMx scripts.
1331
1332 Piecewise Style
1333 ^^^^^^^^^^^^^^^
1334
1335 Here we will illustrate the various approaches with the bivariate example.  For the piecewise approach, we'll show both the path specification and the matrix specification.  The other two approaches, stepwise and classic, will just be shown for the matrix example as specifying models using matrix algebra allows for greater flexibility and variety of models to be built.
1336
1337 Path Method
1338 ++++++++++++
1339
1340 In the path specification, we will use three ``mxPath`` commands to specify (i) the variance paths, (ii) the covariance path, and (iii) the mean paths.  We first specify the number of variables *nv* and which variables are selected for analysis *selVars*.
1341
1342
1343 .. cssclass:: input
1344 ..
1345
1346    OpenMx Code
1347    
1348    .. code-block:: r
1349        
1350         nv      <- 2
1351         selVars <- c('X','Y')
1352         
1353 We start with the two-headed paths for the variances and covariances.  The first one specifies two-headed arrows from **X** and **Y** to themselves - the ``to`` argument is omitted - to represent the variances.  This command now generates two free parameters, each with start value of 1 and lower bound of .01, but with a different label indicating that these are separate free parameters.  Note that we could test whether the variances are equal by specifying a model with the same label for the two variances and comparing it with the current model.  The second ``mxPath`` command specifies a two-headed arrow from **X** to **Y** for the covariance, which is also assigned 'free' and given a start value of .2 and a label.
1354
1355 .. cssclass:: input
1356 ..
1357
1358    OpenMx Code
1359    
1360    .. code-block:: r
1361        
1362          expVars <- mxPath( from=c("X", "Y"), arrows=2, 
1363                             free=TRUE, values=1, lbound=.01, labels=c("varX","varY") )
1364          expCovs <- mxPath( from="X", to="Y", arrows=2, 
1365                             free=TRUE, values=.2, lbound=.01, labels="covXY" )
1366
1367 The resulting MxPath objects 'expVars' and 'expCovs' are as follows:
1368
1369     ..  code-block:: r
1370     
1371        > mxPath( from=c("X", "Y"), arrows=2, 
1372                  free=TRUE, values=1, lbound=.01, labels=c("varX","varY") )
1373         mxPath 
1374         $from:  'X' and 'Y' 
1375         $to:  'X' and 'Y' 
1376         $arrows:  2 
1377         $values:  1 
1378         $free:  TRUE 
1379         $labels:  varX varY 
1380         $lbound:  0.01 
1381         $ubound:  NA 
1382         >  mxPath( from="X", to="Y", arrows=2, 
1383                    free=TRUE, values=.2, lbound=.01, labels="covXY" )
1384         mxPath 
1385         $from:  'X' 
1386         $to:  'Y' 
1387         $arrows:  2 
1388         $values:  0.2 
1389         $free:  TRUE 
1390         $labels:  covXY 
1391         $lbound:  0.01 
1392         $ubound:  NA 
1393         
1394 When observed means are included in addition to the observed covariance matrix, as is necessary when fitting to raw data, we add an ``mxPath`` command with single-headed arrows from ``one`` to the variables to represent the two means.
1395
1396 .. cssclass:: input
1397 ..
1398
1399    OpenMx Code
1400    
1401    .. code-block:: r
1402        
1403          expMeans <- mxPath( from="one", to=c("X", "Y"), arrows=1, 
1404                              free=TRUE, values=.01, labels=c("meanX","meanY") )
1405
1406 The "one" argument in the ``from`` argument is used exclusively for means objects, here called *expMeans*.
1407
1408     ..  code-block:: r
1409     
1410        > mxPath( from="one", to=c("X", "Y"), arrows=1, 
1411                  free=TRUE, values=.01, labels=c("meanX","meanY") )
1412         mxPath 
1413         $from:  'one' 
1414         $to:  'X' and 'Y' 
1415         $arrows:  1 
1416         $values:  0.01 
1417         $free:  TRUE 
1418         $labels:  meanX meanY 
1419         $lbound:  NA 
1420         $ubound:  NA 
1421         
1422
1423 To fit this bivariate model to the simulated data, we have to combine the data and model statements in a MxModel objects.
1424
1425 .. cssclass:: input
1426 ..
1427
1428    OpenMx Code
1429    
1430    .. code-block:: r
1431        
1432          bivSatModel1 <- mxModel("bivSat1", manifestVars=selVars, obsBivData, 
1433                                  expVars, expCovs, expMeans, type="RAM" )
1434          bivSatFit1   <- mxRun(bivSatModel1)
1435          bivSatSumm1  <- summary(bivSatFit1)
1436     
1437 As you can see below, the maximum likelihood (ML) estimates are very close to the summary statistics of the simulated data.
1438
1439     ..  code-block:: r
1440     
1441         > bivSatSumm1
1442         data:
1443         $bivSat1.data
1444                X                   Y            
1445          Min.   :-3.296159   Min.   :-2.942561  
1446          1st Qu.:-0.596177   1st Qu.:-0.633711  
1447          Median :-0.010538   Median :-0.004139  
1448          Mean   :-0.004884   Mean   : 0.032116  
1449          3rd Qu.: 0.598326   3rd Qu.: 0.739236  
1450          Max.   : 4.006771   Max.   : 4.173841  
1451
1452         free parameters:
1453            name matrix row col     Estimate  Std.Error Std.Estimate     Std.SE lbound ubound
1454         1  varX      S   X   X  0.993537344 0.04443221    1.0000000 0.04472123   0.01       
1455         2 covXY      S   X   Y  0.481348846 0.03513471    0.4806856 0.03508630   0.01       
1456         3  varY      S   Y   Y  1.009283953 0.04513849    1.0000000 0.04472328   0.01       
1457         4 meanX      M   1   X -0.004884421 0.03152067           NA         NA              
1458         5 meanY      M   1   Y  0.032116307 0.03177008           NA         NA              
1459
1460         observed statistics:  0 
1461         estimated parameters:  5 
1462         degrees of freedom:  -5 
1463         -2 log likelihood:  5415.772 
1464         saturated -2 log likelihood:  -2 
1465         number of observations:  1000 
1466         chi-square:  5417.772 
1467         p:  NaN 
1468         Information Criteria: 
1469              df Penalty Parameters Penalty Sample-Size Adjusted
1470         AIC:   5425.772           5425.772                   NA
1471         BIC:   5450.311           5450.311             5434.431
1472
1473
1474 Matrix Method
1475 ++++++++++++++
1476
1477 If we use matrices instead of paths to specify the bivariate model, we need to generate matrices to represent the expected covariance matrix and the means.  The ``mxMatrix`` command for the expected covariance matrix now specifies a **2x2** symmetric matrix with all elements free.  Start values have to be given only for the unique elements (diagonal elements plus upper or lower diagonal elements), in this case we provide a list with values of 1 for the variances and 0.5 for the covariance.
1478
1479 .. cssclass:: input
1480 ..
1481
1482    OpenMx Code
1483    
1484    .. code-block:: r
1485        
1486          expCovM <- mxMatrix( type="Symm", nrow=2, ncol=2, 
1487                               free=TRUE, values=c(1,0.5,1), 
1488                               labels=c('V1','Cov','V2'), name="expCov" )
1489     
1490 By specifying labels, we can tell that the two covariance elements, expCovM[1,2] and expCovM[2,1] are constrained to be equal - implied by the fact that the label is the same.  This of course automatically happens when you specify the matrix to be symmetric.
1491
1492     ..  code-block:: r
1493     
1494         > expCovM
1495         SymmMatrix 'expCov' 
1496     
1497         $labels
1498              [,1]  [,2] 
1499         [1,] "V1"  "Cov"
1500         [2,] "Cov" "V2" 
1501     
1502         $values
1503              [,1] [,2]
1504         [1,]  1.0  0.5
1505         [2,]  0.5  1.0
1506     
1507         $free
1508              [,1] [,2]
1509         [1,] TRUE TRUE
1510         [2,] TRUE TRUE
1511     
1512         $lbound: No lower bounds assigned.
1513     
1514         $ubound: No upper bounds assigned.
1515     
1516 When fitting to raw data, we also use a ``mxMatrix`` command to specify the expected means as **1x2** row vector with two free parameters, each given a 0 start value.
1517
1518 .. cssclass:: input
1519 ..
1520
1521    OpenMx Code
1522    
1523    .. code-block:: r
1524        
1525          expMeanM <- mxMatrix( type="Full", nrow=1, ncol=2, 
1526                                free=TRUE, values=c(0,0), labels=c('M1','M2'), name="expMean" )
1527
1528 Similarly to above, the elements in this matrix can also be given labels, although this is entirely optional for both matrices.  However, as soon as we want to change the model to e.g. test equality of means or variances, the most efficient way to do this is by using labels.  Given fitting alternative models to test hypotheses is very common, we highly recommend to use labels at all times.  Note that we truncated the output below as no bounds had been assigned.
1529
1530     ..  code-block:: r
1531     
1532        > expMeanM
1533         FullMatrix 'expMean' 
1534     
1535         $labels
1536              [,1] [,2]
1537         [1,] "M1" "M2"
1538     
1539         $values
1540              [,1] [,2]
1541         [1,]    0    0
1542     
1543         $free
1544              [,1] [,2]
1545         [1,] TRUE TRUE
1546     
1547 So far, we have specified the expected covariance matrix directly as a symmetric matrix.  However, this may cause optimization problems as the matrix could become not positive-definite which would prevent the likelihood to be evaluated.  To overcome this problem, we can use a Cholesky decomposition of the expected covariance matrix instead, by multiplying a lower triangular matrix with its transpose.  To obtain this, we use a ``mxMatrix`` command and define a **2x2** lower triangular matrix using ``type`` ="Lower", declare all elements free with 0.5 starting values.  We name this matrix "Chol" and save the object as *lowerTriM*.
1548
1549 .. cssclass:: input
1550 ..
1551
1552    OpenMx Code
1553    
1554    .. code-block:: r
1555        
1556          lowerTriM <- mxMatrix( type="Lower", nrow=2, ncol=2, 
1557                                 free=TRUE, values=.5, name="Chol" )
1558     
1559 Given we specified the matrix as lower triangular, the start values and free assignments are only applied to the elements on the diagonal and below the diagonal.
1560
1561     ..  code-block:: r
1562     
1563         > lowerTriM
1564         LowerMatrix 'Chol' 
1565     
1566         $labels: No labels assigned.
1567     
1568         $values
1569              [,1] [,2]
1570         [1,]  0.5  0.0
1571         [2,]  0.5  0.5
1572     
1573         $free
1574              [,1]  [,2]
1575         [1,] TRUE FALSE
1576         [2,] TRUE  TRUE
1577     
1578 We then use an ``mxAlgebra`` command to multiply this matrix with its transpose (R function ``t()``).  The ``mxAlgebra`` command is a very useful command to apply any operation or function to matrices.  It only has two arguments, the first for the ``expression`` you intend to generate, the second the name of the resulting matrix.  Note that although the matrix object for the lower triangular matrix was saved as 'lowerTwiM', the matrices in the ``expression`` are referred to by the name given to them within the MxMatrix object.  This is similar to referring to the names of the expected covariance matrices and means when they are needed in the arguments of the ``mxExpectationNormal`` function.
1579
1580 .. cssclass:: input
1581 ..
1582
1583    OpenMx Code
1584    
1585    .. code-block:: r
1586        
1587          expCovMA <- mxAlgebra( expression=Chol %*% t(Chol), name="expCov" )
1588     
1589 So far, we've only specified the algebra, but not computed it yet as shown when we look at the *expCovMA* object.  We need to combine all elements in an ``mxModel`` prior to ``mxRun`` ning the model to compute the algebra.
1590
1591     ..  code-block:: r
1592     
1593         > expCovMA
1594         mxAlgebra 'expCov' 
1595         $formula:  Chol %*% t(Chol) 
1596         $result: (not yet computed) <0 x 0 matrix>
1597         dimnames: NULL
1598     
1599 Given we used the same names for the resulting matrices for the expected covariances and means as in the univariate example, the ``mxExpectationNormal`` command looks identical.  Note that we have redefined *selVars* when starting the bivariate examples.  When you use the piecewise style and you're running more than one job, make sure you're not accidentally using an object from a previous job, especially if you've made an error in a newly specified object with the same name.
1600
1601 .. cssclass:: input
1602 ..
1603
1604    OpenMx Code
1605    
1606    .. code-block:: r
1607        
1608          expectBiv <- mxExpectationNormal( covariance="expCov", means="expMean", 
1609                                            dimnames=selVars )
1610  
1611 Combining these two ``mxMatrix`` and the ``mxAlgebra`` objects with the raw data, specified in the ``mxData`` object 'obsBivData' created earlier and the ``mxExpectationNormal`` command with the appropriate arguments is all that is needed to fit a saturated bivariate model.
1612
1613 .. cssclass:: input
1614 ..
1615
1616    OpenMx Code
1617    
1618    .. code-block:: r
1619        
1620           bivSatModel2 <- mxModel("bivSat2", obsBivData, lowerTriM, 
1621                                   expCovMA, expMeanM, expectBiv, funML )
1622           bivSatFit2   <- mxRun(bivSatModel2)
1623           bivSatSumm2  <- summary(bivSatFit2)
1624     
1625 The goodness-of-fit statistics in the output from the path and matrix specification appear identical.  However, in the latter model, we do not model the variances and covariance directly but parameterize them using a Cholesky decomposition. 
1626
1627     ..  code-block:: r
1628     
1629          > bivSatSumm2
1630          data:
1631          $bivSat2.data
1632                 X                   Y            
1633           Min.   :-3.296159   Min.   :-2.942561  
1634           1st Qu.:-0.596177   1st Qu.:-0.633711  
1635           Median :-0.010538   Median :-0.004139  
1636           Mean   :-0.004884   Mean   : 0.032116  
1637           3rd Qu.: 0.598326   3rd Qu.: 0.739236  
1638           Max.   : 4.006771   Max.   : 4.173841  
1639
1640          free parameters:
1641                         name  matrix row col     Estimate  Std.Error lbound ubound
1642          1 bivSat2.Chol[1,1]    Chol   1   1  0.996763911 0.02228823              
1643          2 bivSat2.Chol[2,1]    Chol   2   1  0.482912544 0.02987720              
1644          3 bivSat2.Chol[2,2]    Chol   2   2  0.880954091 0.01969863              
1645          4                M1 expMean   1   X -0.004883967 0.03151918              
1646          5                M2 expMean   1   Y  0.032116277 0.03176869              
1647
1648          observed statistics:  2000 
1649          estimated parameters:  5 
1650          degrees of freedom:  1995 
1651          -2 log likelihood:  5415.772 
1652          saturated -2 log likelihood:  -2 
1653          number of observations:  1000 
1654          chi-square:  5417.772 
1655          p:  2.595415e-313 
1656          Information Criteria: 
1657               df Penalty Parameters Penalty Sample-Size Adjusted
1658          AIC:   1425.772           5425.772                   NA
1659          BIC:  -8365.200           5450.311             5434.431
1660          
1661          
1662 We can obtain the predicted variances and covariances by printing the *expCov* matrix which can be done with the ``mxEval`` command - either by recalculating or by just printing the calculated algebra, or by grabbing the predicted covariance matrix from the fitted object *bivSatFit2*
1663
1664 .. cssclass:: input
1665 ..
1666
1667    OpenMx Code
1668    
1669    .. code-block:: r
1670        
1671           mxEval(Chol %*% t(Chol), bivSatFit2 )
1672           mxEval(expCov, bivSatFit2 )
1673           bivSatFit2$expCov$result
1674
1675 So far, we have presented the bivariate model (path or matrix method) using the piecewise approach.  As a result, we end up with a series of OpenMx objects, each of which we can check for syntax correctness.  As such, this is a great way to build new models.  An alternative approach is to start the mxModel with one argument, and then add another argument step by step.  We will show the various steps of building the bivariate model with matrices (and algebras).
1676
1677 Here, we simply repeat all the lines that make up the model.
1678
1679 .. cssclass:: input
1680 ..
1681
1682    OpenMx Code
1683    
1684    .. code-block:: r
1685        
1686           obsBivData   <- mxData( observed=bivData, type="raw" )
1687           expMeanM     <- mxMatrix( type="Full", nrow=1, ncol=2, 
1688                                     free=TRUE, values=0, labels=c('M1','M2'), name="expMean" )
1689           lowerTriM    <- mxMatrix( type="Lower", nrow=2, ncol=2, 
1690                                     free=TRUE, values=.5, name="Chol" )
1691           expCovMA     <- mxAlgebra( expression=Chol %*% t(Chol), name="expCov" )
1692           expectBiv    <- mxExpectationNormal( covariance="expCov", means="expMean", 
1693                                                dimnames=selVars )
1694           funML        <- mxFitFunctionML()
1695           bivSatModel2 <- mxModel("bivSat2", obsBivData, lowerTriM, 
1696                                   expCovMA, expMeanM, expectBiv, funML )
1697           bivSatFit2   <- mxRun(bivSatModel2)
1698           bivSatSumm2  <- summary(bivSatFit2)
1699     
1700
1701 Stepwise Style
1702 ^^^^^^^^^^^^^^
1703
1704 Looking back at the MxModel we just built (*bivSatModel2*), the first argument after the name (in quotes) was the MxData object.  Let's now build a new model that just has a new name and the data object to start, specified from scratch - assuming we had not built the object before.  Note we need to close both the ``mxData`` command which resides within the ``mxModel`` command.  We can execute it in R, to check for syntax errors.
1705
1706 .. cssclass:: input
1707 ..
1708
1709    OpenMx Code
1710    
1711    .. code-block:: r
1712        
1713           bivSatModel3 <- mxModel("bivSat3", mxData( observed=bivData, type="raw" ) )
1714
1715 If it checks out to be syntactically correct, we can add another argument, i.e. the ``mxMatrix`` command to define the lower triangular matrix.  Given we now want to build upon the previous model, we use that as the first argument (without quotes).  As the previous model already has a name argument we don't need to include that and can go straight to the new argument.
1716
1717 .. cssclass:: input
1718 ..
1719
1720    OpenMx Code
1721    
1722    .. code-block:: r
1723        
1724           bivSatModel3 <- mxModel(bivSatModel3, 
1725                                   mxMatrix( type="Lower", nrow=2, ncol=2, 
1726                                             free=TRUE, values=.5, name="Chol" ) )
1727
1728 Note that we used the same name for the MxModel object to be generated, thus it will overwrite the previous one.  One might choose a new name if uncertain about the syntax of the new model, to avoid having to rerun the previous step to correct the original model.  As everyone checks out OK after step two, let us add another argument, this time the ``mxAlgebra`` object for the expected covariance matrix. 
1729
1730 .. cssclass:: input
1731 ..
1732
1733    OpenMx Code
1734    
1735    .. code-block:: r
1736        
1737           bivSatModel3 <- mxModel(bivSatModel3, 
1738                                   mxAlgebra( expression=Chol %*% t(Chol), name="expCov" ) )
1739      
1740 As everything still appears OK, we continue to add arguments.  It's not necessary to do them one at the time, but if you're just learning to build a model, it might be the safest bet.  Next, we add the ``mxMatrix`` command for the expected means.
1741
1742 .. cssclass:: input
1743 ..
1744
1745    OpenMx Code
1746    
1747    .. code-block:: r
1748        
1749           bivSatModel3 <- mxModel(bivSatModel3, 
1750                                   mxMatrix( type="Full", nrow=1, ncol=2, 
1751                                   free=TRUE, values=0, name="expMean" ) )
1752      
1753
1754 The only argument left to add is the ``mxExpectationNormal`` and the ``mxFitFunctionML()`` to specify what function to use on test the fit between covariances and means predicted by the observed data and those expected by the model.
1755
1756 .. cssclass:: input
1757 ..
1758
1759    OpenMx Code
1760    
1761    .. code-block:: r
1762        
1763           bivSatModel3 <- mxModel(bivSatModel3, 
1764                                   mxExpectationNormal( covariance="expCov", means="expMean", 
1765                                   dimnames=selVars ) )
1766           bivSatModel3 <- mxModel(bivSatModel3, mxFitFunctionML() )
1767
1768 Now that we have the complete model built, we can evaluate it using the ``mxRun`` command.
1769
1770 .. cssclass:: input
1771 ..
1772
1773    OpenMx Code
1774    
1775    .. code-block:: r
1776        
1777           bivSatFit3  <- mxRun(bivSatModel3)
1778           bivSatSumm3 <- summary(bivSatFit3)
1779
1780 You can verify for yourself that the results (goodness-of-fit statistics and parameter estimates) are entirely equivalent between this stepwise approach and the piecewise approach.
1781
1782 We here combine all the separate lines to see the full picture.
1783
1784 .. cssclass:: input
1785 ..
1786
1787    OpenMx Code
1788    
1789    .. code-block:: r
1790        
1791           bivSatModel3 <- mxModel("bivSat3", mxData( observed=bivData, type="raw" ) )
1792           bivSatModel3 <- mxModel(bivSatModel3, 
1793                                   mxMatrix( type="Lower", nrow=2, ncol=2, 
1794                                   free=TRUE, values=.5, name="Chol" ) )
1795           bivSatModel3 <- mxModel(bivSatModel3, 
1796                                   mxAlgebra( expression=Chol %*% t(Chol), name="expCov" ) )
1797           bivSatModel3 <- mxModel(bivSatModel3, 
1798                                   mxMatrix( type="Full", nrow=1, ncol=2, 
1799                                   free=TRUE, values=0, name="expMean" )
1800           bivSatModel3 <- mxModel(bivSatModel3, 
1801                                   mxExpectationNormal( covariance="expCov", means="expMean", 
1802                                   dimnames=selVars ) )
1803           bivSatFit3   <- mxRun(bivSatModel3)
1804           bivSatSumm3  <- summary(bivSatFit3)
1805
1806
1807 Classic Style
1808 ^^^^^^^^^^^^^
1809
1810 If you are fairly confident that you can specify each of the arguments of the model syntactically correct, there is no need to build the objects one by one and combine them, as we did using the piecewise approach, or to build models step by step by adding one argument at a time, as we did in the stepwise approach.  Instead, we can generate the complete syntax at once.  As a result, this will be the most concise way to write and run models.  The disadvantage, however, is that if you make changes to the model, and they include a syntax error, it is less evident to find the error.  The advantage, on the other hand, is that some models do not require any changes, but maybe you just want to apply them to different data sets in which case this approach works fine.
1811
1812 Here, we present the complete bivariate saturated model, with each argument printed on a different line for clarity of presentation.  To not overwrite the previous objects, we'll start with a new name for the MxModel object.  Remember that  arguments have to be separated by comma's, and that we need a double bracket after the last argument to close both that argument and the full model.
1813
1814 .. cssclass:: input
1815 ..
1816
1817    OpenMx Code
1818    
1819    .. code-block:: r
1820        
1821           bivSatModel4 <- mxModel("bivSat4", 
1822                                    mxData( observed=bivData, type="raw" ), 
1823                                    mxMatrix( type="Lower", nrow=2, ncol=2, 
1824                                              free=TRUE, values=.5, name="Chol" ),
1825                                    mxAlgebra( expression=Chol %*% t(Chol), name="expCov" ),
1826                                    mxMatrix( type="Full", nrow=1, ncol=2, 
1827                                              free=TRUE, values=c(0,0), name="expMean" ),
1828                                    mxExpectationNormal( covariance="expCov", means="expMean", 
1829                                                         dimnames=selVars )
1830                                    mxFitFunctionML() )
1831            bivSatFit4 <- mxRun(bivSatModel4)
1832            bivSatSumm4 <- summary(bivSatFit4)
1833
1834
1835 Again, as you might expect by now, the output of this model run will be identical to that of both thez piecewise and the stepwise approach.  Given their equivalence, it is really up to the OpenMx user to decide which method (path or matrix) and which approach (piecewise, stepwise or classic) is preferred.  It is also not necessary to pick just one of these approaches, as they can be 'mixed and matched'.  For didactic purposes, we recommend the piecewise approach, which we will use in the majority of this documentation.  We will, however, provide some parallel classic scripts.  Furthermore, given some people do better with path diagrams and others with matrix algebra, we present some models both ways, in so far that this is doable.
1836
1837 The following sections describe OpenMx examples in detail beginning with regression, factor analysis, time series analysis, multiple group models, including twin models, and analysis using definition variables.  Each is presented in both path and matrix styles and where relevant, contrasting data input from covariance matrices versus raw data input are also illustrated.  Additional examples will be added as they are implemented in OpenMx.
1838
1839
1840