SEM Example 07 - Comparing ML and WLSMV Estimators for Ordinal, Likert-Type Items: Demonstrations with the Health Behaviour in School-Aged Children Dataset

Chong Xing, Center for Research Methods and Data Analysis, University of Kansas <cxing@ku.edu>

Paul Johnson, Center for Research Methods and Data Analysis, University of Kansas <pauljohn@ku.edu>

Meghan Sullivan, Center for Research Methods and Data Analysis, University of Kansas <meg.sullivan@ku.edu>

Please visit http://crmda.ku.edu/guides

Keywords: Structural Equation Modeling, Ordinal Data, Weighted Least Squares Estimators, Maximum Likelihood, HBSC"

Abstract

Likert scale data can be treated as numbers, say 1 through 5, or they can be treated as ordinal labels, say “Strongly Disagree” to “Strongly Agree”. Methods to estimate the ordinal version have been available for some time, and yet it is still common to treat the ordinal data as if it were numeric. Here we explore the question by presenting examples of structural equation models estimated in both ways. The Maximum Likelihood (ML) method is used when the data are treated as numbers, while Weighted Least Squares (WLSMV) estimator is employed when the data is treated as categorical. Example data was obtained from the Health Behavior in School-Aged Children study (HBSC; Iannotti, 2005-2006). Our data set includes students in grades 6 and 7 (N = 4,284).

Data Importation

A subset of the HBSC data, hbsc-subset2.dat, was created and prepared for usage in Mplus. The Mplus exercises that investigate this data can be found in Mplus.

Here we import the same data file in R, assign column names, and then use the kutils variable key to create additional categorical variables.

hbsc.colnames <- scan(file.path(ddir, "hbsc-subset2-colnames.txt"), what = "character")
hbsc.orig <- read.table(file.path(ddir, "hbsc-subset2.dat"),
                        header = FALSE, col.names = hbsc.colnames,
                        na.strings = c("-999", ".", "NA"))
key <- keyImport("hbsc-subset2-key2.csv")
keyImport guessed that is a wide format key.
hbsc <- keyApply(hbsc.orig, key, drop = FALSE, ignoreCase = FALSE, debug = FALSE,
                 diagnostic = FALSE)

The file named hbsc-subset2.dat is that “raw” text data file that was used in our companion Mplus exercises.

We restrict our attention to the students in the 6th and 7th grades.

hbsc <- hbsc[!is.na(hbsc$Grade) & hbsc$Grade %in% c("6", "7"), ]

The variables in our R analysis have the same names that are used in the Mplus files. However, to demonstrate the different methods available in R for dealing with ordered categorical variables, we insert ordinal versions of the same variables with the letters “_o" appended.

Ordinal.Versions Integer.Versions
body1r_o body1r
body2_o body2
body3r_o body3r
body4_o body4
body5r_o body5r
phyhlth1_o phyhlth1
phyhlth2_o phyhlth2
phyhlth3_o phyhlth3
phyhlth4_o phyhlth4
phyhlth5_o phyhlth5
phyhlth6_o phyhlth6
phyhlth7_o phyhlth7
phyhlth8_o phyhlth8
depress1_o depress1
depress2_o depress2
depress3_o depress3
depress4_o depress4
depress5_o depress5
depress6_o depress6
bullied1_o bullied1
bullied2_o bullied2
bullied3_o bullied3
bullied4_o bullied4
bullied5_o bullied5
bullied6_o bullied6
bullied7_o bullied7
bullied8_o bullied8
bullied9_o bullied9
bullier1_o bullier1
bullier2_o bullier2
bullier3_o bullier3
bullier4_o bullier4
bullier5_o bullier5
bullier6_o bullier6
bullier7_o bullier7
bullier8_o bullier8
bullier9_o bullier9
alc1_o alc1
alc2_o alc2
alc3_o alc3
alc4_o alc4
alc5_o alc5

The following demonstrates that the R summary function is somewhat befuddled by the categorical variables. When we provide the same information in 2 variables, one as an integer and one as a factor, the numeric version is displayed quite differently (Mean and Median), while the categorical varsion is treated as tabular information.

vars <- c("Grade", "depress1", "depress2", "depress3", "depress1_o", "depress2_o", "depress3_o")
summary(hbsc[ , vars])
     Grade          depress1        depress2        depress3    
 Min.   :6.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
 1st Qu.:6.000   1st Qu.:1.000   1st Qu.:2.000   1st Qu.:1.000  
 Median :6.000   Median :2.000   Median :3.000   Median :1.000  
 Mean   :6.439   Mean   :2.332   Mean   :2.684   Mean   :1.817  
 3rd Qu.:7.000   3rd Qu.:3.000   3rd Qu.:3.000   3rd Qu.:3.000  
 Max.   :7.000   Max.   :5.000   Max.   :5.000   Max.   :5.000  
                 NA's   :78      NA's   :87      NA's   :96     
     depress1_o       depress2_o       depress3_o  
 Never    :1212   Never    : 800   Never    :2558  
 Seldom   :1113   Seldom   : 977   Seldom   : 566  
 Sometimes:1307   Sometimes:1449   Sometimes: 558  
 Often    : 420   Often    : 693   Often    : 286  
 Always   : 154   Always   : 278   Always   : 220  
 NA's     :  78   NA's     :  87   NA's     :  96  
                                                   

The values represented depress1 and depress1_o are the same. The difference is that most R functions will treat depress1 as an integer (or numeric) score, while the ordinal variant is an R factor. For example:

                  depress1 as integer
depress1 as factor    1    2    3    4    5
         Never     1212    0    0    0    0
         Seldom       0 1113    0    0    0
         Sometimes    0    0 1307    0    0
         Often        0    0    0  420    0
         Always       0    0    0    0  154

Complication: Missing Values (methodological moving target)

Historically, cases in which there were missing values were subjected to the listwise deletion. That was the standard through the 1970s and 1980s.

In the 1970s, would-be SEM practitioners had no choice. They were forced to treat categorical “Likert-type” variables as if the scores were numbers. Although there were several different estimators considered during that time, lets refer to them as “ML of covariance structures,” or ML for short.

Early in the 1980s, a chain of technical developments results in new methods that explicitly treated the ordinal variables as categorical methods. Again, there were several similar implementations with slightly different terminology, but lets call these “Weighted Least Squares” methods, or WLS for short. Those methods were developed with the assumption of listwise deletion of missing data.

There were efforts to find methods that avoided deletion of cases on which there were missing values. The most important analytical breakthrough was full information maximum likelihood, a method suitable for numeric indicator variables. Colloquilly referred to as “FIML,” this technique was originally developed for factor analysis. Its technique was readily “grafted onto” estimators for confirmatory factor analysis and structural equation models. FIML is regarded as an analytic breakthrough and there are many studies which suggest it should be preferred to listwise deletion. It is possibly as good, or even more desirable than, multiple imputation for missing values.

FIML requires individual-level data. In the days before FIML, we could exchange covariance matrices and conduct replication analysis because the covariance structure contained all of the information necessary to conduct SEM calculations. That is not true of FIML, which uses the individual scores to work its magic. In a real sense, the adoption of FIML was a methodological watershed. It was not meaningful, anymore, to refer to the field as “the analysis of covariance structures” (as it was known in the 1990s).

In the time when SEM was thought of as the analysis of covariance structures, methods of dealing with missing values in the construction of the observed covariance matrix were considered. The “pairwise complete” method does what its name implies, it takes variables in pairs and uses the available observations. There were mixed results about this method, some suggested it was less desirable even than listwise deletion of missing values.

Now the bad news: FIML was not applicable to SEM estimators for ordinal data. FIML succeeded in large part because the indicators were assumed to be draws from a multivariate normal distribution. Since that assumption does not apply when the observed variables are discrete (ordinal) or non-normal, the researcher was forced into an unpleasant choice. The options seemed to be

  1. Assert that the data is actually numeric (normal), in order to use FIML, or
  2. Apply the ordinal estimator, which uses listwise deletion of cases with missing values.

Those were the options, at least until very recently. This is the point at which we hit a confusing problem. In order to compare results between Mplus and R, it is necessary to keep track of at least three key terms

  1. The estimator: ML, or WLS, or others
  2. The treatment of missing values: listwise deletion, pairwise-complete correlation matrix
  3. Corrections for standard errors and test diagnostics, often referred to as “robust” methods.

In order to align the results from R’s lavaan package and Mplus, it is necessary for us to make sure all 3 elements are in line. Estimates can also be returned in different scalings, depending on the parameters.

Once this is understood, we hit a problem. Mplus implements some estimators that are not available in R.

We are able to align the R and Mplus results for data that is treated as if it were numeric. The estimator full information maximum likelihood with robust standard errors (MLR) offers the same results in either. Both use FIML.

However, the Mplus estimators for categorical variables have changed over the years. Mplus WLSMV is no longer using listwise deletion of missing values. It now uses a “pairwise complete” correlation matrix. We were “bitten” by the fact that WLSMV estimates from R’s lavaan and Mplus are no longer equivalent. In order to align the lavaan estimates for categorical data using WLSMV with Mplus, we need to run the Mplus estimator with the argument LISTWISE = ON;. At the current time, lavaan uses listwise deletion by default, but it is not harmful to explicitly state missing = "listwise"

However, it now becomes tricky for us to compare the MLR estimates for the numeric version with the categorical estimates because they are based on samples of different sizes. We made it possible to align the categorical data parameter estimates between Mplus and lavaan, but within lavaan, it is more difficult to compare the numeric ML with categorical WLSMV because the sample sizes differ. We ought to compare WLSMV-listwise with ML-listwise, that is to say.

Putting that issue aside, in Mplus 7, a new maximum likelihood estimator for the categorical data model was introduced. This is an ML estimator for categorical indicators. It uses the estimator familiar to researchers in education, the marginal maximum likelihood estimator that had risen to prominence in item response theory (IRT). The IRT-based ML estimator in Mplus 7 is very time consuming and, at the current time, the lavaan implementation of that estimator is under development.

SECTION-I: One-Factor CFA Models

This section presents three single-factor CFA analyses of depression, alcohol use, and bullying victimization. We provide two sets of estimates, one that treats the 1-2-3-4-5 scores as if they were numeric (MLR, Maximum Likelihood with Robust standard errors) and ordinal (WLSMV, Weighted Least Squares with Mean and Variance adjustment). The WLSMV estimates employ listwise deletion of missing values.

