Note: This is just a copy of my course notes from last year. Given that snow days prevented covering the topic in more depth, I opted to keep my notes unchanged, rather than trying to compress the material, which would be tricky.

Note that this Rmd file depends on having Mplus installed. You can obtain a fully functioning student version for $350 here: https://www.statmodel.com/orderonline/categories.php?category=Mplus-Software/Student-Pricing.

Or you can use the demo version for free here: https://www.statmodel.com/demo.shtml

The demo version has the following constraints:

  1. Maximum number of dependent variables: 6
  2. Maximum number of independent variables: 2
  3. Maximum number of between variables in two-level analysis: 2
  4. Maximum number of continuous latent variables in time series analysis: 2

At this time, Yves Rosseel, the main developer of lavaan, has a prototype of multilevel SEM working for the package, but this has not been released to the general public. I suspect that this will be released in the next 6 months and should provide the functionality to run all of the examples in this RMarkdown file. You can see some of the basic setup here: http://users.ugent.be/~yrosseel/lavaan/tubingen2017/

In the meantime, Mplus is probably the most user-friendly program for multilevel SEM, though there is similar functionality in EQS and LISREL. In addition, the OpenMx package in R is free and supports multilevel analyses, but requires a substantially different approach to syntax and specification.

This RMarkdown also leverages a package called MplusAutomation, which allows for effective cross-talk between Mplus and R. In particular, Mplus outputs all information into a text-based file, which makes it difficult to compare outputs across models or to extract specific subsections. If you’re interested in how this package works, please see Hallquist and Wiley (in press), Structural Equation Modeling.

1 What is multilevel analysis?

In the standard general linear model (GLM), we assume that observations are independent of each other. In particular, we assume that the residuals are distributed iid normal:

\[ \begin{align*} Y &= B_0 + B_1 X_1 + B_2 X_2 + ... B_p X_p + \varepsilon \\ \varepsilon &\sim \mathcal{N}(0, \sigma^2) \end{align*} \]

The assumption of independence is violated when:

  1. There are multiple observations per person (e.g., longitudinal data).
  2. There are known groups in the data that introduce similarity among group members.

Examples of multilevel/clustered data:

  1. Trials of an experimental task nested within individual (e.g., reaction times or accuracy data)
  2. A set of neurons (e.g., a population) measured in each animal, but with the goal of describing the pattern in a group of animals
  3. Longitudinal measurements of questionnaire or interview-based data
  4. Students within classrooms within schools (a so-called three-level model)
  5. Patients within therapy groups

In each case, the goal of multilevel analyses, including multilevel SEM, is to account for both within-cluster and between-cluster variability explicitly. If we were to analyze the data ignoring the group structure, we would overestimate the degrees of freedom (since our observations are non-independent). And we could not separate within- versus betewen-cluster effects; instead, our parameters would reflect an unknown mixture of the two.

1.1 L1 versus L2 effects

Multilevel analysis also allows the researcher to ask specific questions at each level of the model. Let’s consider the case of students (L1) nested within schools (L2). More specifically, let’s say we’re interested in what student- and school-level factors influence test performance. At the school level (aka L2), we hypothesize that the student:teacher ratio influences reading performance such that more students per teacher is associated with poorer performance. By definition, the student:teacher ratio is only meaningful as a school-level variable (setting aside the possibility of classroom variation) and has no variation at the person level. This might look like:

tt <- dplyr::tribble(~student, ~school, ~stratio, ~readtest, 1, "Radio Park", 20, 15, 2, "Radio Park", 20, 12, 3, "Radio Park", 20, 19, 4, "Corl Street", 16, 24, 5, "Corl Street", 16, 
    27, 6, "Corl Street", 16, 21)

kable(tt)
student school stratio readtest
1 Radio Park 20 15
2 Radio Park 20 12
3 Radio Park 20 19
4 Corl Street 16 24
5 Corl Street 16 27
6 Corl Street 16 21

Notice that there is no student-specific variation in the stratio. This is the simplest case for a predictor because there is no concern that it could represent a mixture of student and school effects. That is, stratio cannot predict any L1 variance.

We wish to represent the reading test performance for person i in school j as a function of reading speed. Some notation:

\[ \begin{align*} \textrm{readtest}_{ij} &= \beta_{0j} + r_{ij} \\ \beta_{0j} &= \gamma_{00} + \gamma_{01} \textrm{stratio}_j + \mu_{0j} \\ \end{align*} \] In this notation, \(\gamma_{01}\) represents the school-level effect of student:teacher ratio on test performance. We also allow for there to be normally distributed school-level variation in test performance, \(\mu_{0j}\) around the sample average \(\gamma_{00}\). To keep the \(\beta_{0j}\) equation easy to interpret, L2 predictors are typically centered around the grand mean. This ensures that \(\gamma_{00}\) retains the interpretation of the sample average test performance (fixed effect), and \(\mu_{0j}\) is school variation around this mean (random effect).

In addition, we might ask: what student characteristics enhance performance on standardized tests? These would be considered L1 predictors because they vary at the level of the student.

For example, students who read more words per minute may perform better on a standardized reading comprehension exam.

tt <- tt %>% add_column(readpermin = c(20, 14, 18, 25, 35, 30), .before = "school")
kable(tt)
student readpermin school stratio readtest
1 20 Radio Park 20 15
2 14 Radio Park 20 12
3 18 Radio Park 20 19
4 25 Corl Street 16 24
5 35 Corl Street 16 27
6 30 Corl Street 16 21

1.2 Disaggregating within- versus between-cluster effects

Importantly, the association between an L1 predictors and the outcome could reflect student and/or school variability. Perhaps students read more quickly because of an individual-level (aka L1) factor such as gray matter volume. Or might some schools have policies that promote better reading such that the average reading speed varies by school (i.e., at L2).

Because reading speed is an L1 variable (i.e., individual/student), it may have systematic variation at both the individual (L1) and school (L2) level. In conventional multilevel modeling (MLM), one can disentangle the variance at each level using appropriate centering (Curran & Bauer, 2011). In particular, the usual approach is to center the predictor (reading speed) around the school mean such that we now have a variable that represents deviation in reading speed relative to the school average. Then, we also compute the school means and include those as a predictor at L2. In this way, the effect of per-school reading speed on test performance can be examined specifically. Likewise, individual variation in reading speed not explainable by the school average can also be examined.

More notation from the MLM world.

\[ \begin{align*} \textrm{readtest}_{ij} &= \beta_{0j} + \beta_{1j} \textrm{readspeed}_{ij} + r_{ij} \\ \beta_{0j} &= \gamma_{00} + \gamma_{01} \textrm{stratio}_j + \mu_{0j} \\ \beta_{1j} &= \gamma_{10} + \mu_{1j} \end{align*} \]

The major difficulty of this representation is that \(\textrm{readspeed}_{ij}\) may have variation at the person and school levels. Thus, the main coefficient of interest for this predictor, \(\gamma_{10}\), (fixed effect) is an unknown combination of indiviual and school effects. Here, we also allow for the association of reading speed with reading test performance to vary randomly by school, \(\mu_{1j}\).

As mentioned above, to handle the problem of \(\textrm{readspeed}\) potentially having variation at L1 and L2, we apply the L1 centering + L2 means approach.

The notation for the school-centered reading speed is:

\[ \dot{\textrm{readspeed}_{ij}} = \textrm{readspeed}_{ij} - \bar{\textrm{readspeed}_{j}} \]

That is, we subtract the school mean reading speed, \(\bar{\textrm{readspeed}_{j}}\), from each student in that school.

We then introduce both of these new variables into the model:

\[ \begin{align*} \textrm{readtest}_{ij} &= \beta_{0j} + \beta_{1j} \dot{\textrm{readspeed}_{ij}} + r_{ij} \\ \beta_{0j} &= \gamma_{00} + \gamma_{01} \textrm{stratio}_j + \gamma_{02} \bar{\textrm{readspeed}_{j}} + \mu_{0j} \\ \beta_{1j} &= \gamma_{10} + \mu_{1j} \end{align*} \]

Thus, the school-average reading speed predicts school-level variation in test performance, whereas school-deviated individual variation predicts the association of individual reading speed on performance.

This is perhaps more evident in the combined equation:

\[ \textrm{readtest}_{ij} = ( \gamma_{00} + \gamma_{01} \textrm{stratio}_j + \gamma_{02} \bar{\textrm{readspeed}_{j}} + \gamma_{10} \dot{\textrm{readspeed}_{ij}} ) + \\ ( r_{ij} + \mu_{0j} + \mu_{1j} \dot{\textrm{readspeed}_{ij}} ) \] Note how \(\textrm{readspeed}\) enters the model in two ways (individual and school level). As in previous examples, the fixed effects are denoted in the first set of parentheses, and the random effects in the second.

2 Multilevel SEM (MSEM) overview

Okay, so we’ve got our head around standard multilevel models. Where does SEM enter the picture? In multilevel SEM, we use a latent variable approach to parcellate variation between and within clusters, rather than applying a cluster-based centering approach. This is the intuition:

In this representation, we view the observed variables as jointly caused by within- and between-cluster variation as latent variables.

Importantly, this framework also allows for the specification of random slopes such that within-person/cluster variation in the strength of an association between two variables can itself be a latent variable at the between-person level. For example, does the within-person association of reading speed with performance depend on student:teacher ratio? Said differently, does student:teacher ratio moderate the relationship between reading speed and performance? This is the multilevel SEM equivalent of cross-level interaction in MLM. Our graphical notation is as follows:

So far, however, the expression of multilevel data in a SEM framework has not revealed any new capabilities beyond standard MLM. Why would one consider multilevel SEM? First, some people find it easier to conceptualize of multilevel data structures and models within the SEM framework where we are used to path diagrams. Second, as in other non-multilevel analyses such as mediation, multilevel SEM offers the ability to test hypotheses involving several variables and paths. For an excellent introduction to the value of SEM for testing multilevel mediation, see Preacher, Zyphur, and Zhang (2010), Psychological Methods. Third, multilevel SEM allows us to test hypotheses involving latent variables, which are hard to implement in standard MLM. In particular, we can consider measurement models in which the covariation among a set of indicators is thought to reflect a latent variable. Fourth, the approach of cluster-mean centering introduces a potential bias in the between-group variance. In particular, having few level 1 observations and a low intraclass correlation will bias estimates of between-cluster effects (Lüdtke et al., 2008). This is especially problematic in cross-level mediation, where it is essential to treat the L1 portion of variance as a latent variable to recover unbiased parameters (Asparouhov & Muthén, 2006). In short, the latent decomposition of within-between in MSEM tends to provide better estimates than the centering method in MLM.

