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).
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
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
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
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.
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.
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.
The first model is estimated with lavaan MLR.
A few highlights in our model code:
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.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
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.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
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
WLSMV
as the estimator. This is WLS with corrections for the standard errors and diagnostic statistics.mimic = "Mplus"
is used in our effort to synchronize these results with the estimates from Mplus in the companion exercise.missing = "listwise"
in the function call. That’s because listwise deletion is the current default.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
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.
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
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
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.
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
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
Factors: Bullying Victimization, Depression, and Alcohol Use
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
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 |
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
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 |
The direct effect of bullying on alcohol is compared with the indirect effect of bullying on depression, which also contributes to alcohol use.
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)
'
ML
. We are asking for bootstrapped standard errors, which are not available with MLR
.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
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 |
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 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