One-Factor CFA for Depression

Depression was measured by six items with five-point Likert scales. Items were framed with the following stem:

"Think about how you have been feeling over the last 30 days.  Mark
the number that goes with how often you have felt or done each of
these." 

Example items include “Were you very sad?” and “Did you feel hopeless about the future?”

Response options were 1) Never, 2) Seldom, 3) Sometimes, 4) Often, and 5) Always.

Method 1 - ML estimates treat data as numeric

The first model is estimated with lavaan MLR.

A few highlights in our model code:

  • Fixed-factor identification sets the factor variance equal to 1. This was done to set the scale of the latent variable.
  • NA* specifies that the item’s factor loading should be freely estimated.
  • depress ~~ 1*depress sets/fixes the variance of the latent variable, depress, to 1.
  • Alternatively, std.lv = TRUE can be included in the cfa() function to achieve the same purpose.
cfa.depress <- '
depress =~ NA*depress1 + depress2 + depress3
            + depress4 + depress5 + depress6
depress ~~ 1*depress
'
cfa.depress.mlr.listwise <-
  cfa(model = cfa.depress, data = hbsc,
      mimic = "Mplus", estimator = "MLR", meanstructure = TRUE, 
      missing = "listwise")
  • meanstructure: The meanstructure parameter is included because there is a quirk in the lavaan output. With missing = "listwise", the Intercepts are not included in the output, but with FIML, they are.
  • std.ov: The output from WLSMV will provide coefficients that are scaled. The most similar scaling for the MLR model is for standardized observed variables. (Note in the output the columns for the usual and standardized coefficients will be the same because of this setting.)

Lavaan’s summary function does not create an output object. It prints direct to the screen. The style of the output is similar to another well-known SEM software program.

summary(cfa.depress.mlr.listwise, fit.measures = TRUE, standardized = TRUE)
lavaan 0.6-3 ended normally after 17 iterations

  Optimization method                           NLMINB
  Number of free parameters                         18

                                                  Used       Total
  Number of observations                          4103        4284

  Estimator                                         ML      Robust
  Model Fit Test Statistic                     257.596     195.851
  Degrees of freedom                                 9           9
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  1.315
    for the Yuan-Bentler correction (Mplus variant)

Model test baseline model:

  Minimum Function Test Statistic             5718.519    4373.368
  Degrees of freedom                                15          15
  P-value                                        0.000       0.000

User model versus baseline model:

  Comparative Fit Index (CFI)                    0.956       0.957
  Tucker-Lewis Index (TLI)                       0.927       0.929

  Robust Comparative Fit Index (CFI)                         0.957
  Robust Tucker-Lewis Index (TLI)                            0.928

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)             -37612.835  -37612.835
  Scaling correction factor                                  1.000
    for the MLR correction
  Loglikelihood unrestricted model (H1)     -37484.037  -37484.037
  Scaling correction factor                                  1.105
    for the MLR correction

  Number of free parameters                         18          18
  Akaike (AIC)                               75261.670   75261.670
  Bayesian (BIC)                             75375.420   75375.420
  Sample-size adjusted Bayesian (BIC)        75318.224   75318.224

Root Mean Square Error of Approximation:

  RMSEA                                          0.082       0.071
  90 Percent Confidence Interval          0.074  0.091       0.064  0.079
  P-value RMSEA <= 0.05                          0.000       0.000

  Robust RMSEA                                               0.082
  90 Percent Confidence Interval                             0.072  0.092

Standardized Root Mean Square Residual:

  SRMR                                           0.030       0.030

Parameter Estimates:

  Information                                 Observed
  Observed information based on                Hessian
  Standard Errors                   Robust.huber.white

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  depress =~                                                            
    depress1          0.737    0.018   41.208    0.000    0.737    0.668
    depress2          0.709    0.018   38.771    0.000    0.709    0.616
    depress3          0.790    0.021   37.152    0.000    0.790    0.656
    depress4          0.864    0.021   42.079    0.000    0.864    0.652
    depress5          0.748    0.023   32.912    0.000    0.748    0.543
    depress6          0.733    0.022   32.907    0.000    0.733    0.547

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .depress1          2.335    0.017  135.555    0.000    2.335    2.116
   .depress2          2.681    0.018  149.184    0.000    2.681    2.329
   .depress3          1.810    0.019   96.308    0.000    1.810    1.504
   .depress4          2.204    0.021  106.615    0.000    2.204    1.664
   .depress5          2.476    0.022  115.058    0.000    2.476    1.796
   .depress6          2.466    0.021  117.910    0.000    2.466    1.841
    depress           0.000                               0.000    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    depress           1.000                               1.000    1.000
   .depress1          0.674    0.021   32.684    0.000    0.674    0.554
   .depress2          0.822    0.023   36.527    0.000    0.822    0.621
   .depress3          0.825    0.026   31.370    0.000    0.825    0.569
   .depress4          1.007    0.030   33.127    0.000    1.007    0.574
   .depress5          1.341    0.034   39.011    0.000    1.341    0.705
   .depress6          1.258    0.033   38.212    0.000    1.258    0.701

After this point, we will not include the voluminous output within our presentation here. Instead, we will divert the output into files that can be reviewed in separate files. For this summary, the lavaan output is saved in the file: output/cfa.depress.mlr.listwise.Rout

The Mplus output for the same example is saved in the file: cfa-depress-mlr-listwise.out

If the listwise option is not specified, ML uses the so-called FIML estimator, which does not result in the omission of as many cases in this CFA. Note the Warning, however. Some cases are missing on all of the variables, so they are removed from the analysis. In order to be kept in the calculation, an observation must have at least one observed value on an indicator (and there must be no missings on exogenous predictor variables, but that is not a concern here).

cfa.depress.mlr.fiml <- cfa(model = cfa.depress, data = hbsc,
                            mimic = "Mplus", estimator = "MLR",
                            meanstructure = TRUE)
Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some cases are empty and will be ignored:
  40 47 121 179 189 233 297 315 325 362 506 526 558 663 779 846 847 1006 1177 1200 1204 1222 1228 1241 1530 1533 1597 1614 1662 1671 1688 1856 1925 1961 2011 2039 2087 2203 2238 2264 2272 2282 2324 2381 2385 2390 2605 2750 2809 2867 3015 3391 3398 3551 3673 3724 3763 3875 3938 3985 4094 4198 4201
summary(cfa.depress.mlr.fiml, fit.measures = TRUE, standardized = TRUE)

The lavaan output is saved in the file: output/cfa.depress.mlr.fiml.Rout

The Mplus output for the same example is saved in the file: cfa-depress-mlr-fiml.out

Method 2 - The ordinal “threshold” model with WLS

A casual inspection of the variables indicates that the 5 point scales are quite skewed. These might be the “right tail” of a normal variable that has been split into categorized sections.

Because our variables that are named with “_o" are declared as ordered factors, the fitting function in lavaan will notice that they are not numeric and estimate them accordingly. We rewrite the lavaan formula as follows

cfa.depress.ord <- '
depress =~ NA*depress1_o + depress2_o + depress3_o
            + depress4_o + depress5_o + depress6_o
depress ~~ 1*depress
'
cfa.depress.ord.wlsmv <-
    cfa(model = cfa.depress.ord, data = hbsc, mimic = "Mplus",
        estimator = "WLSMV", missing = "listwise")
summary(cfa.depress.ord.wlsmv, fit.measures = TRUE, standardized = TRUE)

The lavaan output is saved in the file: output/cfa.depress.ord.wlsmv.listwise.Rout

The Mplus output for the same example is saved in the file: cfa-depress-wlsmv-listwise.out

  • We specify WLSMV as the estimator. This is WLS with corrections for the standard errors and diagnostic statistics.
  • The option mimic = "Mplus" is used in our effort to synchronize these results with the estimates from Mplus in the companion exercise.
  • The estimates are unchanged if we insert missing = "listwise" in the function call. That’s because listwise deletion is the current default.
  • The summary method includes standardized = TRUE, which does not alter the main parameter estimates, but prints supplementary information in the form of standardized coefficient estimates on the right side of the table.

It is not truly necessary to use ordinal factor variables to estimate this model. We could obtain the same result by declaring the variables explicitly as ordered. Recall that in our formula cfa.depress we used the names of the numeric variables. The lavaan cfa function has an argument ordered that allows us to specify that some variables should be treated as categorical:

cfa.depress.wlsmv2 <-
  cfa(model = cfa.depress, data = hbsc,
    mimic = "Mplus", estimator = "WLSMV",
    ordered = c("depress1", "depress2", "depress3",
                "depress4", "depress5", "depress6"))
  • Please note that the “Used Number of observations” in the WLSMV model is 4,103 compared to 4,221 in the ML FIML model where we did not enforce listwise deletion.

  • Based on the listwise-deleted data (N = 4,103), the WLSMV estimation produced better model fit on CFI, TLI, RMSEA, and SRMR.

A summary table comparing the 2 ML models, along with the WLSMV estimates for the ordinal model, is provided next.