2.1 Sampling perspective

In MSEM, we have the option of whether to approach clustering from a (population) sampling perspective or a multilevel perspective. The sampling perspective is more typical in epidemiological or demographically representative studies in which we use stratified random sampling to cover the population of interest. This often involves computing sampling weights in which we might try to weight each set of observations so that the overall sample matches the composition of the population. For example, if we undersample African-American males, but wish to include them in the population to which we generalize, we might overweight those observations accordingly so that the parameters are more representative of the general population.

This problem also occurs in situations where individuals are nested in neighborhoods, cities, or states, but we wish to generalize to the overall population (e.g., adults in the USA). Again, we might reweight observations in the estimation in order to estimate a model that is representative of the population of interest. Thus, if we oversample Colorado and undersample Idaho, reweighting could help us recover unbiased parameters that correct for this problem.

According to the sampling perspective, clustering is an important feature of the data, but primarily because we wish to ‘uncluster’ the data in a representative fashion. That is, we are not interested in whether Colorado or Idaho has higher prevalence rates of depression, but instead we wish to estimate the prevalence of depression in the USA, accounting for the design of the epidemiological survey. If any/all of these questions pertain to your dataset or interest, see the lavaan.survey package, which is an add-on to lavaan that provides the functionality to handling design weights, clustering, and stratification. For example, we might wish to remove the effect of interviewer on prevalence estimates, rather than examine within-interviewer versus between-interviewer predictors. This could be done by setting up a svydesign object in which we account for clustering by interviewer.

FYI, the sampling perspective on clustered data is implements in Mplus using the TYPE=COMPLEX specification in the ANALYSIS section.

2.2 Multilevel perspective

By contrast, the multilevel perspective on clustered data is that we are specifically interested in parcelling variance at each level of the data structure so that we could examine predictors that explain variance at each level. For example, some researchers in our department (e.g., Dever Carney and Louis Castonguay) are interested in what are the client (L1), therapist (L2), and center/clinic (L3) factors that predict treatment outcomes in psychotherapy. As we talked about above, if a predictor exists at a lower level of the hierarchy (e.g., L1), it may nevertheless have systematic variance at higher levels. Thus, a client factor such as perceived therapeutic alliance may explain client-specific outcomes (L1), but alliance may also vary at the therapist level such that the clients of some therapists (L2) rate alliance higher, on average, relative to other therapists. Thus, even though we view alliance as the experience of a client, therapist-average alliance may explain why clients of some therapists do better than others. This is why the disaggregation approach mentioned above is crucial to resolving such questions: within-cluster centering, inclusion of cluster means.

In MSEM, we have the option of using the conventional centering strategy or use latent decomposition (as mentioned above) to disaggregate variance for a lower-level variable at each level.

In Mplus, we use TYPE=TWOLEVEL or TYPE=THREELEVEL in the ANALYSIS section to specify a multilevel model.

3 Multilevel SEM specification

The formal specification of multilevel SEM (Muthén & Asparouhov, 2008) is:

\[ \boldsymbol{Y}_{ij} = \boldsymbol{\nu_j} + \boldsymbol{\Lambda}_j\boldsymbol{\eta}_{ij} + \boldsymbol{K_j}\boldsymbol{X}_{ij} + \boldsymbol{\varepsilon}_{ij} \]

As in non-multilevel SEM, \(\boldsymbol{\nu}_j\) captures the intercepts of the variables, but allowing these to vary by cluster. And as usual, we have factor loadings (\(\boldsymbol{\Lambda}_j\)), factor scores (_{ij}), and the influence of exogenous covariates \(\boldsymbol{X}_{ij}\). The important point is that all of these can potentially vary by cluster. I’ll spare you the details on how MSEM parcels the structural and measurement models at the within and between levels… But the MSEM can be estimated using maximum likelihood in the typical fashion, with variance being assigned to the appropriate level by latent variable estimation.

4 A bit about Mplus MSEM syntax

Unlike lavaan, which is a package for R that follows most R conventions, Mplus is a standalone program that has its own syntax. Although Mplus cannot compete with the all-purpose functionality of R for data management and visualization, Mplus has developed a host of useful features for manipulating data, such as within-cluster centering or computing interaction terms.

The basic paradigm of Mplus is that a model consists of three parts: data, model syntax, and model output. By convention, the raw data are provided as a tab-separated file with no column headers, typically using the extension .dat. The Mplus model syntax has its own conventions, such as the use of ON to denote the regression of Y ON X. Likewise, WITH denotes the undirected (residual) association of X WITH Y. And factors are specified using BY, as in ‘measured by’: F1 BY X1 X2 X3. Mplus syntax also has major sections that define important components of the model, such as the estimator, what data are categorical, what additional outputs to provide, and so on.

By convention, Mplus syntax files are stored as text with the extension .inp. And the output provided by Mplus is also text, stored as .out. There is an impressive array of information in the output file, but its storage as text makes it time-intensive to sift through each section. This is where the readModels() command from MplusAutomation can help index a number of output files by storing the information in R-friendly data structures such as lists and data.frames.

4.1 WITHIN

If a variable is specified as WITHIN in the VARIABLE section, then it is assumed to vary only within clusters and to have no between component. In the latent decomposition, the between component of the variable is zero. Thus, this statement is used when 1) the variable varies at the within-cluster level, and 2) we wish to assign the totality of its variance at the within-cluster level, assuming that there is no between component.

4.2 BETWEEN

If a variable is specified as BETWEEN in the VARIABLE section, it is a cluster-level variable that has no variation within-cluster. These are typically called L2 variables. For example, student:teacher ratio would be a BETWEEN variable if we are examining students nested in schools.

4.3 Neither WITHIN nor BETWEEN

If we do not specify that a variable is WITHIN or BETWEEN then Mplus will model its variance at both levels. That is, it will use latent decomposition to parcellate the variance into within- and between-cluster variance. At the between-cluster level, this means we are modeling the random intercept of that variable.

5 Getting started with MSEM in Mplus

Let’s start with a simple multilevel model that we could fit using conventional MLM software. This is Example 9.1 from the Mplus user’s guide.

5.1 Centering approach to disaggregation

The dot on the \(y\) box denotes a random intercept of y. That is, there is between-cluster variation in the average level of \(y\) that can be modeled at the between level. \(x\) is the person-centered L1 predictor. \(xm\) is the person mean on \(x\). \(w\) is a cluster-level covariate (i.e., no within variance).

TITLE:  this is an example of a two-level 
    regression analysis for a continuous 
    dependent variable with a random intercept and an observed covariate
DATA:   FILE = ex9.1a.dat;
VARIABLE:   NAMES = y x w xm clus;
    WITHIN = x;
    BETWEEN = w xm;
    CLUSTER = clus;
DEFINE: CENTER x (GRANDMEAN);
ANALYSIS:   TYPE = TWOLEVEL;
MODEL:
    %WITHIN%    
    y ON x;
    %BETWEEN%
    y ON w xm;
