Expectation-Fit revision to multigroup User's Guide pages.
[openmx:openmx.git] / docs / source / MultipleGroups_Matrix.rst
1 .. _multiplegroups-matrix-specification:
2
3 Multiple Groups, Matrix Specification
4 =====================================
5
6 An important aspect of structural equation modeling is the use of multiple groups to compare means and covariances structures between any two (or more) data groups, for example males and females, different ethnic groups, ages etc.  Other examples include groups which have different expected covariances matrices as a function of parameters in the model, and need to be evaluated together for the parameters to be identified.
7
8 The example includes the heterogeneity model as well as its submodel, the homogeneity model and is available in the following file:
9
10 * http://openmx.psyc.virginia.edu/svn/trunk/demo/BivariateHeterogeneity_MatrixRaw.R
11
12 A parallel version of this example, using paths specification of models rather than matrices, can be found here:
13
14 * http://openmx.psyc.virginia.edu/svn/trunk/demo/BivariateHeterogeneity_PathRaw.R
15
16
17 Heterogeneity Model
18 -------------------
19
20 We will start with a basic example here, building on modeling means and variances in a saturated model.  Assume we have two groups and we want to test whether they have the same mean and covariance structure.  
21
22 The path diagram of the heterogeneity model for a set of variables :math:`x` and :math:`y` are given below.
23
24 .. math::
25 ..   :nowrap:
26    
27 ..   \begin{eqnarray*} 
28 ..   x = \mu_{x1} + \sigma_{x1}
29 ..   \end{eqnarray*}
30
31 .. image:: graph/HeterogeneityModel.png
32     :height: 2.5in
33
34 Data
35 ^^^^
36
37 For this example we simulated two datasets (``xy1`` and ``xy2``) each with zero means and unit variances, one with a correlation of 0.5, and the other with a correlation of 0.4 with 1000 subjects each.  We use the ``mvrnorm`` function in the ``MASS`` package, which takes three arguments: ``Sample Size``, ``Means``, ``Covariance Matrix``).  We check the means and covariance matrix in R and provide ``dimnames`` for the dataframe.  See attached R code for simulation and data summary.
38
39 .. code-block:: r
40
41     #Simulate Data
42     require(MASS)
43     #group 1
44     set.seed(200)
45     rs=0.5
46     xy1 <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
47     set.seed(200)
48     #group 2
49     rs=0.4
50     xy2 <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
51
52     #Print Descriptive Statistics
53     selVars <- c('X','Y')
54     summary(xy1)
55     cov(xy1)
56     dimnames(xy1) <- list(NULL, selVars)
57     summary(xy2)
58     cov(xy2)
59     dimnames(xy2) <- list(NULL, selVars)
60     
61     
62 Model Specification
63 ^^^^^^^^^^^^^^^^^^^
64
65 As before, we include the OpenMx package using a ``require`` statement.
66 We first fit a heterogeneity model, allowing differences in both the mean and covariance structure of the two groups.  As we are interested whether the two structures can be equated, we have to specify the models for the two groups, named ``group1`` and ``group2`` within another model, named ``bivHet``.  The structure of the job thus look as follows, with two ``mxModel`` commands as arguments of another ``mxModel`` command.  Note that ``mxModel`` commands are unlimited in the number of arguments.
67
68 .. code-block:: r
69
70     require(OpenMx)
71
72     bivHetModel <- mxModel("bivHet",
73          mxModel("group1"),
74          mxModel("group2"),
75          mxFitFunctionMultigroup(c("group1.fitfunction", "group2.fitfunction"))
76     )
77      
78 For each of the groups, we fit a saturated model, using a Cholesky decomposition to generate the expected covariance matrix and a row vector for the expected means.  Note that we have specified different labels for all the free elements, in the two ``mxModel`` statements.  For more details, see example 1.
79
80 .. code-block:: r
81
82     #Fit Heterogeneity Model
83     bivHetModel <- mxModel("bivHet",
84         mxModel("group1",
85             mxMatrix(
86                 type="Lower", 
87                 nrow=2, 
88                 ncol=2, 
89                 free=T, 
90                 values=.5,
91                 labels=c("Ch11", "Ch21", "Ch31"),
92                 name="Chol1"
93             ), 
94             mxAlgebra(
95                 Chol1 %*% t(Chol1), 
96                 name="EC1" 
97             ), 
98             mxMatrix(
99                 type="Full", 
100                 nrow=1, 
101                 ncol=2, 
102                 free=T, 
103                 values=c(0,0), 
104                 labels=c("mX1", "mY1"), 
105                 name="EM1"
106             ), 
107             mxData(
108                 xy1, 
109                 type="raw"
110             ), 
111             mxExpectationNormal(
112                 covariance="EC1", 
113                 means="EM1",
114                 dimnames=selVars
115             ),
116             mxFitFunctionML()
117         ),
118         mxModel("group2",
119             mxMatrix(
120                 type="Lower", 
121                 nrow=2, 
122                 ncol=2, 
123                 free=T, 
124                 values=.5,
125                 labels=c("Ch12", "Ch22", "Ch32"),
126                 name="Chol2"
127             ), 
128             mxAlgebra(
129                 Chol2 %*% t(Chol2), 
130                 name="EC2"
131             ), 
132             mxMatrix(
133                 type="Full", 
134                 nrow=1, 
135                 ncol=2, 
136                 free=T, 
137                 values=c(0,0), 
138                 labels=c("mX2", "mY2"), 
139                 name="EM2"
140             ), 
141             mxData(
142                 xy2, 
143                 type="raw"
144             ), 
145             mxExpectationNormal(
146                 covariance="EC2", 
147                 means="EM2",
148                 dimnames=selVars
149             ),
150             mxFitFunctionML()
151         ),
152
153
154 We estimate five parameters (two means, two variances, one covariance) per group for a total of 10 free parameters.  We cut the ``Labels matrix:`` parts from the output generated with ``bivHetModel$group1@matrices`` and ``bivHetModel$group2@matrices``::
155
156     in group1
157         $S
158                 X      Y     
159         X  "Ch11"     NA
160         Y  "Ch21"  "Ch22" 
161
162         $M
163                 X      Y    
164         [1,] "mX1" "mY1"
165
166     in group2
167         $S
168                 X      Y     
169         X  "Ch12"     NA
170         Y  "Ch22" "Ch32" 
171
172         $M
173                 X      Y    
174         [1,] "mX2" "mY2"
175
176 To evaluate both models together, we use an ``mxFitFunctionMultigroup`` command that adds up the values of the fit functions of the two groups.
177
178 .. code-block:: r
179
180         mxFitFunctionMultigroup(c("group1.fitfunction", "group2.fitfunction"))
181     )
182
183 Model Fitting
184 ^^^^^^^^^^^^^
185
186 The ``mxRun`` command is required to actually evaluate the model.  Note that we have adopted the following notation of the objects.  The result of the ``mxModel`` command ends in "Model"; the result of the ``mxRun`` command ends in "Fit".  Of course, these are just suggested naming conventions.
187
188 .. code-block:: r
189
190     bivHetFit <- mxRun(bivHetModel)
191
192 A variety of output can be printed.  We chose here to print the expected means and covariance matrices for the two groups and the likelihood of data given the model.  The ``mxEval`` command takes any R expression, followed by the fitted model name.  Given that the model ``bivHetFit`` included two models (group1 and group2), we need to use the two level names, i.e. ``group1.EM1`` to refer to the objects in the correct model.
193
194 .. code-block:: r
195
196     EM1Het <- mxEval(group1.EM1, bivHetFit)
197     EM2Het <- mxEval(group2.EM2, bivHetFit)
198     EC1Het <- mxEval(group1.EC1, bivHetFit)
199     EC2Het <- mxEval(group2.EC2, bivHetFit)
200     LLHet <- summary(bivHetFit)$Minus2LogLikelihood
201
202
203 Homogeneity Model: a Submodel
204 -----------------------------
205
206 Next, we fit a model in which the mean and covariance structure of the two groups are equated to one another, to test whether there are significant differences between the groups.  Rather than having to specify the entire model again, we copy the previous model ``bivHetModel`` into a new model ``bivHomModel`` to represent homogeneous structures.
207
208 .. code-block:: r
209
210     #Fit Homogeneity Model
211     bivHomModel <- bivHetModel
212
213 As elements in matrices can be equated by assigning the same label, we now have to equate the labels of the free parameters in group 1 to the labels of the corresponding elements in group 2.  This can be done by referring to the relevant matrices using the ``ModelName$MatrixName`` syntax, followed by ``@labels``.  Note that in the same way, one can refer to other arguments of the objects in the model.  Here we assign the labels from group1 to the labels of group2, separately for the Cholesky matrices used for the expected covariance matrices and for the expected means vectors.
214
215 .. code-block:: r
216
217     bivHomModel$group2.Chol2@labels <- bivHomModel$group1.Chol1@labels
218     bivHomModel$group2.EM2@labels <- bivHomModel$group1.EM1@labels
219
220 The specification for the submodel is reflected in the names of the labels which are now equal for the corresponding elements of the mean and covariance matrices, as below::
221
222     in group1
223         $S
224                 X      Y     
225         X  "Ch11"     NA
226         Y  "Ch21" "CH31" 
227
228         $M
229                 X      Y    
230         [1,] "mX1" "mY1"
231     
232     in group2
233         $S
234                 X      Y     
235         X  "Ch11"     NA
236         Y  "Ch21" "Ch31" 
237
238         $M
239                 X      Y     
240         [1,] "mX1" "mY1"
241
242 We can produce similar output for the submodel, i.e. expected means and covariances and likelihood, the only difference in the code being the model name.  Note that as a result of equating the labels, the expected means and covariances of the two groups should be the same.
243
244 .. code-block:: r
245
246     bivHomFit <- mxRun(bivHomModel)
247         EM1Hom <- mxEval(group1.EM1, bivHomFit)
248         EM2Hom <- mxEval(group2.EM2, bivHomFit)
249         EC1Hom <- mxEval(group1.EC1, bivHomFit)
250         EC2Hom <- mxEval(group2.EC2, bivHomFit)
251         LLHom <- summary(bivHomFit)$Minus2LogLikelihood
252
253 Finally, to evaluate which model fits the data best, we generate a likelihood ratio test as the difference between -2 times the log-likelihood of the homogeneity model and -2 times the log-likelihood of the heterogeneity model.  This statistic is asymptotically distributed as a Chi-square, which can be interpreted with the difference in degrees of freedom of the two models.
254
255 .. code-block:: r
256
257         Chi <- LLHom-LLHet
258         LRT <- rbind(LLHet,LLHom,Chi)
259         LRT
260
261 These models may also be specified using paths instead of matrices. See :ref:`multiplegroups-path-specification` for path specification of these models.