ML (FIML) ML (listwise) WLSMV
Estimate(Std.Err.) Estimate(Std.Err.) Estimate(Std.Err.)
Factor Loadings
depress
depress1 0.73(0.02)*** 0.74(0.02)*** 0.71(0.01)***
depress2 0.71(0.02)*** 0.71(0.02)*** 0.66(0.01)***
depress3 0.79(0.02)*** 0.79(0.02)*** 0.74(0.01)***
depress4 0.86(0.02)*** 0.86(0.02)*** 0.72(0.01)***
depress5 0.75(0.02)*** 0.75(0.02)*** 0.60(0.01)***
depress6 0.73(0.02)*** 0.73(0.02)*** 0.59(0.01)***
Intercepts
depress1 2.33(0.02)*** 2.33(0.02)*** 0.00+
depress2 2.69(0.02)*** 2.68(0.02)*** 0.00+
depress3 1.82(0.02)*** 1.81(0.02)*** 0.00+
depress4 2.21(0.02)*** 2.20(0.02)*** 0.00+
depress5 2.49(0.02)*** 2.48(0.02)*** 0.00+
depress6 2.47(0.02)*** 2.47(0.02)*** 0.00+
Residual Variances
depress1 0.69(0.02)*** 0.67(0.02)*** 0.50+
depress2 0.82(0.02)*** 0.82(0.02)*** 0.57+
depress3 0.84(0.03)*** 0.82(0.03)*** 0.45+
depress4 1.01(0.03)*** 1.01(0.03)*** 0.49+
depress5 1.35(0.03)*** 1.34(0.03)*** 0.64+
depress6 1.27(0.03)*** 1.26(0.03)*** 0.65+
Latent Intercepts
depress 0.00+ 0.00+ 0.00+
Latent Variances
depress 1.00+ 1.00+ 1.00+
Thresholds
depress1(1) -0.56(0.02)***
depress1(2) 0.13(0.02)***
depress1(3) 1.10(0.02)***
depress1(4) 1.79(0.04)***
depress2(1) -0.87(0.02)***
depress2(2) -0.19(0.02)***
depress2(3) 0.74(0.02)***
depress2(4) 1.51(0.03)***
depress3(1) 0.29(0.02)***
depress3(2) 0.67(0.02)***
depress3(3) 1.18(0.03)***
depress3(4) 1.62(0.03)***
depress4(1) -0.11(0.02)***
depress4(2) 0.26(0.02)***
depress4(3) 0.89(0.02)***
depress4(4) 1.44(0.03)***
depress5(1) -0.37(0.02)***
depress5(2) 0.07(0.02)***
depress5(3) 0.67(0.02)***
depress5(4) 1.23(0.03)***
depress6(1) -0.43(0.02)***
depress6(2) 0.08(0.02)***
depress6(3) 0.74(0.02)***
depress6(4) 1.25(0.03)***
Fit Indices
Scaled χ2 198.93(9)*** 195.85(9)*** 315.34(9)***
CFI 0.96 0.96 0.99
TLI 0.93 0.93 0.98
RMSEA 0.08 0.08 0.07
+Fixed parameter
p<0.05, p<0.01, p<0.001


The ordinal model is, of course, “longer”, in the sense that it estimates ordinal thresholds between categories. The coefficients are thus not directly comparable. However, it is interesting to compare the them and wonder if treating the variables as numeric or ordinal really makes “such a big difference.”

The difference between ML with and without listwise deletion is a very interesting topic that we have worked on in other reports. In these examples, we are intereted in comparing the WLSMV estimates against the ML estimates based on the exact same (listwise deleted) data. After this point, whenever we refer to ML estimates, we refer to ML with listwise deletion.

Mplus offers a FIML option for categorical data (not available in lavaan yet). Please see cfa-depress-mlr-cat.html

One-Factor CFA for Bully Victimization

The next construct we present is bullying victimization, which was measured by nine Likert-type items. Items were framed with the following item stem:

“Here are some questions about bullying. We say a student is BEING BULLIED when another student, or a group of students, say or do nasty or unpleasant things to him or her. It is also bullying when a student is teased repeatedly in a way he or she does not like or when they are deliberately left out of things. But it is NOT BULLYING when two students of about the same strength or power argue or fight. It is also not bullying when a student is teased in a friendly and playful way.

How often have you been bullied at school in the past couple of months in the ways listed below?"

Example items include “I was called mean names, was made fun of, or teased in a hurtful way” and “I was hit, kicked, pushed, shoved around, or locked indoors.”

Item response options include: 1) I haven’t been bullied in this way, 2) Only once or twice, 3) 2 or 3 times a month, 4) About once a week, and 5) Several times a week.

The first model is estimated with ML treating items as integers in R. This is a common yet inappropriate practice with ordinal, Likert response types.

Method 1 - Treat ordinal data as continuous with ML

Items Coded as Integers

cfa.bullied <- '
bullied =~ NA*bullied1 + bullied2 + bullied3
            + bullied4 + bullied5 + bullied6 
            + bullied7 + bullied8 + bullied9
bullied ~~ 1*bullied
'
cfa.bullied.mlr.listwise <-
    cfa(model = cfa.bullied, data = hbsc, mimic = "Mplus",
        estimator = "MLR", missing = "listwise", meanstructure = TRUE)
summary(cfa.bullied.mlr.listwise, fit.measures = TRUE, standardized = TRUE)

The lavaan output is saved in the file: output/cfa.bullied.mlr.listwise.Rout

The Mplus output for the same example is saved in the file: cfa-bullied-mlr-listwise.out

Method 2 - The ordinal “threshold” model with WLS

Items Coded as Ordered-Factors

cfa.bullied.ord <- '
bullied =~ NA*bullied1_o + bullied2_o + bullied3_o
            + bullied4_o + bullied5_o + bullied6_o 
            + bullied7_o + bullied8_o + bullied9_o
bullied ~~ 1*bullied
'
cfa.bullied.ord.wlsmv <-
    cfa(model = cfa.bullied.ord, data = hbsc,  mimic = "Mplus",
        estimator = "WLSMV", meanstructure = TRUE)

The model summary table can be printed by running:

summary(cfa.bullied.ord.wlsmv, fit.measures = TRUE, standardized = TRUE)

The lavaan output is saved in output/cfa.bullied.ord.wlsmv.listwise.Rout.

The Mplus output for the same example is saved in the file: cfa-bullied-wlsmv-listwise.out

Rather than discuss the raw output, we prefer to align the output as a more pleasant table, which can be done with the semTable function in kutils.

ML (listwise) WLSMV
Estimate(Std.Err.) Estimate(Std.Err.)
Factor Loadings
bullied
bullied1 0.81(0.03)*** 0.76(0.01)***
bullied2 0.75(0.03)*** 0.78(0.01)***
bullied3 0.65(0.03)*** 0.79(0.01)***
bullied4 0.83(0.03)*** 0.80(0.01)***
bullied5 0.65(0.03)*** 0.83(0.01)***
bullied6 0.54(0.03)*** 0.87(0.01)***
bullied7 0.74(0.03)*** 0.76(0.01)***
bullied8 0.47(0.04)*** 0.87(0.01)***
bullied9 0.44(0.04)*** 0.89(0.01)***
Intercepts
bullied1 1.83(0.02)*** 0.00+
bullied2 1.64(0.02)*** 0.00+
bullied3 1.38(0.02)*** 0.00+
bullied4 1.75(0.02)*** 0.00+
bullied5 1.34(0.02)*** 0.00+
bullied6 1.22(0.01)*** 0.00+
bullied7 1.56(0.02)*** 0.00+
bullied8 1.19(0.01)*** 0.00+
bullied9 1.16(0.01)*** 0.00+
Residual Variances
bullied1 1.01(0.05)*** 0.43+
bullied2 0.77(0.04)*** 0.39+
bullied3 0.48(0.03)*** 0.37+
bullied4 0.77(0.04)*** 0.36+
bullied5 0.41(0.03)*** 0.30+
bullied6 0.27(0.02)*** 0.24+
bullied7 0.73(0.04)*** 0.42+
bullied8 0.27(0.02)*** 0.25+
bullied9 0.26(0.02)*** 0.20+
Latent Intercepts
bullied 0.00+ 0.00+
Latent Variances
bullied 1.00+ 1.00+
Thresholds
bullied1(1) 0.26(0.02)***
bullied1(2) 0.85(0.03)***
bullied1(3) 1.06(0.03)***
bullied1(4) 1.32(0.03)***
bullied2(1) 0.47(0.02)***
bullied2(2) 1.03(0.03)***
bullied2(3) 1.24(0.03)***
bullied2(4) 1.54(0.04)***
bullied3(1) 0.90(0.03)***
bullied3(2) 1.31(0.03)***
bullied3(3) 1.53(0.04)***
bullied3(4) 1.78(0.04)***
bullied4(1) 0.31(0.02)***
bullied4(2) 0.91(0.03)***
bullied4(3) 1.19(0.03)***
bullied4(4) 1.46(0.04)***
bullied5(1) 1.00(0.03)***
bullied5(2) 1.37(0.03)***
bullied5(3) 1.56(0.04)***
bullied5(4) 1.83(0.05)***
bullied6(1) 1.25(0.03)***
bullied6(2) 1.57(0.04)***
bullied6(3) 1.79(0.04)***
bullied6(4) 2.03(0.05)***
bullied7(1) 0.65(0.03)***
bullied7(2) 1.08(0.03)***
bullied7(3) 1.29(0.03)***
bullied7(4) 1.57(0.04)***
bullied8(1) 1.35(0.03)***
bullied8(2) 1.68(0.04)***
bullied8(3) 1.86(0.05)***
bullied8(4) 2.06(0.05)***
bullied9(1) 1.47(0.04)***
bullied9(2) 1.73(0.04)***
bullied9(3) 1.87(0.05)***
bullied9(4) 2.08(0.06)***
Fit Indices
Scaled χ2 536.90(27)*** 551.13(27)***
CFI 0.85 0.99
TLI 0.80 0.99
RMSEA 0.15 0.06
+Fixed parameter
p<0.05, p<0.01, p<0.001


Mplus offers a FIML option for categorical data (not available in lavaan yet). Please see cfa-bullied-mlr-cat.html

One-Factor CFA for Alcohol Use

The final construct we present of our single-factor CFA models is Alcohol Use. This construct was measured by five Likert-type items. Example items are as follows:

At present, how often do you drink anything alcoholic, such as beer,
wine, or hard liquor like Vodka or rum?

A. Beer

Try to include even those times when you only drink a small amount
(e.g., one or two sips)."

Items differed in the type of alcohol, including Beer, Wine, Liquor/Spirits, Pre-mixed drinks (for example, Smirnoff Ice, Bacardi Breezer, Mike’s Hard Lemonade), and Any other drink that contains alcohol.

Response options indicated the frequency of alcohol use, such as 1) Everyday, 2) Every week, 3) Every month, 4) Rarely, 5) Never.

Method 1 - Treat ordinal data as continuous with ML

Items Coded as Integers

cfa.alc <- '
alcohol =~ NA*alc1 + alc2 + alc3
            + alc4 + alc5
alcohol ~~ 1*alcohol
'
cfa.alc.mlr.listwise <-
    cfa(model = cfa.alc, data = hbsc, mimic = "Mplus",
        estimator = "MLR", missing = "listwise", meanstructure = TRUE)
summary(cfa.alc.mlr.listwise, fit.measures = TRUE)

The lavaan output is saved in the file: output/cfa.alc.mlr.listwise.Rout

The Mplus output for the same example is saved in the file: cfa-alcohol-mlr-listwise.out

Method 2 - The ordinal “threshold” model with WLS

Items Coded as Ordered-Factors

cfa.alc.ord <- '
alcohol =~ NA*alc1_o + alc2_o + alc3_o
            + alc4_o + alc5_o
alcohol ~~ 1*alcohol
'