ex9.1a <- read.table("ex9.1a.dat", col.names = c("y", "x", "w", "xm", "clus"))
kable(head(ex9.1a, n = 15))
y x w xm clus
1.221 0.340 -2.558 -1.368 1
-0.606 -0.797 -2.558 -1.368 1
-0.290 0.014 -2.558 -1.368 1
-1.605 0.429 -2.558 -1.368 1
0.787 0.813 -2.558 -1.368 1
2.391 -2.649 0.971 0.043 2
2.175 -0.647 0.971 0.043 2
3.916 0.303 0.971 0.043 2
3.836 0.471 0.971 0.043 2
2.284 0.382 0.971 0.043 2
0.864 0.695 -0.945 0.209 3
0.490 -0.552 -0.945 0.209 3
1.420 -0.618 -0.945 0.209 3
0.668 -1.122 -0.945 0.209 3
2.481 0.598 -0.945 0.209 3 
ex9.1a_mobj <- mplusObject(TITLE = "this is an example of a two-level 
\t  regression analysis for a continuous 
\t  dependent variable with a random intercept and an observed covariate", 
    VARIABLE = "WITHIN = x;
\t  BETWEEN = w xm;
\t  CLUSTER = clus;", DEFINE = "CENTER x (GRANDMEAN);", ANALYSIS = "TYPE = TWOLEVEL;
    ESTIMATOR=MLR;", MODEL = "%WITHIN%\t
\t  y ON x;
\t  %BETWEEN%
\t  y ON w xm;", 
    OUTPUT = "STDYX RESIDUAL;", rdata = ex9.1a)
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
## If any issues, suggest explicitly specifying USEVARIABLES.
## A starting point may be:
## USEVARIABLES = y x w xm clus;
ex9.1a_fit <- mplusModeler(ex9.1a_mobj, modelout = "ex9.1a.inp", run = TRUE, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: ex9.1a.inp
## Wrote data to: ex9.1a.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'ex9.1a.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: ex9.1a.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "ex9.1a.inp" 
## Reading model:  ex9.1a.out
screenreg(ex9.1a_fit$results, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"))
## 
## ===================================
##                  Model 1           
## -----------------------------------
## W Y<-X              0.72 (0.03) ***
## B Y<-W              0.57 (0.11) ***
## B Y<-XM             0.98 (0.16) ***
## B Y<-Intercepts     1.99 (0.08) ***
## W Y<->Y             1.02 (0.04) ***
## B Y<->Y             0.57 (0.09) ***
## -----------------------------------
## AICC             3063.96           
## CFI                 1.00           
## SRMR.Within         0.00           
## SRMR.Between        0.00           
## ===================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
# examine w/i and b/w data structure ex9.1a <- ex9.1a %>% group_by(clus) %>% mutate(xm_mine=mean(x), xwi=x - xm_mine) round(cor(ex9.1), 3)

Let’s compare that to the lmer function in R:

ex9.1a_lmer <- lmer(y ~ x + w + xm + (1 | clus), ex9.1a, REML = FALSE)
summary(ex9.1a_lmer)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: y ~ x + w + xm + (1 | clus)
##    Data: ex9.1a
## 
##      AIC      BIC   logLik deviance df.resid 
##     3064     3093    -1526     3052      994 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7739 -0.6851  0.0013  0.6841  2.6850 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  clus     (Intercept) 0.57     0.755   
##  Residual             1.02     1.011   
## Number of obs: 1000, groups:  clus, 110
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   1.9753     0.0799   24.72
## x             0.7241     0.0337   21.49
## w             0.5703     0.1066    5.35
## xm            0.9761     0.1594    6.12
## 
## Correlation of Fixed Effects:
##    (Intr) x      w     
## x  -0.008              
## w  -0.061 -0.003       
## xm  0.030  0.026 -0.709

You can see that the results are essentially identical.

5.2 Latent decomposition approach to disaggregation

Unlike MLM, MSEM also offers the potential to consider variation in an observed variable as joint function of within- and between-cluster latent components. Let’s re-run this analysis, using the decomposition approach.

TITLE:  this is an example of a two-level 
    regression analysis for a continuous 
    dependent variable with a random intercept and a latent covariate
DATA:   FILE = ex9.1b.dat;
VARIABLE:   NAMES = y x w clus;
    BETWEEN = w;
    CLUSTER = clus;
DEFINE: CENTER x (GRANDMEAN);
ANALYSIS:   TYPE = TWOLEVEL;
MODEL:
    %WITHIN%    
    y ON x (gamma10);
    %BETWEEN%
    y ON w 
    x (gamma01);
MODEL CONSTRAINT:
    NEW(betac);
    betac = gamma01 - gamma10;
ex9.1b <- read.table("ex9.1b.dat", col.names = c("y", "x", "w", "clus"))
kable(head(ex9.1b, n = 15))
y x w clus
0.995 0.109 -2.36 1
-0.341 -0.129 -2.36 1
-0.512 -1.646 -2.36 1
2.095 -0.542 -2.36 1
-0.167 -0.678 -2.36 1
1.807 -0.694 -1.22 2
1.675 0.784 -1.22 2
3.047 0.582 -1.22 2
0.964 -0.482 -1.22 2
1.014 -0.928 -1.22 2
1.136 0.980 1.29 3
2.402 0.063 1.29 3
0.519 -0.705 1.29 3
2.939 -0.529 1.29 3
0.913 0.353 1.29 3 
ex9.1b_mobj <- mplusObject(TITLE = "this is an example of a two-level 
\t  regression analysis for a continuous 
\t  dependent variable with a random intercept and a latent covariate", 
    VARIABLE = "BETWEEN = w;
\t  CLUSTER = clus;", DEFINE = "CENTER x (GRANDMEAN);", ANALYSIS = "TYPE = TWOLEVEL;
    ESTIMATOR=MLR;", MODEL = "
    %WITHIN%
\t    y ON x (gamma10);
\t  %BETWEEN%
\t    y ON w 
\t         x (gamma01);", 
    OUTPUT = "STDYX RESIDUAL;", rdata = ex9.1b)
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
## If any issues, suggest explicitly specifying USEVARIABLES.
## A starting point may be:
## USEVARIABLES = y x w clus;
ex9.1b_fit <- mplusModeler(ex9.1b_mobj, modelout = "ex9.1b.inp", run = TRUE, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: ex9.1b.inp
## Wrote data to: ex9.1b.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'ex9.1b.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: ex9.1b.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "ex9.1b.inp" 
## Reading model:  ex9.1b.out
screenreg(ex9.1b_fit$results, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"))
## 
## ===================================
##                  Model 1           
## -----------------------------------
## W Y<-X              0.78 (0.03) ***
## B Y<-W              0.38 (0.12) ** 
## B Y<-X              1.19 (0.18) ***
## B Y<-Intercepts     2.00 (0.08) ***
## W Y<->Y             0.98 (0.04) ***
## B Y<->Y             0.52 (0.08) ***
## -----------------------------------
## AICC             5955.25           
## CFI                 1.00           
## SRMR.Within         0.00           
## SRMR.Between        0.00           
## ===================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
# ex9.1 <- ex9.1 %>% group_by(clus) %>% mutate(xm_mine=mean(x), xwi=x - xm_mine) ex9.1 %>% group_by(clus) %>%

5.3 Random slopes in MSEM

As in standard MLM, we might expect there to be between-person variation in the strength of an association. For example, are there meaningful interindividual differences in the rate of change in age-related declines in memory? In MSEM, we need to specify these random slopes as latent variables that vary at L2.

Here is example 9.2 from the users guide using the within-between latent disaggregation for \(x\), an L1 predictor.

TITLE:  this is an example of a two-level 
    regression analysis for a continuous 
    dependent variable with a random slope and a latent covariate
DATA:   FILE = ex9.2c.dat;
VARIABLE:   NAMES = y x w clus;
    BETWEEN = w;
    CLUSTER = clus;
ANALYSIS:   TYPE = TWOLEVEL RANDOM;
MODEL:
    %WITHIN%    
    s | y ON x;     
    %BETWEEN%   
    y s ON w x;
    y WITH s;
ex9.2c <- read.table("ex9.2c.dat", col.names = c("y", "x", "w", "clus"))
kable(head(ex9.2c, n = 15))
y x w clus
-0.601 -0.215 -2.364 1
-2.188 0.587 -2.364 1
-1.053 0.232 -2.364 1
-1.161 -0.094 -2.364 1
-1.003 -0.126 -2.364 1
0.798 -0.883 -0.663 2
-2.493 -2.123 -0.663 2
1.207 -0.600 -0.663 2
0.067 -1.884 -0.663 2
1.510 -1.500 -0.663 2
1.731 -0.875 -0.183 3
0.357 -1.643 -0.183 3
2.661 -1.467 -0.183 3
0.050 -0.585 -0.183 3
0.909 0.409 -0.183 3 
ex9.2c_mobj <- mplusObject(TITLE = "this is an example of a two-level 
\tregression analysis for a continuous 
\tdependent variable with a random slope and a latent covariate", VARIABLE = "BETWEEN = w;
\t  CLUSTER = clus;", 
    DEFINE = "CENTER x (GRANDMEAN);", ANALYSIS = "TYPE = TWOLEVEL RANDOM;
    ESTIMATOR=MLR;", MODEL = "
    %WITHIN%\t
 \t    s | y ON x;\t\t
\t  %BETWEEN%\t
\t    y s ON w x;
\t    y WITH s;", 
    OUTPUT = "STDYX RESIDUAL;", rdata = ex9.2c)
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
## If any issues, suggest explicitly specifying USEVARIABLES.
## A starting point may be:
## USEVARIABLES = y x w clus;
ex9.2c_fit <- mplusModeler(ex9.2c_mobj, modelout = "ex9.2c.inp", run = TRUE, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: ex9.2c.inp
## Wrote data to: ex9.2c.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'ex9.2c.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: ex9.2c.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "ex9.2c.inp" 
## Reading model:  ex9.2c.out
screenreg(ex9.2c_fit$results, single.row = TRUE, summaries = c("AICC"))
## 
## ===================================
##                  Model 1           
## -----------------------------------
## B S<-W              0.57 (0.09) ***
## B S<-X              0.32 (0.18)    
## B Y<-W              1.13 (0.11) ***
## B Y<-X              0.99 (0.21) ***
## B Y<->S             0.23 (0.06) ***
## B Y<-Intercepts     1.90 (0.08) ***
## B S<-Intercepts     0.99 (0.07) ***
## W Y<->Y             1.03 (0.05) ***
## B Y<->Y             0.44 (0.09) ***
## B S<->S             0.37 (0.06) ***
## -----------------------------------
## AICC             6197.21           
## ===================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6 Hox Family IQ MSEM example

These data are from Joop Hox’s multilevel book. The basic structure is that we have six intelligence measures for 400 children from 60 families. Thus, we have children nested within families. If intelligence has a genetic or heritable component, the between-family differences may be of particular interest. The measures are word list, cards, matrices, figures, animals, and occupations. One possibility is that the intelligence tests all form a single g factor at the between-family level such that some families have higher average scores on intelligence measures than others. Another possibility is that the subscales form two factor: numeric (word list, cards, matrices) versus perceptual (figures, animals, occupations) intelligence.

As in the case of standard MLM, in MSEM, we often wish to examine the amount of variation at each level in an unconditional model. For SEM purposes, examining the independence and saturated models is informative. As we’ve talked about in standard SEM, the independence model gives us a worst case fit where we do not permit item covariation. And the saturated model gives us perfect fit without any parsimony or prediction. This is true of MSEM as well, but we can examine the spectrum between independence and saturated at both the within and between levels.

6.1 Independence model

To get a sense of whether there is structured variance in IQ due to family, we can estimate an independence model for the data in which the IQ variables are assumed to be unrelated at both the within- and between-family levels. This gives us a benchmark for variation at these levels.

Here is the raw Mplus syntax

TITLE:  Within independence; Between independence;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
 %BETWEEN%
OUTPUT: STDYX SAMPSTAT;

Note that I’ve printed out the unstandardized estimates here to give a sense of variation in each variable in the original unit. Otherwise, we’d see variances of 1.0 at both levels for the standardized estimates. For other models, however, I primarily report completely standardized effects (std.all in lavaan terms) for interpretability.

# import data from spss
hoxdata <- haven::read_spss("FamIQData.sav")

mobj1 <- mplusObject(TITLE = "Within independence; Between independence;", VARIABLE = "CLUSTER = family;
    USEVARIABLES ARE wordlist cards matrices figures animals occupats;", ANALYSIS = "TYPE IS TWOLEVEL;
    ESTIMATOR IS MLR;", 
    MODEL = "%WITHIN%
    %BETWEEN%", OUTPUT = "STDYX SAMPSTAT;", rdata = hoxdata)
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit1 <- mplusModeler(mobj1, modelout = "hox_m1.inp", run = 1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m1.inp
## Wrote data to: hox_m1.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m1.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m1.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m1.inp" 
## Reading model:  hox_m1.out
## Warning in detectColumnNames(filename, modelSection, sectionType): Unable to determine column names for section model_results.
##   hox_m1.out
screenreg(mfit1, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"))
## 
## ==========================================
##                        Model 1            
## ------------------------------------------
## B WORDLIST<-Means         29.86 (0.47) ***
## B CARDS<-Means            29.89 (0.46) ***
## B MATRICES<-Means         29.75 (0.44) ***
## B FIGURES<-Means          30.06 (0.45) ***
## B ANIMALS<-Means          30.18 (0.47) ***
## B OCCUPATS<-Means         29.98 (0.51) ***
## W WORDLIST<->WORDLIST     16.16 (1.98) ***
## W CARDS<->CARDS           15.26 (1.45) ***
## W MATRICES<->MATRICES     15.61 (1.40) ***
## W FIGURES<->FIGURES       16.55 (1.25) ***
## W ANIMALS<->ANIMALS       15.08 (1.18) ***
## W OCCUPATS<->OCCUPATS     13.24 (1.04) ***
## B WORDLIST<->WORDLIST     10.59 (2.11) ***
## B CARDS<->CARDS           10.47 (2.47) ***
## B MATRICES<->MATRICES      8.94 (1.76) ***
## B FIGURES<->FIGURES        9.78 (2.04) ***
## B ANIMALS<->ANIMALS       10.82 (2.40) ***
## B OCCUPATS<->OCCUPATS     13.44 (2.76) ***
## ------------------------------------------
## AICC                   14007.16           
## CFI                        0.00           
## SRMR.Within                0.36           
## SRMR.Between               0.64           
## ==========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.2 Saturated model

At the other end of the spectrum, we could estimate a saturated model for both within- and between-cluster variation.

Here is the raw Mplus syntax:

TITLE:  Within saturated, Between saturated;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  wordlist cards matrices figures animals occupats WITH
  wordlist cards matrices figures animals occupats;
 %BETWEEN%
  wordlist cards matrices figures animals occupats WITH
  wordlist cards matrices figures animals occupats; 
OUTPUT: STDYX SAMPSTAT;

These are the standardized covariances, aka correlations.

mobj2 <- update(mobj1, TITLE = ~"W/i saturated, B/w saturated;", MODEL = ~"
    %WITHIN%
      wordlist cards matrices figures animals occupats WITH
        wordlist cards matrices figures animals occupats;
    %BETWEEN%
      wordlist cards matrices figures animals occupats WITH
        wordlist cards matrices figures animals occupats;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit2 <- mplusModeler(mobj2, modelout = "hox_m2.inp", run = 1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m2.inp
## Wrote data to: hox_m2.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m2.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m2.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m2.inp" 
## Reading model:  hox_m2.out
## Warning in detectColumnNames(filename, modelSection, sectionType): Unable to determine column names for section model_results.
##   hox_m2.out
screenreg(mfit2, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type = "stdyx")
## 
## ==========================================
##                        Model 1            
## ------------------------------------------
## W WORDLIST<->CARDS         0.63 (0.04) ***
## W WORDLIST<->MATRICES      0.60 (0.04) ***
## W WORDLIST<->FIGURES       0.23 (0.05) ***
## W WORDLIST<->ANIMALS       0.30 (0.06) ***
## W WORDLIST<->OCCUPATS      0.24 (0.05) ***
## W CARDS<->MATRICES         0.62 (0.04) ***
## W CARDS<->FIGURES          0.18 (0.05) ** 
## W CARDS<->ANIMALS          0.29 (0.04) ***
## W CARDS<->OCCUPATS         0.24 (0.05) ***
## W MATRICES<->FIGURES       0.18 (0.06) ** 
## W MATRICES<->ANIMALS       0.24 (0.05) ***
## W MATRICES<->OCCUPATS      0.19 (0.06) ** 
## W FIGURES<->ANIMALS        0.63 (0.04) ***
## W FIGURES<->OCCUPATS       0.60 (0.04) ***
## W ANIMALS<->OCCUPATS       0.63 (0.03) ***
## B WORDLIST<->CARDS         0.87 (0.05) ***
## B WORDLIST<->MATRICES      0.86 (0.04) ***
## B WORDLIST<->FIGURES       0.82 (0.08) ***
## B WORDLIST<->ANIMALS       0.89 (0.06) ***
## B WORDLIST<->OCCUPATS      0.87 (0.06) ***
## B CARDS<->MATRICES         0.82 (0.06) ***
## B CARDS<->FIGURES          0.84 (0.07) ***
## B CARDS<->ANIMALS          0.92 (0.04) ***
## B CARDS<->OCCUPATS         0.88 (0.06) ***
## B MATRICES<->FIGURES       0.83 (0.07) ***
## B MATRICES<->ANIMALS       0.81 (0.08) ***
## B MATRICES<->OCCUPATS      0.82 (0.06) ***
## B FIGURES<->ANIMALS        0.86 (0.05) ***
## B FIGURES<->OCCUPATS       0.83 (0.05) ***
## B ANIMALS<->OCCUPATS       0.91 (0.04) ***
## B WORDLIST<-Means          9.12 (0.96) ***
## B CARDS<-Means             9.21 (1.11) ***
## B MATRICES<-Means          9.85 (0.96) ***
## B FIGURES<-Means           9.57 (1.02) ***
## B ANIMALS<-Means           9.14 (1.04) ***
## B OCCUPATS<-Means          8.19 (0.88) ***
## W WORDLIST<->WORDLIST      1.00 (0.00)    
## W CARDS<->CARDS            1.00 (0.00)    
## W MATRICES<->MATRICES      1.00 (0.00)    
## W FIGURES<->FIGURES        1.00 (0.00)    
## W ANIMALS<->ANIMALS        1.00 (0.00)    
## W OCCUPATS<->OCCUPATS      1.00 (0.00)    
## B WORDLIST<->WORDLIST      1.00 (0.00)    
## B CARDS<->CARDS            1.00 (0.00)    
## B MATRICES<->MATRICES      1.00 (0.00)    
## B FIGURES<->FIGURES        1.00 (0.00)    
## B ANIMALS<->ANIMALS        1.00 (0.00)    
## B OCCUPATS<->OCCUPATS      1.00 (0.00)    
## ------------------------------------------
## AICC                   12898.10           
## CFI                        1.00           
## SRMR.Within                0.00           
## SRMR.Between               0.00           
## ==========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.3 One-factor between-family, one-factor within-family

What about a simple/parsimonious model in which we assume a one-factor structure at the within and between levels?

TITLE:  Within one factor; Between one factor;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  gen_wi BY wordlist* cards matrices figures animals occupats;
  gen_wi@1;
 %BETWEEN%
  gen_bw BY wordlist* cards matrices figures animals occupats;
  gen_bw@1;
OUTPUT: STDYX SAMPSTAT;

Note the Heywood case here: there’s a negative residual on the animals variable at the between level. One strategy would be to add animals@0; to the between syntax to fix the residual variance to 0, but the case is pretty bad in the sense that the negative variance is significant, not trivially below zero. Given how badly the model fits at the within-person level (see SRMR.Within), I would probably give up on this and move on to other models.

mobj3 <- update(mobj1, TITLE = ~"W/i one factor; B/w one factor;", MODEL = ~"
    %WITHIN%
      gen_wi BY wordlist* cards matrices figures animals occupats;
      gen_wi@1;
    %BETWEEN%
      gen_bw BY wordlist* cards matrices figures animals occupats;
      gen_bw@1;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit3 <- mplusModeler(mobj3, modelout = "hox_m3.inp", run = 1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m3.inp
## Wrote data to: hox_m3.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m3.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m3.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m3.inp" 
## Reading model:  hox_m3.out
## Warning in rbind(character(0), c("WORDLIST", "0.836", "0.054", "15.582", : number of columns of result is not a multiple of vector length (arg 6)
screenreg(mfit3, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type = "stdyx")
## 
## =============================================
##                         Model 1              
## ---------------------------------------------
## W WORDLIST<-GEN_WI          0.79   (0.03) ***
## W CARDS<-GEN_WI             0.80   (0.03) ***
## W MATRICES<-GEN_WI          0.76   (0.04) ***
## W FIGURES<-GEN_WI           0.35   (0.07) ***
## W ANIMALS<-GEN_WI           0.42   (0.06) ***
## W OCCUPATS<-GEN_WI          0.38   (0.07) ***
## B WORDLIST<-GEN_BW          0.91   (0.03) ***
## B CARDS<-GEN_BW             0.92   (0.03) ***
## B MATRICES<-GEN_BW          0.86   (0.04) ***
## B FIGURES<-GEN_BW           0.94   (0.01) ***
## B ANIMALS<-GEN_BW           1.01   (0.00) ***
## B OCCUPATS<-GEN_BW          0.98   (0.00) ***
## B WORDLIST<-Intercepts      9.52   (0.63) ***
## B CARDS<-Intercepts         9.52   (0.61) ***
## B MATRICES<-Intercepts     10.37   (0.72) ***
## B FIGURES<-Intercepts       9.40   (0.59) ***
## B ANIMALS<-Intercepts       9.08   (0.17) ***
## B OCCUPATS<-Intercepts      8.19   (0.32) ***
## W GEN_WI<->GEN_WI           1.00   (0.00)    
## W WORDLIST<->WORDLIST       0.38   (0.05) ***
## W CARDS<->CARDS             0.37   (0.05) ***
## W MATRICES<->MATRICES       0.43   (0.05) ***
## W FIGURES<->FIGURES         0.88   (0.05) ***
## W ANIMALS<->ANIMALS         0.82   (0.05) ***
## W OCCUPATS<->OCCUPATS       0.86   (0.05) ***
## B GEN_BW<->GEN_BW           1.00   (0.00)    
## B WORDLIST<->WORDLIST       0.16   (0.05) ** 
## B CARDS<->CARDS             0.16   (0.05) ** 
## B MATRICES<->MATRICES       0.26   (0.07) ***
## B FIGURES<->FIGURES         0.11   (0.02) ***
## B ANIMALS<->ANIMALS        -0.02 (999.00)    
## B OCCUPATS<->OCCUPATS       0.05   (0.00) ***
## ---------------------------------------------
## AICC                    13188.09             
## CFI                         0.78             
## SRMR.Within                 0.19             
## SRMR.Between                0.05             
## =============================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.4 One-factor between-family, saturated within-family model

TITLE:  Within saturated--Between one factor; Table 14.7, p. 308;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  wordlist cards matrices figures animals occupats WITH
  wordlist cards matrices figures animals occupats;
 %BETWEEN%
  general BY wordlist* cards matrices figures animals occupats;
  general@1;
OUTPUT: STDYX SAMPSTAT;
mobj4 <- update(
  mobj1, 
  TITLE = ~ "W/i saturated; B/w one factor;", # Table 14.7, p. 308
  MODEL = ~ "%WITHIN%
  wordlist cards matrices figures animals occupats WITH
  wordlist cards matrices figures animals occupats;
 %BETWEEN%
  general BY wordlist* cards matrices figures animals occupats;
  general@1;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit4 <- mplusModeler(mobj4, modelout = "hox_m4.inp", run=1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m4.inp
## Wrote data to: hox_m4.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m4.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m4.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m4.inp" 
## Reading model:  hox_m4.out
screenreg(mfit4, single.row=TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type="stdyx")
## 
## ===========================================
##                         Model 1            
## -------------------------------------------
## B WORDLIST<-GENERAL         0.94 (0.03) ***
## B CARDS<-GENERAL            0.94 (0.04) ***
## B MATRICES<-GENERAL         0.89 (0.04) ***
## B FIGURES<-GENERAL          0.89 (0.04) ***
## B ANIMALS<-GENERAL          0.97 (0.03) ***
## B OCCUPATS<-GENERAL         0.94 (0.02) ***
## W WORDLIST<->CARDS          0.63 (0.04) ***
## W WORDLIST<->MATRICES       0.61 (0.04) ***
## W WORDLIST<->FIGURES        0.22 (0.05) ***
## W WORDLIST<->ANIMALS        0.29 (0.06) ***
## W WORDLIST<->OCCUPATS       0.23 (0.05) ***
## W CARDS<->MATRICES          0.62 (0.04) ***
## W CARDS<->FIGURES           0.18 (0.05) ** 
## W CARDS<->ANIMALS           0.29 (0.04) ***
## W CARDS<->OCCUPATS          0.24 (0.05) ***
## W MATRICES<->FIGURES        0.18 (0.06) ** 
## W MATRICES<->ANIMALS        0.23 (0.05) ***
## W MATRICES<->OCCUPATS       0.19 (0.06) ** 
## W FIGURES<->ANIMALS         0.63 (0.04) ***
## W FIGURES<->OCCUPATS        0.60 (0.04) ***
## W ANIMALS<->OCCUPATS        0.63 (0.03) ***
## B WORDLIST<-Intercepts      9.17 (0.98) ***
## B CARDS<-Intercepts         9.18 (1.11) ***
## B MATRICES<-Intercepts      9.94 (0.98) ***
## B FIGURES<-Intercepts       9.53 (0.98) ***
## B ANIMALS<-Intercepts       9.17 (1.06) ***
## B OCCUPATS<-Intercepts      8.20 (0.88) ***
## W WORDLIST<->WORDLIST       1.00 (0.00)    
## W CARDS<->CARDS             1.00 (0.00)    
## W MATRICES<->MATRICES       1.00 (0.00)    
## W FIGURES<->FIGURES         1.00 (0.00)    
## W ANIMALS<->ANIMALS         1.00 (0.00)    
## W OCCUPATS<->OCCUPATS       1.00 (0.00)    
## B GENERAL<->GENERAL         1.00 (0.00)    
## B WORDLIST<->WORDLIST       0.12 (0.06) *  
## B CARDS<->CARDS             0.12 (0.07)    
## B MATRICES<->MATRICES       0.21 (0.07) ** 
## B FIGURES<->FIGURES         0.21 (0.07) ** 
## B ANIMALS<->ANIMALS         0.07 (0.05)    
## B OCCUPATS<->OCCUPATS       0.12 (0.05) *  
## -------------------------------------------
## AICC                    12880.04           
## CFI                         1.00           
## SRMR.Within                 0.00           
## SRMR.Between                0.02           
## ===========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.5 Within two-factor, Between saturated

TITLE: Within two factor--Between saturated;
 Table 14.1, p. 301 & Table 14.7, p. 308;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  numeric BY wordlist* cards matrices;
  percept BY figures* animals occupats;
  numeric@1 percept@1;
 %BETWEEN%
  wordlist cards matrices figures animals occupats WITH
  wordlist cards matrices figures animals occupats; 
OUTPUT: STDYX SAMPSTAT;
mobj5 <- update(
  mobj1, 
  TITLE = ~ "Within two factor--Between saturated;", # Table 14.1, p. 301 & Table 14.7, p. 308
  MODEL = ~ "%WITHIN%
    numeric BY wordlist* cards matrices;
    percept BY figures* animals occupats;
    numeric@1 percept@1;
   %BETWEEN%
    wordlist cards matrices figures animals occupats WITH
    wordlist cards matrices figures animals occupats;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit5 <- mplusModeler(mobj5, modelout = "hox_m5.inp", run=1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m5.inp
## Wrote data to: hox_m5.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m5.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m5.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m5.inp" 
## Reading model:  hox_m5.out
screenreg(mfit5, single.row=TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type="stdyx")
## 
## ==========================================
##                        Model 1            
## ------------------------------------------
## W WORDLIST<-NUMERIC        0.78 (0.04) ***
## W CARDS<-NUMERIC           0.81 (0.03) ***
## W MATRICES<-NUMERIC        0.77 (0.03) ***
## W FIGURES<-PERCEPT         0.77 (0.04) ***
## W ANIMALS<-PERCEPT         0.82 (0.02) ***
## W OCCUPATS<-PERCEPT        0.77 (0.03) ***
## W PERCEPT<->NUMERIC        0.38 (0.05) ***
## B WORDLIST<->CARDS         0.87 (0.05) ***
## B WORDLIST<->MATRICES      0.86 (0.04) ***
## B WORDLIST<->FIGURES       0.82 (0.08) ***
## B WORDLIST<->ANIMALS       0.90 (0.06) ***
## B WORDLIST<->OCCUPATS      0.87 (0.06) ***
## B CARDS<->MATRICES         0.82 (0.06) ***
## B CARDS<->FIGURES          0.82 (0.07) ***
## B CARDS<->ANIMALS          0.93 (0.04) ***
## B CARDS<->OCCUPATS         0.88 (0.06) ***
## B MATRICES<->FIGURES       0.82 (0.07) ***
## B MATRICES<->ANIMALS       0.81 (0.08) ***
## B MATRICES<->OCCUPATS      0.82 (0.06) ***
## B FIGURES<->ANIMALS        0.86 (0.05) ***
## B FIGURES<->OCCUPATS       0.84 (0.05) ***
## B ANIMALS<->OCCUPATS       0.91 (0.04) ***
## B WORDLIST<-Means          9.12 (0.96) ***
## B CARDS<-Means             9.20 (1.11) ***
## B MATRICES<-Means          9.86 (0.96) ***
## B FIGURES<-Means           9.55 (1.01) ***
## B ANIMALS<-Means           9.13 (1.03) ***
## B OCCUPATS<-Means          8.19 (0.88) ***
## W NUMERIC<->NUMERIC        1.00 (0.00)    
## W PERCEPT<->PERCEPT        1.00 (0.00)    
## W WORDLIST<->WORDLIST      0.38 (0.06) ***
## W CARDS<->CARDS            0.35 (0.05) ***
## W MATRICES<->MATRICES      0.41 (0.05) ***
## W FIGURES<->FIGURES        0.41 (0.05) ***
## W ANIMALS<->ANIMALS        0.32 (0.04) ***
## W OCCUPATS<->OCCUPATS      0.40 (0.05) ***
## B WORDLIST<->WORDLIST      1.00 (0.00)    
## B CARDS<->CARDS            1.00 (0.00)    
## B MATRICES<->MATRICES      1.00 (0.00)    
## B FIGURES<->FIGURES        1.00 (0.00)    
## B ANIMALS<->ANIMALS        1.00 (0.00)    
## B OCCUPATS<->OCCUPATS      1.00 (0.00)    
## ------------------------------------------
## AICC                   12884.68           
## CFI                        1.00           
## SRMR.Within                0.02           
## SRMR.Between               0.00           
## ==========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

Detailed comparison of within two-factor model against saturated model, including nested chi-square test.

compareModels(mfit2$results, mfit5$results, diffTest = TRUE)
## 
## ==============
## 
## Mplus model comparison
## ----------------------
## 
## ------
## Model 1:  hox_m2.out 
## Model 2:  hox_m5.out 
## ------
## 
## Model Summary Comparison
## ------------------------
## 
##                m1                            m2                                   
## Title          W/i saturated, B/w saturated; Within two factor--Between saturated;
## Observations   400                           400                                  
## Estimator      MLR                           MLR                                  
## Parameters     48                            40                                   
## LL             -6394.348                     -6397.77                             
## AIC            12884.697                     12875.54                             
## BIC            13076.287                     13035.198                            
## ChiSqM_Value   0                             6.832                                
## ChiSqM_DF      0                             8                                    
## CFI            1                             1                                    
## TLI            1                             1.004                                
## RMSEA_Estimate 0                             0                                    
## 
##   MLR Chi-Square Difference Test for Nested Models Based on Loglikelihood
##   -----------------------------------------------------------------------
## 
##   Difference Test Scaling Correction:  1.01 
##   Chi-square difference:  6.8 
##   Diff degrees of freedom:  8 
##   P-value:  0.559 
## 
##   Note: The chi-square difference test assumes that these models are nested.
##   It is up to you to verify this assumption.
## 
##   MLR Chi-Square Difference test for nested models
##   --------------------------------------------
## 
##   Difference Test Scaling Correction: 1 
##   Chi-square difference: 6.83 
##   Diff degrees of freedom: 8 
##   P-value: 0.555 
## 
## Note: The chi-square difference test assumes that these models are nested.
##   It is up to you to verify this assumption.
## 
## =========
## 
## Model parameter comparison
## --------------------------
##   Parameters present in both models
## =========
## 
##   Approximately equal in both models (param. est. diff <= 1e-04)
##   ----------------------------------------------
## None
## 
## 
##   Parameter estimates that differ between models (param. est. diff > 1e-04)
##   ----------------------------------------------
##    paramHeader    param BetweenWithin m1_est m2_est . m1_se m2_se . m1_est_se m2_est_se . m1_pval m2_pval
##   ANIMALS.WITH OCCUPATS       Between  11.01  11.01 | 2.373 2.373 |      4.64      4.64 |       0       0
##     CARDS.WITH  ANIMALS       Between   9.89  10.02 | 2.322 2.319 |      4.26      4.32 |       0       0
##     CARDS.WITH  FIGURES       Between   8.57   8.44 | 2.095 2.115 |      4.09      3.99 |       0       0
##     CARDS.WITH MATRICES       Between   7.98   7.99 | 1.992 1.996 |      4.01      4.00 |       0       0
##     CARDS.WITH OCCUPATS       Between  10.41  10.42 | 2.503 2.517 |      4.16      4.14 |       0       0
##   FIGURES.WITH  ANIMALS       Between   8.88   8.88 | 2.003 2.001 |      4.43      4.44 |       0       0
##   FIGURES.WITH OCCUPATS       Between   9.59   9.63 | 2.155 2.149 |      4.45      4.48 |       0       0
##  MATRICES.WITH  ANIMALS       Between   8.04   8.05 | 1.862 1.843 |      4.32      4.37 |       0       0
##  MATRICES.WITH  FIGURES       Between   7.86   7.76 | 1.675 1.682 |      4.70      4.62 |       0       0
##  MATRICES.WITH OCCUPATS       Between   9.11   9.03 | 1.860 1.862 |      4.90      4.85 |       0       0
##          Means  ANIMALS       Between  30.15  30.15 | 0.472 0.472 |     63.86     63.92 |       0       0
##          Means    CARDS       Between  29.88  29.89 | 0.465 0.464 |     64.33     64.38 |       0       0
##          Means  FIGURES       Between  30.07  30.06 | 0.458 0.459 |     65.69     65.51 |       0       0
##          Means MATRICES       Between  29.73  29.73 | 0.439 0.439 |     67.64     67.68 |       0       0
##          Means OCCUPATS       Between  29.99  29.99 | 0.513 0.512 |     58.50     58.58 |       0       0
##          Means WORDLIST       Between  29.88  29.89 | 0.470 0.470 |     63.55     63.56 |       0       0
##      Variances  ANIMALS       Between  10.88  10.89 | 2.422 2.420 |      4.49      4.50 |       0       0
##      Variances    CARDS       Between  10.54  10.56 | 2.492 2.494 |      4.23      4.24 |       0       0
##      Variances  FIGURES       Between   9.88   9.91 | 2.050 2.051 |      4.82      4.83 |       0       0
##      Variances MATRICES       Between   9.11   9.10 | 1.762 1.763 |      5.17      5.16 |       0       0
##      Variances OCCUPATS       Between  13.41  13.42 | 2.765 2.764 |      4.85      4.86 |       0       0
##      Variances WORDLIST       Between  10.74  10.73 | 2.147 2.149 |      5.00      4.99 |       0       0
##  WORDLIST.WITH  ANIMALS       Between   9.58   9.72 | 2.119 2.106 |      4.52      4.61 |       0       0
##  WORDLIST.WITH    CARDS       Between   9.29   9.28 | 2.283 2.282 |      4.07      4.06 |       0       0
##  WORDLIST.WITH  FIGURES       Between   8.43   8.43 | 1.996 2.004 |      4.23      4.21 |       0       0
##  WORDLIST.WITH MATRICES       Between   8.53   8.51 | 1.773 1.775 |      4.81      4.79 |       0       0
##  WORDLIST.WITH OCCUPATS       Between  10.41  10.43 | 2.290 2.296 |      4.55      4.54 |       0       0
## 
## 
##   P-values that differ between models (p-value diff > 1e-04)
##   -----------------------------------
## None
## 
## 
##   Parameters unique to model 1: 21
##   -----------------------------
## 
##    paramHeader    param m1_est m1_se m1_est_se m1_pval BetweenWithin
##  WORDLIST.WITH    CARDS   9.86 1.527      6.46   0.000        Within
##  WORDLIST.WITH MATRICES   9.54 1.383      6.90   0.000        Within
##  WORDLIST.WITH  FIGURES   3.69 0.881      4.19   0.000        Within
##  WORDLIST.WITH  ANIMALS   4.61 1.143      4.03   0.000        Within
##  WORDLIST.WITH OCCUPATS   3.45 0.717      4.81   0.000        Within
##     CARDS.WITH MATRICES   9.62 1.147      8.39   0.000        Within
##     CARDS.WITH  FIGURES   2.87 0.880      3.25   0.001        Within
##     CARDS.WITH  ANIMALS   4.43 0.833      5.32   0.000        Within
##     CARDS.WITH OCCUPATS   3.37 0.755      4.46   0.000        Within
##  MATRICES.WITH  FIGURES   2.85 1.006      2.84   0.005        Within
##  MATRICES.WITH  ANIMALS   3.68 0.868      4.23   0.000        Within
##  MATRICES.WITH OCCUPATS   2.74 0.848      3.23   0.001        Within
##   FIGURES.WITH  ANIMALS   9.95 1.069      9.31   0.000        Within
##   FIGURES.WITH OCCUPATS   8.91 1.061      8.40   0.000        Within
##   ANIMALS.WITH OCCUPATS   8.90 0.863     10.32   0.000        Within
##      Variances WORDLIST  16.14 1.974      8.18   0.000        Within
##      Variances    CARDS  15.28 1.455     10.50   0.000        Within
##      Variances MATRICES  15.59 1.393     11.19   0.000        Within
##      Variances  FIGURES  16.54 1.247     13.26   0.000        Within
##      Variances  ANIMALS  15.10 1.181     12.78   0.000        Within
##      Variances OCCUPATS  13.25 1.040     12.74   0.000        Within
## 
##   Parameters unique to model 2: 15
##   -----------------------------
## 
##         paramHeader    param m2_est m2_se m2_est_se m2_pval BetweenWithin
##          NUMERIC.BY WORDLIST  3.151 0.302     10.41       0        Within
##          NUMERIC.BY    CARDS  3.153 0.231     13.64       0        Within
##          NUMERIC.BY MATRICES  3.029 0.217     13.98       0        Within
##          PERCEPT.BY  FIGURES  3.117 0.223     13.98       0        Within
##          PERCEPT.BY  ANIMALS  3.200 0.163     19.66       0        Within
##          PERCEPT.BY OCCUPATS  2.810 0.177     15.87       0        Within
##        PERCEPT.WITH  NUMERIC  0.379 0.050      7.62       0        Within
##  Residual.Variances WORDLIST  6.225 0.785      7.93       0        Within
##  Residual.Variances    CARDS  5.329 0.654      8.15       0        Within
##  Residual.Variances MATRICES  6.416 0.785      8.17       0        Within
##  Residual.Variances  FIGURES  6.816 0.773      8.82       0        Within
##  Residual.Variances  ANIMALS  4.859 0.633      7.68       0        Within
##  Residual.Variances OCCUPATS  5.358 0.712      7.53       0        Within
## 
## 
##   2 filtered from output (fixed and/or n.s.)
## 
##     Variances.NUMERIC.Within, Variances.PERCEPT.Within 
## 
## 
## ==============

6.6 Within two-factor, between one-factor

TITLE: Within two factor; Between one factor;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  numeric BY wordlist* cards matrices;
  percept BY figures* animals occupats;
  numeric@1 percept@1;
 %BETWEEN%
  general BY wordlist* cards matrices figures animals occupats;
  general@1;
OUTPUT: STDYX SAMPSTAT;
mobj6 <- update(mobj1, TITLE = ~"W/i two factor; B/w one factor;", MODEL = ~" %WITHIN%
    numeric BY wordlist* cards matrices;
    percept BY figures* animals occupats;
    numeric@1 percept@1;
   %BETWEEN%
    general BY wordlist* cards matrices figures animals occupats;
    general@1;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit6 <- mplusModeler(mobj6, modelout = "hox_m6.inp", run = 1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m6.inp
## Wrote data to: hox_m6.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m6.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m6.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m6.inp" 
## Reading model:  hox_m6.out
screenreg(mfit6, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type = "stdyx")
## 
## ===========================================
##                         Model 1            
## -------------------------------------------
## W WORDLIST<-NUMERIC         0.79 (0.04) ***
## W CARDS<-NUMERIC            0.80 (0.03) ***
## W MATRICES<-NUMERIC         0.77 (0.03) ***
## W FIGURES<-PERCEPT          0.77 (0.03) ***
## W ANIMALS<-PERCEPT          0.82 (0.02) ***
## W OCCUPATS<-PERCEPT         0.77 (0.03) ***
## B WORDLIST<-GENERAL         0.94 (0.03) ***
## B CARDS<-GENERAL            0.94 (0.04) ***
## B MATRICES<-GENERAL         0.88 (0.04) ***
## B FIGURES<-GENERAL          0.89 (0.04) ***
## B ANIMALS<-GENERAL          0.97 (0.03) ***
## B OCCUPATS<-GENERAL         0.94 (0.02) ***
## W PERCEPT<->NUMERIC         0.38 (0.05) ***
## B WORDLIST<-Intercepts      9.17 (0.98) ***
## B CARDS<-Intercepts         9.14 (1.09) ***
## B MATRICES<-Intercepts      9.97 (0.99) ***
## B FIGURES<-Intercepts       9.54 (1.00) ***
## B ANIMALS<-Intercepts       9.14 (1.05) ***
## B OCCUPATS<-Intercepts      8.20 (0.88) ***
## W NUMERIC<->NUMERIC         1.00 (0.00)    
## W PERCEPT<->PERCEPT         1.00 (0.00)    
## W WORDLIST<->WORDLIST       0.38 (0.06) ***
## W CARDS<->CARDS             0.35 (0.05) ***
## W MATRICES<->MATRICES       0.41 (0.05) ***
## W FIGURES<->FIGURES         0.41 (0.05) ***
## W ANIMALS<->ANIMALS         0.32 (0.04) ***
## W OCCUPATS<->OCCUPATS       0.40 (0.05) ***
## B GENERAL<->GENERAL         1.00 (0.00)    
## B WORDLIST<->WORDLIST       0.12 (0.06) *  
## B CARDS<->CARDS             0.12 (0.07)    
## B MATRICES<->MATRICES       0.22 (0.07) ** 
## B FIGURES<->FIGURES         0.21 (0.07) ** 
## B ANIMALS<->ANIMALS         0.06 (0.05)    
## B OCCUPATS<->OCCUPATS       0.12 (0.05) ** 
## -------------------------------------------
## AICC                    12868.03           
## CFI                         1.00           
## SRMR.Within                 0.02           
## SRMR.Between                0.02           
## ===========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.7 Within saturated, between independence

TITLE: Within saturated; Between independence; Table 14.7, p. 308;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  wordlist cards matrices figures animals occupats WITH
  wordlist cards matrices figures animals occupats;
 %BETWEEN%
OUTPUT: STDYX SAMPSTAT;
mobj7 <- update(
  mobj1, 
  TITLE = ~ "W/i saturated; B/w independence;", # Table 14.7, p. 308
  MODEL = ~ "
    %WITHIN%
      wordlist cards matrices figures animals occupats WITH
        wordlist cards matrices figures animals occupats;
    %BETWEEN%")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit7 <- mplusModeler(mobj7, modelout = "hox_m7.inp", run=1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m7.inp
## Wrote data to: hox_m7.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m7.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m7.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m7.inp" 
## Reading model:  hox_m7.out
## Warning in detectColumnNames(filename, modelSection, sectionType): Unable to determine column names for section model_results.
##   hox_m7.out
screenreg(mfit7, single.row=TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type="stdyx")
## 
## ==========================================
##                        Model 1            
## ------------------------------------------
## W WORDLIST<->CARDS         0.74 (0.03) ***
## W WORDLIST<->MATRICES      0.72 (0.03) ***
## W WORDLIST<->FIGURES       0.45 (0.05) ***
## W WORDLIST<->ANIMALS       0.52 (0.05) ***
## W WORDLIST<->OCCUPATS      0.48 (0.06) ***
## W CARDS<->MATRICES         0.73 (0.03) ***
## W CARDS<->FIGURES          0.42 (0.06) ***
## W CARDS<->ANIMALS          0.53 (0.04) ***
## W CARDS<->OCCUPATS         0.49 (0.06) ***
## W MATRICES<->FIGURES       0.40 (0.06) ***
## W MATRICES<->ANIMALS       0.46 (0.04) ***
## W MATRICES<->OCCUPATS      0.43 (0.06) ***
## W FIGURES<->ANIMALS        0.74 (0.03) ***
## W FIGURES<->OCCUPATS       0.72 (0.03) ***
## W ANIMALS<->OCCUPATS       0.76 (0.02) ***
## B WORDLIST<-Means         24.30 (4.00) ***
## B CARDS<-Means            22.75 (3.93) ***
## B MATRICES<-Means         20.66 (2.96) ***
## B FIGURES<-Means          20.38 (2.97) ***
## B ANIMALS<-Means          29.28 (9.62) ** 
## B OCCUPATS<-Means         17.77 (2.66) ***
## W WORDLIST<->WORDLIST      1.00 (0.00)    
## W CARDS<->CARDS            1.00 (0.00)    
## W MATRICES<->MATRICES      1.00 (0.00)    
## W FIGURES<->FIGURES        1.00 (0.00)    
## W ANIMALS<->ANIMALS        1.00 (0.00)    
## W OCCUPATS<->OCCUPATS      1.00 (0.00)    
## B WORDLIST<->WORDLIST      1.00 (0.00)    
## B CARDS<->CARDS            1.00 (0.00)    
## B MATRICES<->MATRICES      1.00 (0.00)    
## B FIGURES<->FIGURES        1.00 (0.00)    
## B ANIMALS<->ANIMALS        1.00 (0.00)    
## B OCCUPATS<->OCCUPATS      1.00 (0.00)    
## ------------------------------------------
## AICC                   13029.19           
## CFI                        0.83           
## SRMR.Within                0.29           
## SRMR.Between               0.76           
## ==========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.8 Within independence, between saturated

TITLE: Within independence--Between saturated; Table 14.7, p. 308;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
 %BETWEEN%
  wordlist cards matrices figures animals occupats WITH
  wordlist cards matrices figures animals occupats; 
OUTPUT: STDYX SAMPSTAT;
mobj8 <- update(
  mobj1, 
  TITLE = ~ "W/i independence; B/w saturated;", # Table 14.7, p. 308
  MODEL = ~ "
    %WITHIN%
    %BETWEEN%
      wordlist cards matrices figures animals occupats WITH
        wordlist cards matrices figures animals occupats;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit8 <- mplusModeler(mobj8, modelout = "hox_m8.inp", run=1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m8.inp
## Wrote data to: hox_m8.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m8.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m8.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m8.inp" 
## Reading model:  hox_m8.out
## Warning in detectColumnNames(filename, modelSection, sectionType): Unable to determine column names for section model_results.
##   hox_m8.out
screenreg(mfit8, single.row=TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type="stdyx")
## 
## ==========================================
##                        Model 1            
## ------------------------------------------
## B WORDLIST<->CARDS         0.98 (0.04) ***
## B WORDLIST<->MATRICES      0.97 (0.04) ***
## B WORDLIST<->FIGURES       0.85 (0.08) ***
## B WORDLIST<->ANIMALS       0.91 (0.09) ***
## B WORDLIST<->OCCUPATS      0.90 (0.06) ***
## B CARDS<->MATRICES         0.96 (0.06) ***
## B CARDS<->FIGURES          0.89 (0.07) ***
## B CARDS<->ANIMALS          0.94 (0.07) ***
## B CARDS<->OCCUPATS         0.92 (0.05) ***
## B MATRICES<->FIGURES       0.85 (0.07) ***
## B MATRICES<->ANIMALS       0.87 (0.08) ***
## B MATRICES<->OCCUPATS      0.85 (0.07) ***
## B FIGURES<->ANIMALS        0.97 (0.08) ***
## B FIGURES<->OCCUPATS       0.95 (0.04) ***
## B ANIMALS<->OCCUPATS       0.99 (0.07) ***
## B WORDLIST<-Means          9.02 (0.98) ***
## B CARDS<-Means             9.12 (1.12) ***
## B MATRICES<-Means          9.71 (0.95) ***
## B FIGURES<-Means           9.39 (1.01) ***
## B ANIMALS<-Means           8.91 (1.12) ***
## B OCCUPATS<-Means          8.18 (0.91) ***
## W WORDLIST<->WORDLIST      1.00 (0.00)    
## W CARDS<->CARDS            1.00 (0.00)    
## W MATRICES<->MATRICES      1.00 (0.00)    
## W FIGURES<->FIGURES        1.00 (0.00)    
## W ANIMALS<->ANIMALS        1.00 (0.00)    
## W OCCUPATS<->OCCUPATS      1.00 (0.00)    
## B WORDLIST<->WORDLIST      1.00 (0.00)    
## B CARDS<->CARDS            1.00 (0.00)    
## B MATRICES<->MATRICES      1.00 (0.00)    
## B FIGURES<->FIGURES        1.00 (0.00)    
## B ANIMALS<->ANIMALS        1.00 (0.00)    
## B OCCUPATS<->OCCUPATS      1.00 (0.00)    
## ------------------------------------------
## AICC                   13673.39           
## CFI                        0.27           
## SRMR.Within                0.36           
## SRMR.Between               0.06           
## ==========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.9 Within independence, between one-factor

TITLE: Within independence--Between one factor
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
 %BETWEEN%
  general BY wordlist* cards matrices figures animals occupats;
  general@1;
OUTPUT: STDYX SAMPSTAT;
mobj9 <- update(mobj1, TITLE = ~"Within independence; Between one factor", MODEL = ~"
    %WITHIN%
    %BETWEEN%
      general BY wordlist* cards matrices figures animals occupats;
      general@1;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit9 <- mplusModeler(mobj9, modelout = "hox_m9.inp", run = 1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m9.inp
## Wrote data to: hox_m9.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m9.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m9.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m9.inp" 
## Reading model:  hox_m9.out
## Warning in rbind(character(0), c("WORDLIST", "0.939", "0.069", "13.619", : number of columns of result is not a multiple of vector length (arg 6)
screenreg(mfit9, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type = "stdyx")
## 
## =============================================
##                         Model 1              
## ---------------------------------------------
## B WORDLIST<-GENERAL         0.97   (0.04) ***
## B CARDS<-GENERAL            0.98   (0.03) ***
## B MATRICES<-GENERAL         0.94   (0.04) ***
## B FIGURES<-GENERAL          0.95   (0.03) ***
## B ANIMALS<-GENERAL          1.00   (0.04) ***
## B OCCUPATS<-GENERAL         0.97   (0.03) ***
## B WORDLIST<-Intercepts      9.07   (0.96) ***
## B CARDS<-Intercepts         9.16   (1.11) ***
## B MATRICES<-Intercepts      9.80   (0.97) ***
## B FIGURES<-Intercepts       9.44   (0.98) ***
## B ANIMALS<-Intercepts       9.04   (1.03) ***
## B OCCUPATS<-Intercepts      8.21   (0.89) ***
## W WORDLIST<->WORDLIST       1.00   (0.00)    
## W CARDS<->CARDS             1.00   (0.00)    
## W MATRICES<->MATRICES       1.00   (0.00)    
## W FIGURES<->FIGURES         1.00   (0.00)    
## W ANIMALS<->ANIMALS         1.00   (0.00)    
## W OCCUPATS<->OCCUPATS       1.00   (0.00)    
## B GENERAL<->GENERAL         1.00   (0.00)    
## B WORDLIST<->WORDLIST       0.06   (0.07)    
## B CARDS<->CARDS             0.03   (0.07)    
## B MATRICES<->MATRICES       0.13   (0.07)    
## B FIGURES<->FIGURES         0.10   (0.06)    
## B ANIMALS<->ANIMALS        -0.00 (999.00)    
## B OCCUPATS<->OCCUPATS       0.06   (0.05)    
## ---------------------------------------------
## AICC                    13671.30             
## CFI                         0.30             
## SRMR.Within                 0.36             
## SRMR.Between                0.06             
## =============================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.10 Within two-factor, between independence

TITLE: Within two factor--Between independence; Table 14.1, p. 301;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  numeric BY wordlist* cards matrices;
  percept BY figures* animals occupats;
  numeric@1 percept@1;
 %BETWEEN%
OUTPUT: STDYX SAMPSTAT;
mobj10 <- update(
  mobj1, 
  TITLE = ~ "W/i two factor; B/w independence;", # Table 14.1, p. 301
  MODEL = ~ "
    %WITHIN%
      numeric BY wordlist* cards matrices;
      percept BY figures* animals occupats;
      numeric@1 percept@1;
    %BETWEEN%")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit10 <- mplusModeler(mobj10, modelout = "hox_m10.inp", run=1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m10.inp
## Wrote data to: hox_m10.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m10.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m10.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m10.inp" 
## Reading model:  hox_m10.out
screenreg(mfit10, single.row=TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type="stdyx")
## 
## ==========================================
##                        Model 1            
## ------------------------------------------
## W WORDLIST<-NUMERIC        0.86 (0.02) ***
## W CARDS<-NUMERIC           0.87 (0.02) ***
## W MATRICES<-NUMERIC        0.83 (0.02) ***
## W FIGURES<-PERCEPT         0.83 (0.02) ***
## W ANIMALS<-PERCEPT         0.89 (0.02) ***
## W OCCUPATS<-PERCEPT        0.86 (0.02) ***
## W PERCEPT<->NUMERIC        0.64 (0.04) ***
## B WORDLIST<-Means         24.05 (3.86) ***
## B CARDS<-Means            22.27 (3.65) ***
## B MATRICES<-Means         20.86 (3.09) ***
## B FIGURES<-Means          20.50 (3.04) ***
## B ANIMALS<-Means          28.19 (9.07) ** 
## B OCCUPATS<-Means         17.80 (2.68) ***
## W NUMERIC<->NUMERIC        1.00 (0.00)    
## W PERCEPT<->PERCEPT        1.00 (0.00)    
## W WORDLIST<->WORDLIST      0.26 (0.04) ***
## W CARDS<->CARDS            0.24 (0.03) ***
## W MATRICES<->MATRICES      0.31 (0.04) ***
## W FIGURES<->FIGURES        0.30 (0.04) ***
## W ANIMALS<->ANIMALS        0.20 (0.03) ***
## W OCCUPATS<->OCCUPATS      0.27 (0.04) ***
## B WORDLIST<->WORDLIST      1.00 (0.00)    
## B CARDS<->CARDS            1.00 (0.00)    
## B MATRICES<->MATRICES      1.00 (0.00)    
## B FIGURES<->FIGURES        1.00 (0.00)    
## B ANIMALS<->ANIMALS        1.00 (0.00)    
## B OCCUPATS<->OCCUPATS      1.00 (0.00)    
## ------------------------------------------
## AICC                   13018.93           
## CFI                        0.85           
## SRMR.Within                0.29           
## SRMR.Between               0.76           
## ==========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.11 Within two-factor, between null

TITLE: Within two factor--Between null; Table 14.1, p. 301;
DATA: FILE IS "FamIQData.dat";
VARIABLE: NAMES ARE
 family child wordlist cards matrices figures animals occupats;
USEVARIABLES ARE wordlist cards matrices figures animals occupats;
CLUSTER IS family;
ANALYSIS: TYPE IS TWOLEVEL;
 ESTIMATOR IS MLR;
MODEL:
 %WITHIN%
  numeric BY wordlist* cards matrices;
  percept BY figures* animals occupats;
  numeric-percept@1;
 %BETWEEN%
 wordlist@0 cards@0 matrices@0 figures@0 animals@0 occupats@0; 
OUTPUT: STDYX SAMPSTAT;
mobj11a <- update(
  mobj1, 
  TITLE = ~ "Within two factor; Between null;", #Table 14.1, p. 301
  MODEL = ~ "
    %WITHIN%
      numeric BY wordlist* cards matrices;
      percept BY figures* animals occupats;
      numeric@1 percept@1;
    %BETWEEN%
      wordlist@0 cards@0 matrices@0 figures@0 animals@0 occupats@0;")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit11a <- mplusModeler(mobj11a, modelout = "hox_m11a.inp", run=1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m11a.inp
## Wrote data to: hox_m11a.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m11a.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m11a.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m11a.inp" 
## Reading model:  hox_m11a.out
screenreg(mfit11a, single.row=TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type="stdyx")
## 
## ===========================================
##                        Model 1             
## -------------------------------------------
## W WORDLIST<-NUMERIC        0.85  (0.02) ***
## W CARDS<-NUMERIC           0.85  (0.02) ***
## W MATRICES<-NUMERIC        0.81  (0.02) ***
## W FIGURES<-PERCEPT         0.80  (0.02) ***
## W ANIMALS<-PERCEPT         0.89  (0.02) ***
## W OCCUPATS<-PERCEPT        0.85  (0.02) ***
## W PERCEPT<->NUMERIC        0.67  (0.04) ***
## B WORDLIST<-Means       2994.25 (44.43) ***
## B CARDS<-Means          2984.00 (44.36) ***
## B MATRICES<-Means       2973.50 (43.21) ***
## B FIGURES<-Means        3008.50 (45.05) ***
## B ANIMALS<-Means        3011.75 (44.48) ***
## B OCCUPATS<-Means       3003.00 (47.50) ***
## W NUMERIC<->NUMERIC        1.00  (0.00)    
## W PERCEPT<->PERCEPT        1.00  (0.00)    
## W WORDLIST<->WORDLIST      0.28  (0.04) ***
## W CARDS<->CARDS            0.27  (0.04) ***
## W MATRICES<->MATRICES      0.34  (0.04) ***
## W FIGURES<->FIGURES        0.35  (0.04) ***
## W ANIMALS<->ANIMALS        0.22  (0.03) ***
## W OCCUPATS<->OCCUPATS      0.27  (0.04) ***
## B WORDLIST<->WORDLIST      1.00  (0.00)    
## B CARDS<->CARDS            1.00  (0.00)    
## B MATRICES<->MATRICES      1.00  (0.00)    
## B FIGURES<->FIGURES        1.00  (0.00)    
## B ANIMALS<->ANIMALS        1.00  (0.00)    
## B OCCUPATS<->OCCUPATS      1.00  (0.00)    
## -------------------------------------------
## AICC                   13152.12            
## CFI                        0.58            
## SRMR.Within                0.40            
## SRMR.Between               0.79            
## ===========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

Different assumption in which we force variance to be at level 1 only, rather than allowing for it to be at L1 and L2, but constraining L2 to zero.

mobj11b <- update(mobj1, TITLE = ~"Within two factor--Between null; Table 14.1, p. 301;", VARIABLE = ~. + "WITHIN ARE wordlist cards matrices figures animals occupats;", MODEL = ~"
    %WITHIN%
      numeric BY wordlist* cards matrices;
      percept BY figures* animals occupats;
      numeric@1 percept@1;
    %BETWEEN%")
## R variables selected automatically as any variable name that
## occurs in the MODEL, VARIABLE, or DEFINE section.
mfit11b <- mplusModeler(mobj11b, modelout = "hox_m11b.inp", run = 1L, Mplus_command = mbin, hashfilename = FALSE)
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## Wrote model to: hox_m11b.inp
## Wrote data to: hox_m11b.dat
## Warning in prepareMplusData(df = data[i, , drop = FALSE], keepCols = object$usevariables, : The file(s)
##  'hox_m11b.dat' 
## currently exist(s) and will be overwritten
## 
## Running model: hox_m11b.inp 
## System command: cd "." && "/Users/mnh5174/Applications/Mplus/mplus" "hox_m11b.inp" 
## Reading model:  hox_m11b.out
screenreg(mfit11b, single.row = TRUE, summaries = c("AICC", "CFI", "SRMR.Within", "SRMR.Between"), type = "stdyx")
## 
## ===========================================
##                         Model 1            
## -------------------------------------------
## W WORDLIST<-NUMERIC         0.85 (0.02) ***
## W CARDS<-NUMERIC            0.85 (0.02) ***
## W MATRICES<-NUMERIC         0.81 (0.02) ***
## W FIGURES<-PERCEPT          0.80 (0.02) ***
## W ANIMALS<-PERCEPT          0.89 (0.02) ***
## W OCCUPATS<-PERCEPT         0.85 (0.02) ***
## W PERCEPT<->NUMERIC         0.67 (0.04) ***
## W WORDLIST<-Intercepts      5.86 (0.36) ***
## W CARDS<-Intercepts         5.95 (0.32) ***
## W MATRICES<-Intercepts      6.03 (0.27) ***
## W FIGURES<-Intercepts       5.89 (0.25) ***
## W ANIMALS<-Intercepts       6.00 (0.31) ***
## W OCCUPATS<-Intercepts      5.94 (0.32) ***
## W NUMERIC<->NUMERIC         1.00 (0.00)    
## W PERCEPT<->PERCEPT         1.00 (0.00)    
## W WORDLIST<->WORDLIST       0.28 (0.04) ***
## W CARDS<->CARDS             0.27 (0.04) ***
## W MATRICES<->MATRICES       0.34 (0.04) ***
## W FIGURES<->FIGURES         0.35 (0.04) ***
## W ANIMALS<->ANIMALS         0.22 (0.03) ***
## W OCCUPATS<->OCCUPATS       0.27 (0.04) ***
## -------------------------------------------
## AICC                    13152.15           
## CFI                         1.00           
## SRMR.Within                 0.01           
## SRMR.Between                0.00           
## ===========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

6.12 Model selection

SummaryTable(list(mfit1, mfit2, mfit3, mfit4, mfit5, mfit6, mfit7, mfit8, mfit9, mfit10, mfit11a), keepCols = c("Title", "LL", "Parameters", "AICC", "BIC", "RMSEA_Estimate", "CFI", 
    "SRMR.Within", "SRMR.Between"), sortBy = "AICC")
##                                          Title    LL Parameters  AICC   BIC RMSEA_Estimate   CFI SRMR.Within SRMR.Between
## 1              W/i two factor; B/w one factor; -6400         31 12868 12986          0.000 1.000       0.022        0.018
## 2               W/i saturated; B/w one factor; -6397         39 12880 13027          0.000 1.000       0.003        0.017
## 3        Within two factor--Between saturated; -6398         40 12885 13035          0.000 1.000       0.022        0.005
## 4                W/i saturated, B/w saturated; -6394         48 12898 13076          0.000 1.000       0.000        0.000
## 5            W/i two factor; B/w independence; -6483         25 13019 13115          0.141 0.848       0.292        0.758
## 6             W/i saturated; B/w independence; -6479         33 13029 13155          0.183 0.833       0.290        0.758
## 7             Within two factor; Between null; -6556         19 13152 13226          0.210 0.577       0.400        0.794
## 8              W/i one factor; B/w one factor; -6562         30 13188 13303          0.191 0.782       0.186        0.051
## 9      Within independence; Between one factor -6810         24 13671 13764          0.297 0.300       0.365        0.063
## 10            W/i independence; B/w saturated; -6801         33 13673 13799          0.383 0.272       0.365        0.060
## 11  Within independence; Between independence; -6985         18 14007 14077          0.317 0.000       0.365        0.638