cfa.alc.ord.wlsmv <-
    cfa(model = cfa.alc.ord, data = hbsc, mimic = "Mplus",
        estimator = "WLSMV", missing = "listwise")
summary(cfa.alc.ord.wlsmv, fit.measures = TRUE, standardized = TRUE)

The lavaan summary output is saved in output/cfa.alc.ord.wlsmv.listwise.Rout.

The Mplus output for the same example is saved in the file: cfa-alcohol-wlsmv-listwise.out

Comparing the MLR and WLSMV estimates side by side, we offer the following table:

ML (listwise) WLSMV
Estimate(Std.Err.) Estimate(Std.Err.)
Factor Loadings
alcohol
alc1 0.45(0.03)*** 0.87(0.01)***
alc2 0.41(0.03)*** 0.79(0.01)***
alc3 0.48(0.03)*** 0.95(0.01)***
alc4 0.53(0.02)*** 0.88(0.01)***
alc5 0.52(0.02)*** 0.92(0.01)***
Intercepts
alc1 1.19(0.01)*** 0.00+
alc2 1.22(0.01)*** 0.00+
alc3 1.14(0.01)*** 0.00+
alc4 1.24(0.01)*** 0.00+
alc5 1.20(0.01)*** 0.00+
Residual Variances
alc1 0.13(0.01)*** 0.24+
alc2 0.20(0.01)*** 0.38+
alc3 0.09(0.01)*** 0.10+
alc4 0.18(0.02)*** 0.22+
alc5 0.11(0.01)*** 0.16+
Latent Intercepts
alcohol 0.00+ 0.00+
Latent Variances
alcohol 1.00+ 1.00+
Thresholds
alc1(1) 1.13(0.03)***
alc1(2) 1.86(0.04)***
alc1(3) 2.17(0.05)***
alc1(4) 2.32(0.06)***
alc2(1) 0.99(0.02)***
alc2(2) 1.84(0.04)***
alc2(3) 2.09(0.05)***
alc2(4) 2.41(0.06)***
alc3(1) 1.37(0.03)***
alc3(2) 1.85(0.04)***
alc3(3) 2.09(0.05)***
alc3(4) 2.35(0.06)***
alc4(1) 1.03(0.02)***
alc4(2) 1.62(0.03)***
alc4(3) 1.96(0.04)***
alc4(4) 2.21(0.05)***
alc5(1) 1.13(0.03)***
alc5(2) 1.76(0.04)***
alc5(3) 2.04(0.05)***
alc5(4) 2.31(0.06)***
Fit Indices
Scaled χ2 12.95(5)* 31.14(5)***
CFI 0.99 1.00
TLI 0.99 1.00
RMSEA 0.06 0.02
+Fixed parameter
p<0.05, p<0.01, p<0.001


Mplus offers a FIML option for categorical data (not available in lavaan yet). Please see cfa-alcohol-mlr-cat.html

SECTION-II: A Three-Factor CFA

Factors: Bullying Victimization, Depression, and Alcohol Use

Method 1 - Treat ordinal data as continuous with ML

cfa.3factor <- '
bullied =~ NA*bullied1 + bullied2 + bullied3
            + bullied4 + bullied5 + bullied6
            + bullied7 + bullied8 + bullied9
bullied ~~ 1*bullied
depress =~ NA*depress1 + depress2 + depress3
            + depress4 + depress5 + depress6
depress ~~ 1*depress
alcohol =~ NA*alc1 + alc2 + alc3
            + alc4 + alc5
alcohol ~~ 1*alcohol
'
cfa.3factor.mlr.listwise <-
    cfa(model = cfa.3factor, data = hbsc,  mimic = "Mplus",
        estimator = "MLR", missing = "listwise", meanstructure = TRUE)
summary(cfa.3factor.mlr.listwise, fit.measures = TRUE, standardized = TRUE)
lavaan 0.6-3 ended normally after 35 iterations

  Optimization method                           NLMINB
  Number of free parameters                         63

                                                  Used       Total
  Number of observations                          2684        4284

  Estimator                                         ML      Robust
  Model Fit Test Statistic                    2201.447    1225.318
  Degrees of freedom                               167         167
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  1.797
    for the Yuan-Bentler correction (Mplus variant)

Model test baseline model:

  Minimum Function Test Statistic            22007.778   11122.450
  Degrees of freedom                               190         190
  P-value                                        0.000       0.000

User model versus baseline model:

  Comparative Fit Index (CFI)                    0.907       0.903
  Tucker-Lewis Index (TLI)                       0.894       0.890

  Robust Comparative Fit Index (CFI)                         0.912
  Robust Tucker-Lewis Index (TLI)                            0.900

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)             -61512.306  -61512.306
  Scaling correction factor                                  3.082
    for the MLR correction
  Loglikelihood unrestricted model (H1)     -60411.583  -60411.583
  Scaling correction factor                                  2.149
    for the MLR correction

  Number of free parameters                         63          63
  Akaike (AIC)                              123150.612  123150.612
  Bayesian (BIC)                            123522.001  123522.001
  Sample-size adjusted Bayesian (BIC)       123321.831  123321.831

Root Mean Square Error of Approximation:

  RMSEA                                          0.067       0.049
  90 Percent Confidence Interval          0.065  0.070       0.047  0.051
  P-value RMSEA <= 0.05                          0.000       0.886

  Robust RMSEA                                               0.065
  90 Percent Confidence Interval                             0.062  0.069

Standardized Root Mean Square Residual:

  SRMR                                           0.046       0.046

Parameter Estimates:

  Information                                 Observed
  Observed information based on                Hessian
  Standard Errors                   Robust.huber.white

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  bullied =~                                                            
    bullied1          0.808    0.032   25.356    0.000    0.808    0.626
    bullied2          0.748    0.030   24.948    0.000    0.748    0.648
    bullied3          0.642    0.031   20.851    0.000    0.642    0.686
    bullied4          0.823    0.028   29.184    0.000    0.823    0.686
    bullied5          0.655    0.032   20.204    0.000    0.655    0.718
    bullied6          0.541    0.035   15.454    0.000    0.541    0.723
    bullied7          0.734    0.030   24.851    0.000    0.734    0.650
    bullied8          0.459    0.037   12.516    0.000    0.459    0.671
    bullied9          0.436    0.038   11.450    0.000    0.436    0.653
  depress =~                                                            
    depress1          0.732    0.022   33.626    0.000    0.732    0.665
    depress2          0.691    0.022   31.109    0.000    0.691    0.609
    depress3          0.784    0.026   29.635    0.000    0.784    0.654
    depress4          0.850    0.025   34.495    0.000    0.850    0.653
    depress5          0.754    0.027   27.577    0.000    0.754    0.551
    depress6          0.751    0.027   28.068    0.000    0.751    0.563
  alcohol =~                                                            
    alc1              0.431    0.029   14.787    0.000    0.431    0.753
    alc2              0.383    0.030   12.961    0.000    0.383    0.640
    alc3              0.483    0.029   16.721    0.000    0.483    0.852
    alc4              0.522    0.029   18.228    0.000    0.522    0.784
    alc5              0.522    0.030   17.684    0.000    0.522    0.831

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  bullied ~~                                                            
    depress           0.371    0.031   12.101    0.000    0.371    0.371
    alcohol           0.236    0.041    5.809    0.000    0.236    0.236
  depress ~~                                                            
    alcohol           0.258    0.029    8.778    0.000    0.258    0.258

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .bullied1          1.828    0.025   73.343    0.000    1.828    1.416
   .bullied2          1.640    0.022   73.567    0.000    1.640    1.420
   .bullied3          1.371    0.018   75.906    0.000    1.371    1.465
   .bullied4          1.747    0.023   75.420    0.000    1.747    1.456
   .bullied5          1.334    0.018   75.728    0.000    1.334    1.462
   .bullied6          1.221    0.014   84.433    0.000    1.221    1.630
   .bullied7          1.556    0.022   71.399    0.000    1.556    1.378
   .bullied8          1.180    0.013   89.430    0.000    1.180    1.726
   .bullied9          1.157    0.013   89.822    0.000    1.157    1.734
   .depress1          2.317    0.021  109.091    0.000    2.317    2.106
   .depress2          2.676    0.022  122.156    0.000    2.676    2.358
   .depress3          1.803    0.023   77.877    0.000    1.803    1.503
   .depress4          2.217    0.025   88.162    0.000    2.217    1.702
   .depress5          2.502    0.026   94.634    0.000    2.502    1.827
   .depress6          2.474    0.026   96.149    0.000    2.474    1.856
   .alc1              1.199    0.011  108.614    0.000    1.199    2.096
   .alc2              1.235    0.012  106.958    0.000    1.235    2.065
   .alc3              1.155    0.011  105.616    0.000    1.155    2.039
   .alc4              1.249    0.013   97.270    0.000    1.249    1.878
   .alc5              1.216    0.012  100.316    0.000    1.216    1.936
    bullied           0.000                               0.000    0.000
    depress           0.000                               0.000    0.000
    alcohol           0.000                               0.000    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    bullied           1.000                               1.000    1.000
    depress           1.000                               1.000    1.000
    alcohol           1.000                               1.000    1.000
   .bullied1          1.014    0.050   20.078    0.000    1.014    0.608
   .bullied2          0.775    0.043   18.218    0.000    0.775    0.581
   .bullied3          0.463    0.029   15.803    0.000    0.463    0.529
   .bullied4          0.763    0.041   18.789    0.000    0.763    0.530
   .bullied5          0.403    0.029   14.034    0.000    0.403    0.485
   .bullied6          0.268    0.022   11.981    0.000    0.268    0.478
   .bullied7          0.737    0.042   17.508    0.000    0.737    0.578
   .bullied8          0.257    0.021   12.232    0.000    0.257    0.549
   .bullied9          0.255    0.022   11.804    0.000    0.255    0.573
   .depress1          0.675    0.024   27.908    0.000    0.675    0.558
   .depress2          0.811    0.027   29.864    0.000    0.811    0.629
   .depress3          0.823    0.032   25.357    0.000    0.823    0.572
   .depress4          0.975    0.036   27.433    0.000    0.975    0.574
   .depress5          1.307    0.041   31.654    0.000    1.307    0.697
   .depress6          1.214    0.039   30.814    0.000    1.214    0.683
   .alc1              0.142    0.013   11.306    0.000    0.142    0.433
   .alc2              0.211    0.016   12.861    0.000    0.211    0.590
   .alc3              0.088    0.011    7.746    0.000    0.088    0.274
   .alc4              0.170    0.017   10.184    0.000    0.170    0.385
   .alc5              0.122    0.013    9.661    0.000    0.122    0.309

The Mplus output for the same example is saved in the file: cfa-3factor-mlr-listwise.out

Method 2 - The ordinal “threshold” model with WLS

Items Coded as Ordered-Factors

cfa.3factor.ord <- '
bullied =~ NA*bullied1_o + bullied2_o + bullied3_o
            + bullied4_o + bullied5_o + bullied6_o
            + bullied7_o + bullied8_o + bullied9_o
bullied ~~ 1*bullied
depress =~ NA*depress1_o + depress2_o + depress3_o
            + depress4_o + depress5_o + depress6_o
depress ~~ 1*depress
alcohol =~ NA*alc1_o + alc2_o + alc3_o
            + alc4_o + alc5_o
alcohol ~~ 1*alcohol
'
cfa.3factor.wlsmv.listwise <-
  cfa(model = cfa.3factor.ord, data = hbsc,
      mimic = "Mplus", estimator = "WLSMV")
summary(cfa.3factor.wlsmv.listwise, fit.measures = TRUE, standardized = TRUE)

The Mplus output for the same example is saved in the file: cfa-3factor-wlsmv-listwise.out

Mplus offers a FIML option for categorical data (not available in lavaan yet). Please see cfa-3factor-mlr-cat.html

ML (listwise) WLSMV
Estimate(Std.Err.) Estimate(Std.Err.)
Factor Loadings
bullied
bullied1 0.81(0.03)*** 0.75(0.01)***
bullied2 0.75(0.03)*** 0.78(0.01)***
bullied3 0.64(0.03)*** 0.78(0.01)***
bullied4 0.82(0.03)*** 0.81(0.01)***
bullied5 0.66(0.03)*** 0.84(0.01)***
bullied6 0.54(0.04)*** 0.87(0.01)***
bullied7 0.73(0.03)*** 0.77(0.01)***
bullied8 0.46(0.04)*** 0.86(0.01)***
bullied9 0.44(0.04)*** 0.89(0.02)***
depress
depress1 0.73(0.02)*** 0.70(0.03)***
depress2 0.69(0.02)*** 0.69(0.03)***
depress3 0.78(0.03)*** 0.78(0.03)***
depress4 0.85(0.02)*** 0.84(0.04)***
depress5 0.75(0.03)*** 0.73(0.04)***
depress6 0.75(0.03)*** 0.82(0.04)***
alcohol
alc1 0.43(0.03)*** 0.86(0.01)***
alc2 0.38(0.03)*** 0.77(0.02)***
alc3 0.48(0.03)*** 0.95(0.01)***
alc4 0.52(0.03)*** 0.88(0.01)***
alc5 0.52(0.03)*** 0.91(0.01)***
Intercepts
bullied1 1.83(0.02)*** 0.00+
bullied2 1.64(0.02)*** 0.00+
bullied3 1.37(0.02)*** 0.00+
bullied4 1.75(0.02)*** 0.00+
bullied5 1.33(0.02)*** 0.00+
bullied6 1.22(0.01)*** 0.00+
bullied7 1.56(0.02)*** 0.00+
bullied8 1.18(0.01)*** 0.00+
bullied9 1.16(0.01)*** 0.00+
depress1 2.32(0.02)*** 2.32(0.02)***
depress2 2.68(0.02)*** 2.68(0.02)***
depress3 1.80(0.02)*** 1.80(0.04)***
depress4 2.22(0.03)*** 2.22(0.03)***
depress5 2.50(0.03)*** 2.50(0.03)***
depress6 2.47(0.03)*** 2.47(0.03)***
alc1 1.20(0.01)*** 0.00+
alc2 1.24(0.01)*** 0.00+
alc3 1.15(0.01)*** 0.00+
alc4 1.25(0.01)*** 0.00+
alc5 1.22(0.01)*** 0.00+
Residual Variances
bullied1 1.01(0.05)*** 0.44+
bullied2 0.77(0.04)*** 0.39+
bullied3 0.46(0.03)*** 0.39+
bullied4 0.76(0.04)*** 0.35+
bullied5 0.40(0.03)*** 0.29+
bullied6 0.27(0.02)*** 0.25+
bullied7 0.74(0.04)*** 0.41+
bullied8 0.26(0.02)*** 0.26+
bullied9 0.26(0.02)*** 0.20+
depress1 0.68(0.02)*** 0.73(0.03)***
depress2 0.81(0.03)*** 0.82(0.03)***
depress3 0.82(0.03)*** 0.83(0.04)***
depress4 0.97(0.04)*** 1.00(0.04)***
depress5 1.31(0.04)*** 1.34(0.06)***
depress6 1.21(0.04)*** 1.10(0.05)***
alc1 0.14(0.01)*** 0.26+
alc2 0.21(0.02)*** 0.40+
alc3 0.09(0.01)*** 0.09+
alc4 0.17(0.02)*** 0.22+
alc5 0.12(0.01)*** 0.18+
Latent Intercepts
bullied 0.00+ 0.00+
depress 0.00+ 0.00+
alcohol 0.00+ 0.00+
Latent Variances
bullied 1.00+ 1.00+
depress 1.00+ 1.00+
alcohol 1.00+ 1.00+
Latent Covariances
bullied w/depress 0.37(0.03)*** 0.43(0.02)***
bullied w/alcohol 0.24(0.04)*** 0.26(0.03)***
depress w/alcohol 0.26(0.03)*** 0.34(0.03)***
Thresholds
bullied1(1) 0.27(0.02)***
bullied1(2) 0.85(0.03)***
bullied1(3) 1.07(0.03)***
bullied1(4) 1.32(0.03)***
bullied2(1) 0.47(0.03)***
bullied2(2) 1.03(0.03)***
bullied2(3) 1.24(0.03)***
bullied2(4) 1.54(0.04)***
bullied3(1) 0.91(0.03)***
bullied3(2) 1.33(0.03)***
bullied3(3) 1.55(0.04)***
bullied3(4) 1.80(0.05)***
bullied4(1) 0.30(0.02)***
bullied4(2) 0.92(0.03)***
bullied4(3) 1.20(0.03)***
bullied4(4) 1.46(0.04)***
bullied5(1) 1.01(0.03)***
bullied5(2) 1.37(0.03)***
bullied5(3) 1.55(0.04)***
bullied5(4) 1.83(0.05)***
bullied6(1) 1.25(0.03)***
bullied6(2) 1.58(0.04)***
bullied6(3) 1.79(0.05)***
bullied6(4) 2.03(0.05)***
bullied7(1) 0.65(0.03)***
bullied7(2) 1.08(0.03)***
bullied7(3) 1.29(0.03)***
bullied7(4) 1.56(0.04)***
bullied8(1) 1.36(0.03)***
bullied8(2) 1.69(0.04)***
bullied8(3) 1.89(0.05)***
bullied8(4) 2.08(0.06)***
bullied9(1) 1.49(0.04)***
bullied9(2) 1.74(0.04)***
bullied9(3) 1.89(0.05)***
bullied9(4) 2.07(0.06)***
alc1(1) 1.06(0.03)***
alc1(2) 1.84(0.05)***
alc1(3) 2.19(0.06)***
alc1(4) 2.43(0.08)***
alc2(1) 0.93(0.03)***
alc2(2) 1.82(0.05)***
alc2(3) 2.12(0.06)***
alc2(4) 2.47(0.08)***
alc3(1) 1.32(0.03)***
alc3(2) 1.82(0.05)***
alc3(3) 2.07(0.06)***
alc3(4) 2.43(0.08)***
alc4(1) 0.98(0.03)***
alc4(2) 1.60(0.04)***
alc4(3) 2.00(0.05)***
alc4(4) 2.35(0.07)***
alc5(1) 1.07(0.03)***
alc5(2) 1.72(0.04)***
alc5(3) 2.04(0.05)***
alc5(4) 2.35(0.07)***
Fit Indices
Scaled χ2 1225.32(167)*** 993.15(167)***
CFI 0.91 0.99
TLI 0.89 0.99
RMSEA 0.07 0.04
+Fixed parameter
p<0.05, p<0.01, p<0.001


One of the questions in CFA methodology that arises from time-to-time is whether one ought to fit CFA models for each factor separately, or if one should fit all of the separate factors in the first step. Some veterans of SEM suggest one ought to dispense with the individual factor estimates and proceed right to this step, where we have 3 sets of observed indicators and 3 factors.

The following table provides 4 sets of estimates side by side. The estimates for the 3 factor model are in the 4th column, and we are interested to know if the estimates of the factor loadings in the individual models are called into question by the larger model.

Alcohol Bullied Depression 3 factor
Estimate(Std.Err.) Estimate(Std.Err.) Estimate(Std.Err.) Estimate(Std.Err.)
Factor Loadings
alcohol
alc1 0.87(0.01)*** 0.86(0.01)***
alc2 0.79(0.01)*** 0.77(0.02)***
alc3 0.95(0.01)*** 0.95(0.01)***
alc4 0.88(0.01)*** 0.88(0.01)***
alc5 0.92(0.01)*** 0.91(0.01)***
bullied
bullied1 0.76(0.01)*** 0.75(0.01)***
bullied2 0.78(0.01)*** 0.78(0.01)***
bullied3 0.79(0.01)*** 0.78(0.01)***
bullied4 0.80(0.01)*** 0.81(0.01)***
bullied5 0.83(0.01)*** 0.84(0.01)***
bullied6 0.87(0.01)*** 0.87(0.01)***
bullied7 0.76(0.01)*** 0.77(0.01)***
bullied8 0.87(0.01)*** 0.86(0.01)***
bullied9 0.89(0.01)*** 0.89(0.02)***
depress
depress1 0.71(0.01)*** 0.70(0.03)***
depress2 0.66(0.01)*** 0.69(0.03)***
depress3 0.74(0.01)*** 0.78(0.03)***
depress4 0.72(0.01)*** 0.84(0.04)***
depress5 0.60(0.01)*** 0.73(0.04)***
depress6 0.59(0.01)*** 0.82(0.04)***
Intercepts
alc1 0.00+ 0.00+
alc2 0.00+ 0.00+
alc3 0.00+ 0.00+
alc4 0.00+ 0.00+
alc5 0.00+ 0.00+
bullied1 0.00+ 0.00+
bullied2 0.00+ 0.00+
bullied3 0.00+ 0.00+
bullied4 0.00+ 0.00+
bullied5 0.00+ 0.00+
bullied6 0.00+ 0.00+
bullied7 0.00+ 0.00+
bullied8 0.00+ 0.00+
bullied9 0.00+ 0.00+
depress1 0.00+ 2.32(0.02)***
depress2 0.00+ 2.68(0.02)***
depress3 0.00+ 1.80(0.04)***
depress4 0.00+ 2.22(0.03)***
depress5 0.00+ 2.50(0.03)***
depress6 0.00+ 2.47(0.03)***
Residual Variances
alc1 0.24+ 0.26+
alc2 0.38+ 0.40+
alc3 0.10+ 0.09+
alc4 0.22+ 0.22+
alc5 0.16+ 0.18+
bullied1 0.43+ 0.44+
bullied2 0.39+ 0.39+
bullied3 0.37+ 0.39+
bullied4 0.36+ 0.35+
bullied5 0.30+ 0.29+
bullied6 0.24+ 0.25+
bullied7 0.42+ 0.41+
bullied8 0.25+ 0.26+
bullied9 0.20+ 0.20+
depress1 0.50+ 0.73(0.03)***
depress2 0.57+ 0.82(0.03)***
depress3 0.45+ 0.83(0.04)***
depress4 0.49+ 1.00(0.04)***
depress5 0.64+ 1.34(0.06)***
depress6 0.65+ 1.10(0.05)***
Latent Intercepts
alcohol 0.00+ 0.00+
bullied 0.00+ 0.00+
depress 0.00+ 0.00+
Latent Variances
alcohol 1.00+ 1.00+
bullied 1.00+ 1.00+
depress 1.00+ 1.00+
Latent Covariances
bullied w/depress 0.43(0.02)***
bullied w/alcohol 0.26(0.03)***
depress w/alcohol 0.34(0.03)***
Thresholds
alc1(1) 1.13(0.03)*** 1.06(0.03)***
alc1(2) 1.86(0.04)*** 1.84(0.05)***
alc1(3) 2.17(0.05)*** 2.19(0.06)***
alc1(4) 2.32(0.06)*** 2.43(0.08)***
alc2(1) 0.99(0.02)*** 0.93(0.03)***
alc2(2) 1.84(0.04)*** 1.82(0.05)***
alc2(3) 2.09(0.05)*** 2.12(0.06)***
alc2(4) 2.41(0.06)*** 2.47(0.08)***
alc3(1) 1.37(0.03)*** 1.32(0.03)***
alc3(2) 1.85(0.04)*** 1.82(0.05)***
alc3(3) 2.09(0.05)*** 2.07(0.06)***
alc3(4) 2.35(0.06)*** 2.43(0.08)***
alc4(1) 1.03(0.02)*** 0.98(0.03)***
alc4(2) 1.62(0.03)*** 1.60(0.04)***
alc4(3) 1.96(0.04)*** 2.00(0.05)***
alc4(4) 2.21(0.05)*** 2.35(0.07)***
alc5(1) 1.13(0.03)*** 1.07(0.03)***
alc5(2) 1.76(0.04)*** 1.72(0.04)***
alc5(3) 2.04(0.05)*** 2.04(0.05)***
alc5(4) 2.31(0.06)*** 2.35(0.07)***
bullied1(1) 0.26(0.02)*** 0.27(0.02)***
bullied1(2) 0.85(0.03)*** 0.85(0.03)***
bullied1(3) 1.06(0.03)*** 1.07(0.03)***
bullied1(4) 1.32(0.03)*** 1.32(0.03)***
bullied2(1) 0.47(0.02)*** 0.47(0.03)***
bullied2(2) 1.03(0.03)*** 1.03(0.03)***
bullied2(3) 1.24(0.03)*** 1.24(0.03)***
bullied2(4) 1.54(0.04)*** 1.54(0.04)***
bullied3(1) 0.90(0.03)*** 0.91(0.03)***
bullied3(2) 1.31(0.03)*** 1.33(0.03)***
bullied3(3) 1.53(0.04)*** 1.55(0.04)***
bullied3(4) 1.78(0.04)*** 1.80(0.05)***
bullied4(1) 0.31(0.02)*** 0.30(0.02)***
bullied4(2) 0.91(0.03)*** 0.92(0.03)***
bullied4(3) 1.19(0.03)*** 1.20(0.03)***
bullied4(4) 1.46(0.04)*** 1.46(0.04)***
bullied5(1) 1.00(0.03)*** 1.01(0.03)***
bullied5(2) 1.37(0.03)*** 1.37(0.03)***
bullied5(3) 1.56(0.04)*** 1.55(0.04)***
bullied5(4) 1.83(0.05)*** 1.83(0.05)***
bullied6(1) 1.25(0.03)*** 1.25(0.03)***
bullied6(2) 1.57(0.04)*** 1.58(0.04)***
bullied6(3) 1.79(0.04)*** 1.79(0.05)***
bullied6(4) 2.03(0.05)*** 2.03(0.05)***
bullied7(1) 0.65(0.03)*** 0.65(0.03)***
bullied7(2) 1.08(0.03)*** 1.08(0.03)***
bullied7(3) 1.29(0.03)*** 1.29(0.03)***
bullied7(4) 1.57(0.04)*** 1.56(0.04)***
bullied8(1) 1.35(0.03)*** 1.36(0.03)***
bullied8(2) 1.68(0.04)*** 1.69(0.04)***
bullied8(3) 1.86(0.05)*** 1.89(0.05)***
bullied8(4) 2.06(0.05)*** 2.08(0.06)***
bullied9(1) 1.47(0.04)*** 1.49(0.04)***
bullied9(2) 1.73(0.04)*** 1.74(0.04)***
bullied9(3) 1.87(0.05)*** 1.89(0.05)***
bullied9(4) 2.08(0.06)*** 2.07(0.06)***
depress1(1) -0.56(0.02)***
depress1(2) 0.13(0.02)***
depress1(3) 1.10(0.02)***
depress1(4) 1.79(0.04)***
depress2(1) -0.87(0.02)***
depress2(2) -0.19(0.02)***
depress2(3) 0.74(0.02)***
depress2(4) 1.51(0.03)***
depress3(1) 0.29(0.02)***
depress3(2) 0.67(0.02)***
depress3(3) 1.18(0.03)***
depress3(4) 1.62(0.03)***
depress4(1) -0.11(0.02)***
depress4(2) 0.26(0.02)***
depress4(3) 0.89(0.02)***
depress4(4) 1.44(0.03)***
depress5(1) -0.37(0.02)***
depress5(2) 0.07(0.02)***
depress5(3) 0.67(0.02)***
depress5(4) 1.23(0.03)***
depress6(1) -0.43(0.02)***
depress6(2) 0.08(0.02)***
depress6(3) 0.74(0.02)***
depress6(4) 1.25(0.03)***
Fit Indices
Scaled χ2 31.14(5)*** 551.13(27)*** 315.34(9)*** 993.15(167)***
CFI 1.00 0.99 0.99 0.99
TLI 1.00 0.99 0.98 0.99
RMSEA 0.02 0.06 0.07 0.04
+Fixed parameter
p<0.05, p<0.01, p<0.001


SECTION-III: A Two-Factor Structural Model: Bullying Victimization Predicts Alcohol Use

Method 1 - Treat ordinal data as continuous with ML

sem.01 <- '
## the measurement model
bullied =~ NA*bullied1 + bullied2 + bullied3
             + bullied4 + bullied5 + bullied6
             + bullied7 + bullied8 + bullied9
bullied ~~ 1*bullied
alcohol =~ NA*alc1 + alc2 + alc3
            + alc4 + alc5
alcohol ~~ 1*alcohol
# regress 
alcohol ~ bullied
'

sem.01.mlr.listwise <-
    sem(model = sem.01, data = hbsc,  mimic = "Mplus",
        estimator = "MLR", missing = "listwise")
summary(sem.01.mlr.listwise, fit.measures = TRUE,
        standardized = TRUE)

The lavaan summary output is saved in output/sem.01.mlr.listwise.Rout.

The Mplus output for the same example is saved in the file: sem-01-mlr-listwise.out

Method 2 - The ordinal “threshold” model with WLS

Items Coded as Ordered-Factors

sem.01.ord <- '
## the measurement model
bullied =~ NA*bullied1_o + bullied2_o + bullied3_o
             + bullied4_o + bullied5_o + bullied6_o
             + bullied7_o + bullied8_o + bullied9_o
bullied ~~ 1*bullied
alcohol =~ NA*alc1_o + alc2_o + alc3_o
            + alc4_o + alc5_o
alcohol ~~ 1*alcohol
## the structural model
alcohol ~ bullied
'
sem.01.ord.wlsmv.listwise <-
    sem(model = sem.01.ord, data = hbsc, mimic = "Mplus",
        estimator = "WLSMV", missing = "listwise")
summary(sem.01.ord.wlsmv, fit.measures = TRUE, standardized = TRUE)

The lavaan summary output is saved in output/sem.01.ord.wlsmv.listwise.Rout.

The Mplus output for the same example is saved in the file: sem-01-ord-wlsmv-listwise.out

ML (listwise) WLSMV
Estimate(Std.Err.) Estimate(Std.Err.)
Factor Loadings
bullied
bullied1 0.79(0.03)*** 0.74(0.01)***
bullied2 0.73(0.03)*** 0.77(0.01)***
bullied3 0.65(0.03)*** 0.79(0.01)***
bullied4 0.81(0.03)*** 0.80(0.01)***
bullied5 0.66(0.03)*** 0.84(0.01)***
bullied6 0.55(0.03)*** 0.88(0.01)***
bullied7 0.73(0.03)*** 0.76(0.01)***
bullied8 0.48(0.04)*** 0.87(0.01)***
bullied9 0.45(0.04)*** 0.90(0.01)***
alcohol
alc1 0.43(0.03)*** 0.83(0.01)***
alc2 0.38(0.03)*** 0.76(0.02)***
alc3 0.47(0.03)*** 0.91(0.01)***
alc4 0.51(0.03)*** 0.85(0.01)***
alc5 0.51(0.03)*** 0.87(0.01)***
Regression Slopes
alcohol
bullied 0.26(0.05)*** 0.29(0.03)***
Intercepts
bullied1 0.00+
bullied2 0.00+
bullied3 0.00+
bullied4 0.00+
bullied5 0.00+
bullied6 0.00+
bullied7 0.00+
bullied8 0.00+
bullied9 0.00+
alc1 0.00+
alc2 0.00+
alc3 0.00+
alc4 0.00+
alc5 0.00+
Residual Variances
bullied1 1.04(0.05)*** 0.45+
bullied2 0.79(0.04)*** 0.40+
bullied3 0.46(0.03)*** 0.38+
bullied4 0.78(0.04)*** 0.36+
bullied5 0.40(0.03)*** 0.29+
bullied6 0.26(0.02)*** 0.23+
bullied7 0.74(0.04)*** 0.42+
bullied8 0.25(0.02)*** 0.24+
bullied9 0.25(0.02)*** 0.18+
alc1 0.14(0.01)*** 0.25+
alc2 0.21(0.02)*** 0.38+
alc3 0.09(0.01)*** 0.10+
alc4 0.17(0.02)*** 0.22+
alc5 0.12(0.01)*** 0.18+
Latent Intercepts
bullied 0.00+
alcohol 0.00+
Latent Variances
bullied 1.00+ 1.00+
alcohol 1.00+ 1.00+
Thresholds
bullied1(1) 0.27(0.02)***
bullied1(2) 0.85(0.03)***
bullied1(3) 1.07(0.03)***
bullied1(4) 1.32(0.03)***
bullied2(1) 0.47(0.02)***
bullied2(2) 1.03(0.03)***
bullied2(3) 1.24(0.03)***
bullied2(4) 1.54(0.04)***
bullied3(1) 0.90(0.03)***
bullied3(2) 1.32(0.03)***
bullied3(3) 1.54(0.04)***
bullied3(4) 1.80(0.05)***
bullied4(1) 0.30(0.02)***
bullied4(2) 0.91(0.03)***
bullied4(3) 1.20(0.03)***
bullied4(4) 1.46(0.04)***
bullied5(1) 1.01(0.03)***
bullied5(2) 1.37(0.03)***
bullied5(3) 1.55(0.04)***
bullied5(4) 1.83(0.05)***
bullied6(1) 1.25(0.03)***
bullied6(2) 1.57(0.04)***
bullied6(3) 1.78(0.04)***
bullied6(4) 2.02(0.05)***
bullied7(1) 0.65(0.03)***
bullied7(2) 1.08(0.03)***
bullied7(3) 1.29(0.03)***
bullied7(4) 1.57(0.04)***
bullied8(1) 1.36(0.03)***
bullied8(2) 1.68(0.04)***
bullied8(3) 1.87(0.05)***
bullied8(4) 2.06(0.06)***
bullied9(1) 1.48(0.04)***
bullied9(2) 1.74(0.04)***
bullied9(3) 1.88(0.05)***
bullied9(4) 2.08(0.06)***
alc1(1) 1.06(0.03)***
alc1(2) 1.82(0.05)***
alc1(3) 2.17(0.06)***
alc1(4) 2.39(0.08)***
alc2(1) 0.92(0.03)***
alc2(2) 1.81(0.05)***
alc2(3) 2.10(0.06)***
alc2(4) 2.44(0.08)***
alc3(1) 1.31(0.03)***
alc3(2) 1.81(0.05)***
alc3(3) 2.07(0.06)***
alc3(4) 2.43(0.08)***
alc4(1) 0.99(0.03)***
alc4(2) 1.60(0.04)***
alc4(3) 1.99(0.05)***
alc4(4) 2.32(0.07)***
alc5(1) 1.07(0.03)***
alc5(2) 1.72(0.04)***
alc5(3) 2.03(0.05)***
alc5(4) 2.35(0.07)***
Fit Indices
Scaled χ2 727.31(76)*** 646.90(76)***
CFI 0.90 0.99
TLI 0.89 0.99
RMSEA 0.09 0.04
+Fixed parameter
p<0.05, p<0.01, p<0.001


SECTION-IV: Introducing Mediation: Bullying leads to Depression, which also affects Alcohol Use

The direct effect of bullying on alcohol is compared with the indirect effect of bullying on depression, which also contributes to alcohol use.

Method 1 - Treat ordinal data as continuous with ML

Items Coded as Integers

sem.02 <- '
## the measurement model
bullied =~ NA*bullied1 + bullied2 + bullied3
             + bullied4 + bullied5 + bullied6 
             + bullied7 + bullied8 + bullied9
bullied ~~ 1*bullied
depress =~ NA*depress1 + depress2 + depress3
            + depress4 + depress5 + depress6
depress ~~ 1*depress
alcohol =~ NA*alc1 + alc2 + alc3
            + alc4 + alc5
alcohol ~~ 1*alcohol
## the structural model
## direct effect (the c path)
alcohol ~ c*bullied
## mediator paths (the a and b path)
depress ~ a*bullied # the a path - IV predicting ME
alcohol ~ b*depress  # the b path - ME predicting DV
## indirect effect (a*b)
## := operator defines new parameters
ab := a*b
## total effect
total := c + (a*b)
'
NBOOT <- 10
sem.02.mlboot.listwise <- sem(model = sem.02, data = hbsc,
                           mimic = "Mplus", estimator = "ML",
                           missing = "listwise", se = "bootstrap",
                           verbose = TRUE, bootstrap = NBOOT)
Quasi-Newton steps using NLMINB:
Objective function  = 0.7610958025234247
Objective function  =                Inf
Objective function  = 0.6261883155717491
Objective function  = 0.5969380824175232
Objective function  = 0.5574657831870571
Objective function  = 0.5237220884801737
Objective function  = 0.5132137359131050
Objective function  = 0.4614986716887337
Objective function  = 0.4450837143337711
Objective function  = 0.4300752430023813
Objective function  = 0.4258803144428498
Objective function  = 0.4201880581488791
Objective function  = 0.4207021470190515
Objective function  = 0.4166934490861394
Objective function  = 0.4167723630244318
Objective function  = 0.4138707843922820
Objective function  = 0.4159779242726813
Objective function  = 0.4133172849470856
Objective function  = 0.4129658435237467
Objective function  = 0.4123862404077556
Objective function  = 0.4113982513162817
Objective function  = 0.4113902754435905
Objective function  = 0.4103205563367389
Objective function  = 0.4104473408485134
Objective function  = 0.4102194724071797
Objective function  = 0.4101694948808436
Objective function  = 0.4101331615003687
Objective function  = 0.4101299250257604
Objective function  = 0.4101136960227336
Objective function  = 0.4101166471343287
Objective function  = 0.4101085990565920
Objective function  = 0.4101071829062093
Objective function  = 0.4101063575019559
Objective function  = 0.4101058652012473
Objective function  = 0.4101058726397326
Objective function  = 0.4101056859304251
Objective function  = 0.4101056342819795
Objective function  = 0.4101056856153278
Objective function  = 0.4101056080853311
Objective function  = 0.4101055944153860
Objective function  = 0.4101055883698059
Objective function  = 0.4101055850551489
Objective function  = 0.4101055845019559
Objective function  = 0.4101055842358470
Objective function  = 0.4101055842049437
Objective function  = 0.4101055841712000
Objective function  = 0.4101055841712000
convergence status (0=ok):  0 
nlminb message says:  relative convergence (4) 
number of iterations:  36 
number of function evaluations [objective, gradient]:  46 37 
Computing VCOV for se = bootstrap ...
  ... bootstrap draw number:    1   OK -- niter =   37  fx =    0.475041704 
  ... bootstrap draw number:    2   OK -- niter =   35  fx =    0.493201822 
  ... bootstrap draw number:    3   OK -- niter =   34  fx =    0.450459562 
  ... bootstrap draw number:    4   OK -- niter =   32  fx =    0.439189091 
  ... bootstrap draw number:    5   OK -- niter =   37  fx =    0.457102997 
  ... bootstrap draw number:    6   OK -- niter =   33  fx =    0.462109375 
  ... bootstrap draw number:    7   OK -- niter =   33  fx =    0.423492563 
  ... bootstrap draw number:    8   OK -- niter =   34  fx =    0.464663698 
  ... bootstrap draw number:    9   OK -- niter =   33  fx =    0.469361398 
  ... bootstrap draw number:   10   OK -- niter =   34  fx =    0.511879195 
Number of successful bootstrap draws: 10 
Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
    The variance-covariance matrix of the estimated parameters (vcov)
    does not appear to be positive definite! The smallest eigenvalue
    (= -2.792746e-18) is smaller than zero. This may be a symptom that
    the model is not identified.
 done.
Computing TEST for test = standard ... done.
summary(sem.02.mlboot.listwise, fit.measures = TRUE, standardized = TRUE)

The lavaan summary output is saved in output/sem.02.mlboot.listwise.Rout.

The Mplus output for the same example is saved in the file: sem-02-mlboot-listwise.out

Method 2 - The ordinal “threshold” model with WLS

Items Coded as Ordered-Factors

sem.02.ord <- '
## the measurement model
bullied =~ NA*bullied1_o + bullied2_o + bullied3_o
             + bullied4_o + bullied5_o + bullied6_o 
             + bullied7_o + bullied8_o + bullied9_o
bullied ~~ 1*bullied
depress =~ NA*depress1_o + depress2_o + depress3_o
            + depress4_o + depress5_o + depress6_o
depress ~~ 1*depress
alcohol =~ NA*alc1_o + alc2_o + alc3_o
            + alc4_o + alc5_o
alcohol ~~ 1*alcohol
## the structural model
## direct effect (the c path)
alcohol ~ c*bullied
## mediator paths (the a and b path)
depress ~ a*bullied # the a path - IV predicting ME
alcohol ~ b*depress  # the b path - ME predicting DV
## indirect effect (a*b)
## := operator defines new parameters
ab := a*b
## total effect
total := c + (a*b)
'
sem.02.ord.wlsmv.listwise <-
  sem(model = sem.02.ord, data = hbsc, mimic = "Mplus",
      estimator = "DWLS", missing = "listwise",
      se = "bootstrap", verbose = TRUE,
      bootstrap = NBOOT)
Warning in lav_options_set(opt): lavaan WARNING: information will be set to
"expected" for estimator = "DWLS"
Estimating sample thresholds and correlations ... done
Quasi-Newton steps using NLMINB:
Objective function  = 1.6304782009635421
Objective function  = 4.3078918156426198
Objective function  = 0.5235747585056084
Objective function  = 0.4648468857816807
Objective function  = 0.2962602989148983
Objective function  = 0.5158296626894420
Objective function  = 0.1797810390665839
Objective function  = 0.1794633005574890
Objective function  = 0.1639437760930223
Objective function  = 0.1646792256351012
Objective function  = 0.1643922647433192
Objective function  = 0.1631750766018704
Objective function  = 0.1627453483709781
Objective function  = 0.1623158721605166
Objective function  = 0.1623687520419041
Objective function  = 0.1623792467608653
Objective function  = 0.1622108420962807
Objective function  = 0.1621514247201504
Objective function  = 0.1621324780519058
Objective function  = 0.1621388401834119
Objective function  = 0.1621310484917360
Objective function  = 0.1621297447100511
Objective function  = 0.1621298352933153
Objective function  = 0.1621290450641717
Objective function  = 0.1621291173556044
Objective function  = 0.1621289404318876
Objective function  = 0.1621289056282594
Objective function  = 0.1621289055209842
Objective function  = 0.1621288871512102
Objective function  = 0.1621288847412473
Objective function  = 0.1621288850877062
Objective function  = 0.1621288846146514
Objective function  = 0.1621288846087745
Objective function  = 0.1621288845890415
Objective function  = 0.1621288845956790
Objective function  = 0.1621288845890415
convergence status (0=ok):  0 
nlminb message says:  relative convergence (4) 
number of iterations:  24 
number of function evaluations [objective, gradient]:  35 24 
Computing VCOV for se = bootstrap ...
  ... bootstrap draw number:    1   OK -- niter =   33  fx =    0.174321098 
  ... bootstrap draw number:    2   OK -- niter =   29  fx =    0.215058527 
  ... bootstrap draw number:    3   OK -- niter =   31  fx =    0.182289438 
  ... bootstrap draw number:    4   OK -- niter =   30  fx =    0.188415868 
  ... bootstrap draw number:    5   OK -- niter =   32  fx =    0.160597206 
  ... bootstrap draw number:    6   OK -- niter =   25  fx =    0.186750075 
  ... bootstrap draw number:    7   OK -- niter =   26  fx =    0.187428097 
  ... bootstrap draw number:    8   OK -- niter =   27  fx =    0.190863004 
  ... bootstrap draw number:    9   OK -- niter =   27  fx =    0.204454806 
  ... bootstrap draw number:   10   OK -- niter =   30  fx =    0.202095154 
Number of successful bootstrap draws: 10 
Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
    The variance-covariance matrix of the estimated parameters (vcov)
    does not appear to be positive definite! The smallest eigenvalue
    (= -1.200669e-17) is smaller than zero. This may be a symptom that
    the model is not identified.
 done.
Computing TEST for test = standard ... done.
summary(sem.02.ord.wlsmv.listwise, fit.measures = TRUE, standardized = TRUE)

The lavaan summary output is saved in output/sem.02.ord.wlsmv.listwise.Rout.

The Mplus output for the same example is saved in the file: sem-02-ord-wlsmv-listwise.out

Following is the summary table comparing the two models.

ML (listwise) DWLS
Estimate(Std.Err.) Estimate(Std.Err.)
Factor Loadings
bullied
bullied1 0.81(0.04)*** 0.75(0.01)***
bullied2 0.75(0.03)*** 0.78(0.01)***
bullied3 0.64(0.01)*** 0.78(0.01)***
bullied4 0.82(0.02)*** 0.81(0.01)***
bullied5 0.66(0.03)*** 0.84(0.02)***
bullied6 0.54(0.04)*** 0.87(0.02)***
bullied7 0.73(0.03)*** 0.77(0.02)***
bullied8 0.46(0.05)*** 0.86(0.02)***
bullied9 0.44(0.04)*** 0.89(0.02)***
depress
depress1 0.68(0.02)*** 0.63(0.03)***
depress2 0.64(0.02)*** 0.62(0.02)***
depress3 0.73(0.02)*** 0.70(0.02)***
depress4 0.79(0.02)*** 0.76(0.03)***
depress5 0.70(0.01)*** 0.66(0.03)***
depress6 0.70(0.02)*** 0.74(0.03)***
alcohol
alc1 0.41(0.03)*** 0.80(0.02)***
alc2 0.37(0.02)*** 0.72(0.02)***
alc3 0.46(0.02)*** 0.88(0.01)***
alc4 0.50(0.02)*** 0.82(0.02)***
alc5 0.50(0.03)*** 0.84(0.02)***
Regression Slopes
alcohol
bullied 0.17(0.05)*** 0.15(0.05)**
depress 0.19(0.02)*** 0.28(0.04)***
depress
bullied 0.40(0.03)*** 0.47(0.02)***
Intercepts
bullied1 0.00+
bullied2 0.00+
bullied3 0.00+
bullied4 0.00+
bullied5 0.00+
bullied6 0.00+
bullied7 0.00+
bullied8 0.00+
bullied9 0.00+
depress1 2.32(0.01)***
depress2 2.68(0.01)***
depress3 1.80(0.02)***
depress4 2.22(0.02)***
depress5 2.50(0.03)***
depress6 2.47(0.03)***
alc1 0.00+
alc2 0.00+
alc3 0.00+
alc4 0.00+
alc5 0.00+
Residual Variances
bullied1 1.01(0.05)*** 0.44+
bullied2 0.77(0.04)*** 0.39+
bullied3 0.46(0.02)*** 0.39+
bullied4 0.76(0.04)*** 0.35+
bullied5 0.40(0.02)*** 0.29+
bullied6 0.27(0.03)*** 0.25+
bullied7 0.74(0.04)*** 0.41+
bullied8 0.26(0.02)*** 0.26+
bullied9 0.26(0.03)*** 0.20+
depress1 0.68(0.02)*** 0.73(0.04)***
depress2 0.81(0.04)*** 0.82(0.03)***
depress3 0.82(0.02)*** 0.83(0.05)***
depress4 0.97(0.05)*** 1.00(0.04)***
depress5 1.31(0.04)*** 1.34(0.04)***
depress6 1.21(0.03)*** 1.10(0.06)***
alc1 0.14(0.01)*** 0.26+
alc2 0.21(0.01)*** 0.40+
alc3 0.09(0.01)*** 0.09+
alc4 0.17(0.01)*** 0.22+
alc5 0.12(0.02)*** 0.18+
Latent Intercepts
bullied 0.00+
depress 0.00+
alcohol 0.00+
Latent Variances
bullied 1.00+ 1.00+
depress 1.00+ 1.00+
alcohol 1.00+ 1.00+
Thresholds
bullied1(1) 0.27(0.02)***
bullied1(2) 0.85(0.02)***
bullied1(3) 1.07(0.02)***
bullied1(4) 1.32(0.03)***
bullied2(1) 0.47(0.02)***
bullied2(2) 1.03(0.02)***
bullied2(3) 1.24(0.02)***
bullied2(4) 1.54(0.03)***
bullied3(1) 0.91(0.02)***
bullied3(2) 1.33(0.02)***
bullied3(3) 1.55(0.01)***
bullied3(4) 1.80(0.03)***
bullied4(1) 0.30(0.02)***
bullied4(2) 0.92(0.02)***
bullied4(3) 1.20(0.03)***
bullied4(4) 1.46(0.03)***
bullied5(1) 1.01(0.02)***
bullied5(2) 1.37(0.03)***
bullied5(3) 1.55(0.04)***
bullied5(4) 1.83(0.04)***
bullied6(1) 1.25(0.03)***
bullied6(2) 1.58(0.03)***
bullied6(3) 1.79(0.03)***
bullied6(4) 2.03(0.05)***
bullied7(1) 0.65(0.02)***
bullied7(2) 1.08(0.03)***
bullied7(3) 1.29(0.02)***
bullied7(4) 1.56(0.03)***
bullied8(1) 1.36(0.03)***
bullied8(2) 1.69(0.04)***
bullied8(3) 1.89(0.05)***
bullied8(4) 2.08(0.06)***
bullied9(1) 1.49(0.04)***
bullied9(2) 1.74(0.05)***
bullied9(3) 1.89(0.05)***
bullied9(4) 2.07(0.06)***
alc1(1) 1.06(0.03)***
alc1(2) 1.84(0.05)***
alc1(3) 2.19(0.08)***
alc1(4) 2.43(0.06)***
alc2(1) 0.93(0.03)***
alc2(2) 1.82(0.03)***
alc2(3) 2.12(0.05)***
alc2(4) 2.47(0.09)***
alc3(1) 1.32(0.03)***
alc3(2) 1.82(0.04)***
alc3(3) 2.07(0.06)***
alc3(4) 2.43(0.09)***
alc4(1) 0.98(0.02)***
alc4(2) 1.60(0.03)***
alc4(3) 2.00(0.06)***
alc4(4) 2.35(0.07)***
alc5(1) 1.07(0.03)***
alc5(2) 1.72(0.03)***
alc5(3) 2.04(0.06)***
alc5(4) 2.35(0.05)***
Constructed
ab 0.08(0.01)*** 0.13(0.02)***
total 0.25(0.05)*** 0.28(0.06)***
Fit Indices
χ2 2201.45(167)*** 745.10(167)***
CFI 0.91 0.99
TLI 0.89 0.99
RMSEA 0.07 0.04
+Fixed parameter
p<0.05, p<0.01, p<0.001


Summary

When the coefficients of the 2 kinds of models are rescaled onto a more-or-less common scale, the loadings and estimates of underlying structures are very similar. In this example, we don’t find coefficients that change from “significant” to “insignificant” when comparing the ML-as-if-numeric and WLSMV estimates.

If we allow the ML estimator to use FIML to avoid dropping cases, the estimates are interestingly different.

R Session Info

R version 3.5.1 (2018-07-02)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] xtable_1.8-3        lavaan_0.6-3        kutils_1.62        
[4] stationery_0.98.5.7

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0      digest_0.6.18   MASS_7.3-51.1   plyr_1.8.4     
 [5] magrittr_1.5    stats4_3.5.1    evaluate_0.12   zip_1.0.0      
 [9] stringi_1.2.4   pbivnorm_0.6.0  openxlsx_4.1.0  rmarkdown_1.11 
[13] tools_3.5.1     stringr_1.3.1   foreign_0.8-71  yaml_2.2.0     
[17] xfun_0.4        compiler_3.5.1  mnormt_1.5-5    htmltools_0.3.6
[21] knitr_1.21     

Available under Created Commons license 3.0 CC BY