EXAMPLES

Brian O’Connor

2026-09-13

This R markdown document provides R code and output for multiple functions and datasets.

Load the EFA.dimensions package

library(EFA.dimensions)
**************************************************************************************************
EFA.dimensions 0.1.9.1

Please contact Brian O'Connor at brian.oconnor@ubc.ca if you have questions or suggestions.
**************************************************************************************************

ESEM examples:

ESEM, ‘efa_blocks’, no target, Holzinger & Swineford 1939 data

ESEM(data = data_HS_1939[,-10], method = 'efa_blocks', Nfactors = 3)


Exploratory Structural Equation Modeling:

The number of useable cases in data: 301

The number of variables: 9

method =  efa_blocks

The lavaan rotation method (rotation_LV): geomin 

Fit coefficients

 RMSR  = 0.02   GFI = 0.99   NFI  = 0.98   BIC = 7601.42   Chisq  = 22.90 
 RMSEA = 0.05   TLI = 0.96   NNFI = 0.96   AIC = 7479.08   df     = 12.00 
 SRMR  = 0.02   CFI = 0.99   IFI  = 0.99   MFI = 0.98       pvalue = 0.03

Standardized loadings

          LV_1     LV_2     LV_3
test1     0.61     0.17     0.01
test2     0.53     0.03    -0.14
test3     0.71    -0.09     0.00
test4     0.01     0.84     0.00
test5    -0.08     0.90     0.01
test6     0.07     0.81    -0.02
test7    -0.26     0.01     0.78
test8     0.01    -0.07     0.74
test9     0.31     0.00     0.48

Factor correlations

        LV_1    LV_2    LV_3
LV_1    1.00                
LV_2    0.37    1.00        
LV_3    0.43    0.31    1.00

ESEM, ‘startvalues’, no target, Holzinger & Swineford 1939 data

ESEM(data = data_HS_1939[,-10], method = 'startvalues', Nfactors = 3)


Exploratory Structural Equation Modeling:

The number of useable cases in data: 301

The number of variables: 9

method =  startvalues

The efa extraction method: ml

The efa rotation_EFA method: oblimin

The kind of correlations for the efa: pearson 

The number of factors:  3

The efa loadings:

      Factor_1 Factor_2 Factor_3
test1    0.031    0.191    0.602
test2   -0.117    0.044    0.505
test3    0.023   -0.069    0.689
test4    0.005    0.840    0.022
test5    0.008    0.888   -0.067
test6   -0.011    0.808    0.078
test7    0.723    0.044   -0.152
test8    0.701   -0.033    0.104
test9    0.463    0.035    0.366


The identified anchors:  test7 test5 test3 


The lavaan rotation method (rotation_LV): geomin 

Fit coefficients

 RMSR  = 0.11   GFI = 0.98   NFI  = 0.95   BIC = 7604.01   Chisq  = 42.61 
 RMSEA = 0.08   TLI = 0.92   NNFI = 0.92   AIC = 7492.80   df     = 15.00 
 SRMR  = 0.08   CFI = 0.97   IFI  = 0.97   MFI = 0.96       pvalue = 0.00

Standardized loadings

            F1       F2       F3
test1     0.02     0.18     0.61
test2    -0.11     0.04     0.50
test3     0.02    -0.06     0.63
test4    -0.01     0.84     0.01
test5     0.01     0.79    -0.06
test6    -0.02     0.79     0.07
test7     0.68     0.04    -0.14
test8     0.72    -0.05     0.11
test9     0.46     0.01     0.38

Factor correlations

        F1      F2      F3
F1    1.00                
F2    0.22    1.00        
F3    0.25    0.30    1.00

ESEM, ‘efa_blocks’, with a target, Holzinger & Swineford 1939 data

ESEM(data = data_HS_1939[,-10], method = 'efa_blocks', target_keys = c(1,1,1,2,2,2,3,3,3))


Exploratory Structural Equation Modeling:

The number of useable cases in data: 301

The number of variables: 9

method =  efa_blocks

The target loading matrix:

      Factor_1 Factor_2 Factor_3
test1        1        0        0
test2        1        0        0
test3        1        0        0
test4        0        1        0
test5        0        1        0
test6        0        1        0
test7        0        0        1
test8        0        0        1
test9        0        0        1


The lavaan rotation method (rotation_LV): target 

Fit coefficients

 RMSR  = 0.02   GFI = 0.99   NFI  = 0.98   BIC = 7601.42   Chisq  = 22.90 
 RMSEA = 0.05   TLI = 0.96   NNFI = 0.96   AIC = 7479.08   df     = 12.00 
 SRMR  = 0.02   CFI = 0.99   IFI  = 0.99   MFI = 0.98       pvalue = 0.03

Standardized loadings

          LV_1     LV_2     LV_3
test1     0.59     0.17     0.08
test2     0.51     0.03    -0.08
test3     0.68    -0.09     0.07
test4     0.02     0.84     0.01
test5    -0.07     0.89     0.00
test6     0.08     0.81     0.00
test7    -0.23     0.04     0.74
test8     0.03    -0.05     0.73
test9     0.31     0.02     0.50

Factor correlations

        LV_1    LV_2    LV_3
LV_1    1.00                
LV_2    0.34    1.00        
LV_3    0.31    0.26    1.00

ESEM, ‘startvalues’, with a target, Holzinger & Swineford 1939 data

ESEM(data = data_HS_1939[,-10], method = 'startvalues', target_keys = c(1,1,1,2,2,2,3,3,3))


Exploratory Structural Equation Modeling:

A target was provided and so rotation_EFA was set to "targetQ" for the EFA.

The number of useable cases in data: 301

The number of variables: 9

method =  startvalues

The efa extraction method: ml

The efa rotation_EFA method: targetQ  (oblique)

The kind of correlations for the efa: pearson 

Nfactors was not specified and so the EMPKC test was
conducted to determine the number of factors to extract: Nfactors =  3

The target loading matrix:

      Factor_1 Factor_2 Factor_3
test1       NA        0        0
test2       NA        0        0
test3       NA        0        0
test4        0       NA        0
test5        0       NA        0
test6        0       NA        0
test7        0        0       NA
test8        0        0       NA
test9        0        0       NA


The efa loadings:

      Factor_1 Factor_2 Factor_3
test1    0.589   -0.174   -0.075
test2    0.510   -0.032    0.084
test3    0.676    0.089   -0.073
test4    0.020   -0.840   -0.008
test5   -0.068   -0.890   -0.004
test6    0.077   -0.805    0.005
test7   -0.229   -0.036   -0.736
test8    0.025    0.047   -0.732
test9    0.309   -0.017   -0.505


The identified anchors:  test3 test5 test7 


The lavaan rotation method (rotation_LV): geomin 

Fit coefficients

 RMSR  = 0.11   GFI = 0.98   NFI  = 0.95   BIC = 7603.99   Chisq  = 42.59 
 RMSEA = 0.08   TLI = 0.92   NNFI = 0.92   AIC = 7492.78   df     = 15.00 
 SRMR  = 0.08   CFI = 0.97   IFI  = 0.97   MFI = 0.96       pvalue = 0.00

Standardized loadings

            F1       F2       F3
test1     0.60    -0.16    -0.07
test2     0.50    -0.03     0.08
test3     0.62     0.08    -0.07
test4     0.01    -0.84     0.01
test5    -0.06    -0.79     0.00
test6     0.07    -0.79     0.01
test7    -0.21    -0.03    -0.69
test8     0.03     0.07    -0.74
test9     0.32     0.01    -0.50

Factor correlations

         F1       F2       F3
F1     1.00                  
F2    -0.31     1.00         
F3    -0.29     0.26     1.00

ESEM, ‘efa_blocks’, with a target, Sanchez-Oliva et al. (2017) data

ESEM(data = data_SDT, method = 'efa_blocks', target_keys = c(1,1,1,1,2,2,2,2,3,3,3,3))


Exploratory Structural Equation Modeling:

The number of useable cases in data: 366

The number of variables: 12

method =  efa_blocks

The target loading matrix:

   Factor_1 Factor_2 Factor_3
R1        1        0        0
R2        1        0        0
R3        1        0        0
R4        1        0        0
C1        0        1        0
C2        0        1        0
C3        0        1        0
C4        0        1        0
A1        0        0        1
A2        0        0        1
A3        0        0        1
A4        0        0        1


The lavaan rotation method (rotation_LV): target 

Fit coefficients

 RMSR  = 0.04   GFI = 0.96   NFI  = 0.94   BIC = 10248.04   Chisq  = 120.41 
 RMSEA = 0.09   TLI = 0.91   NNFI = 0.91   AIC = 10072.42   df     = 33.00 
 SRMR  = 0.03   CFI = 0.95   IFI  = 0.96   MFI = 0.89       pvalue = 0.00

Standardized loadings

       LV_1     LV_2     LV_3
R1     0.80     0.06    -0.14
R2     0.82     0.03    -0.03
R3     0.69    -0.07     0.13
R4     0.60    -0.06     0.12
C1    -0.05     0.92    -0.08
C2    -0.02     0.89    -0.08
C3     0.08     0.18     0.67
C4     0.06     0.47     0.30
A1     0.30     0.24     0.18
A2     0.18     0.21     0.34
A3    -0.03     0.25     0.70
A4     0.28     0.10     0.38

Factor correlations

        LV_1    LV_2    LV_3
LV_1    1.00                
LV_2    0.55    1.00        
LV_3    0.46    0.49    1.00

ESEM, ‘efa_blocks’, MPlus 4.1b data

MPlus_4.1b_data <- read.table("https://www.statmodel.com/usersguide/chap4/ex4.1b.dat", header=FALSE)
colnames(MPlus_4.1b_data) <- paste0("y", 1:ncol(MPlus_4.1b_data))

ESEM(data = MPlus_4.1b_data, method = 'efa_blocks', Nfactors = 4, estimator = 'ML')


Exploratory Structural Equation Modeling:

The number of useable cases in data: 500

The number of variables: 12

method =  efa_blocks

The lavaan rotation method (rotation_LV): geomin 

Fit coefficients

 RMSR  = 0.01   GFI = 1.00   NFI  = 0.98   BIC = 16074.34   Chisq  = 25.80 
 RMSEA = 0.01   TLI = 1.00   NNFI = 1.00   AIC = 15846.75   df     = 24.00 
 SRMR  = 0.01   CFI = 1.00   IFI  = 1.00   MFI = 1.00       pvalue = 0.36

Standardized loadings

        LV_1     LV_2     LV_3     LV_4
y1      0.64     0.00     0.08     0.00
y2      0.81     0.01     0.00     0.06
y3      0.63    -0.05    -0.05    -0.01
y4      0.03     0.65     0.00    -0.02
y5     -0.03     0.76    -0.02     0.02
y6      0.01     0.67     0.03    -0.01
y7      0.00     0.00     0.73     0.02
y8     -0.03     0.00     0.73    -0.02
y9      0.05    -0.01     0.71     0.00
y10    -0.04     0.00    -0.01     0.69
y11     0.00     0.01     0.00     0.79
y12     0.03    -0.04     0.01     0.66

Factor correlations

         LV_1     LV_2     LV_3     LV_4
LV_1     1.00                           
LV_2    -0.02     1.00                  
LV_3    -0.01     0.02     1.00         
LV_4    -0.02    -0.12    -0.03     1.00
# compare with the MPlus output at  
# https://www.statmodel.com/usersguide/chap4/ex4.1part2.html

BIFACTOR examples:

BIFACTOR using ESEM, Holzinger & Swineford 1939 data

BIFACTOR(rawdata = data_HS_1939[,-10],
         bifactor_kind = 'ESEM',
         Nfactors = 4,
         LV_options = list(group_keys = NULL,
                           estimator = 'MLR',
                           rotation = 'bigeomin',
                           resid_correls = NULL,
                           LV_names = NULL,
                           ordered = FALSE))
Warning: lavaan->lav_object_post_check():  
   some estimated ov variances are negative


BIFACTOR output:

The number of factors =  4
The number of group factors =  3
The number of cases =  301 

bifactor_kind =  ESEM

LV_options group_keys = none provided

LV_options estimator =  MLR

LV_options rotation =  bigeomin

LV_options ordered =  FALSE

lavaan model syntax:

efa("BI-ESEM")*G + 
efa("BI-ESEM")*LV_1 + 
efa("BI-ESEM")*LV_2 + 
efa("BI-ESEM")*LV_3 =~ 
 test1 + test2 + test3 + test4 + test5 + test6 + test7 + test8 + test9

Fit coefficients

 RMSR  = 0.01   GFI = 1.00   NFI  = 0.99   BIC = 7618.15   Chisq  = 5.39 
 RMSEA = 0.00   TLI = 1.00   NNFI = 1.00   AIC = 7473.58   df     = 6.00 
 SRMR  = 0.01   CFI = 1.00   IFI  = 1.00   MFI = 1.00       pvalue = 0.49

Robust fit coefficients

 CFI robust   = 1.00   TLI robust = 1.00   NNFI robust = 1.00 
 RMSEA robust = 0.00   RNI robust = 1.00   GFI robust  = 1.00


WARNING: The bifactor loadings do not seem consistent with a bifactor structure:
         There are items that have more than two loadings that are >= the
         min_loading value of 0.2


Bifactor Loadings

        General   Group 1   Group 2   Group 3
test1      0.69      0.16      0.00     -0.09
test2      0.48      0.00     -0.14     -0.07
test3      0.63     -0.01      0.00     -0.16
test4      0.38      0.96      0.00      0.00
test5      0.44      0.00      0.00      1.23
test6      0.43      0.46      0.00      0.21
test7      0.09      0.14      0.71     -0.01
test8      0.33     -0.03      0.64      0.01
test9      0.54      0.00      0.41      0.00

Factor correlations

           General    Group 1    Group 2    Group 3
General          1       0.00       0.00       0.00
Group 1          0       1.00       0.01       0.48
Group 2          0       0.01       1.00      -0.01
Group 3          0       0.48      -0.01       1.00


Bifactor Model Statistics

   Omega total = 0.88 

   Omega hierarchical = 0.62

   Cronbach alpha = 0.76

   Average relative parameter bias (ARPB) = 0.72

   Percent of uncontaminated correlations (PUC) = NA

   Root mean square of the residuals using the general & group factors = 0.11

   Root mean square of the residuals using only the general factor = 0.2


(min_loading = 0.2, which is important for PUC and for omega and ECV group factor statistics)


Bifactor Factor Statistics
                        General   Group 1   Group 2   Group 3
   Omega total             0.88      0.85      0.73      1.05
   Omega hierarchical      0.62      0.64      0.57      0.77
   ECV - SS                0.38      0.75      0.72      0.72
   ECV - SG                0.34      0.19      0.19      0.26
   ECV - GS                0.38      0.22      0.27      0.18
   Factor Determinacy      0.88      1.02      0.83      1.33
   coefficient H           0.75      0.93      0.66      1.53


Bifactor Item Statistics

         Communalities    Uniquenesses     IECV    Rel. Param. Bias
test1            0.508           0.492    0.937               0.365
test2            0.256           0.744    0.911               0.544
test3            0.419           0.581    0.941               0.645
test4            1.069          -0.069    0.135               1.231
test5            1.699          -0.699    0.115               0.900
test6            0.444           0.556    0.424               0.930
test7            0.538           0.462    0.015               1.023
test8            0.524           0.476    0.206               0.387
test9            0.465           0.535    0.632               0.434



Eigenvalues and Proportions of Total Variance Explained

               Initial            Bifactor
           Eig  prop   cum   Eig    prop    cum 
Factor 1  3.22  0.36  0.36   2.04    0.23   0.23
Factor 2  1.64  0.18  0.54   1.22    0.14   0.36
Factor 3  1.37  0.15  0.69   1.12    0.12   0.49
Factor 4  0.70  0.08  0.77   1.63    0.18   0.67
Factor 5  0.58  0.06  0.83                      
Factor 6  0.50  0.06  0.89                      
Factor 7  0.47  0.05  0.94                      
Factor 8  0.29  0.03  0.97                      
Factor 9  0.24  0.03  1.00

BIFACTOR using CFA, Sanchez-Oliva et al. (2017) data

BIFACTOR(rawdata = data_SDT,   
         bifactor_kind = 'CFA', 
         LV_options = list(group_keys = c(1,1,1,1,2,2,2,2,3,3,3,3), 
                           estimator='MLR', 
                           resid_correls=NULL, 
                           LV_names=NULL,
                           ordered = FALSE))


BIFACTOR output:

The number of factors =  4
The number of group factors =  3
The number of cases =  366 

bifactor_kind =  CFA

LV_options group_keys:

R1 R2 R3 R4 C1 C2 C3 C4 A1 A2 A3 A4 
 1  1  1  1  2  2  2  2  3  3  3  3 


LV_options estimator =  MLR

LV_options ordered =  FALSE

lavaan model syntax:


Gen =~  R1 + R2 + R3 + R4 + C1 + C2 + C3 + C4 + A1 + A2 + A3 + A4
LV_1 =~  R1 + R2 + R3 + R4
LV_2 =~  C1 + C2 + C3 + C4
LV_3 =~  A1 + A2 + A3 + A4

Fit coefficients

 RMSR  = 0.05   GFI = 0.95   NFI  = 0.93   BIC = 10221.65   Chisq  = 147.15 
 RMSEA = 0.08   TLI = 0.91   NNFI = 0.91   AIC = 10081.16   df     = 42.00 
 SRMR  = 0.05   CFI = 0.95   IFI  = 0.95   MFI = 0.87       pvalue = 0.00

Robust fit coefficients

 CFI robust   = 0.94   TLI robust = 0.91   NNFI robust = 0.91 
 RMSEA robust = 0.08   RNI robust = 0.94   GFI robust  = 0.95


Bifactor Loadings

     General   Group 1   Group 2   Group 3
R1      0.43      0.64      0.00      0.00
R2      0.53      0.60      0.00      0.00
R3      0.51      0.52      0.00      0.00
R4      0.44      0.46      0.00      0.00
C1      0.57      0.00      0.71      0.00
C2      0.56      0.00      0.57      0.00
C3      0.80      0.00     -0.05      0.00
C4      0.67      0.00      0.24      0.00
A1      0.57      0.00      0.00      0.16
A2      0.61      0.00      0.00      0.46
A3      0.81      0.00      0.00     -0.17
A4      0.63      0.00      0.00      0.25

Factor correlations

           General    Group 1    Group 2    Group 3
General          1          0          0          0
Group 1          0          1          0          0
Group 2          0          0          1          0
Group 3          0          0          0          1


Bifactor Model Statistics

   Omega total = 0.92 

   Omega hierarchical = 0.8

   Cronbach alpha = 0.89

   Average relative parameter bias (ARPB) = 0.13

   Percent of uncontaminated correlations (PUC) = 0.85

   Root mean square of the residuals using the general & group factors = 0.05

   Root mean square of the residuals using only the general factor = 0.12


(min_loading = 0.2, which is important for PUC and for omega and ECV group factor statistics)


Bifactor Factor Statistics
                        General   Group 1   Group 2   Group 3
   Omega total             0.92      0.82      0.84      0.68
   Omega hierarchical      0.80      0.48      0.35      0.17
   ECV - SS                0.64      0.58      0.45      0.27
   ECV - SG                0.64      0.18      0.13      0.04
   ECV - GS                0.64      0.42      0.55      0.73
   Factor Determinacy      0.94      0.84      0.86      0.65
   coefficient H           0.89      0.66      0.60      0.29


Bifactor Item Statistics

      Communalities    Uniquenesses     IECV    Rel. Param. Bias
R1            0.595           0.405    0.304               0.346
R2            0.641           0.359    0.437               0.244
R3            0.529           0.471    0.492               0.194
R4            0.403           0.597    0.477               0.215
C1            0.822           0.178    0.394               0.129
C2            0.638           0.362    0.497               0.131
C3            0.640           0.360    0.996               0.104
C4            0.502           0.498    0.885               0.031
A1            0.351           0.649    0.926               0.061
A2            0.592           0.408    0.635               0.014
A3            0.681           0.319    0.959               0.140
A4            0.455           0.545    0.858               0.008



Eigenvalues and Proportions of Total Variance Explained

               Initial            Bifactor
            Eig  prop   cum   Eig    prop    cum 
Factor 1   5.42  0.45  0.45   4.38    0.37   0.37
Factor 2   1.39  0.12  0.57   1.25     0.1   0.47
Factor 3   0.98  0.08  0.65   0.88    0.07   0.54
Factor 4   0.83  0.07  0.72   0.33    0.03   0.57
Factor 5   0.63  0.05  0.77                      
Factor 6   0.60  0.05  0.82                      
Factor 7   0.45  0.04  0.86                      
Factor 8   0.43  0.04  0.90                      
Factor 9   0.39  0.03  0.93                      
Factor 10  0.33  0.03  0.96                      
Factor 11  0.27  0.02  0.98                      
Factor 12  0.26  0.02  1.00

BIFACTOR using ESEM, Sanchez-Oliva et al. (2017) data

BIFACTOR(rawdata = data_SDT,   
         bifactor_kind = 'ESEM', 
         EFA_options   = list(extraction = 'minres', rotation = 'oblimin', Nfactors = 3), 
         LV_options = list(group_keys = c(1,1,1,1,2,2,2,2,3,3,3,3), 
                           estimator='MLR', 
                           resid_correls=NULL, 
                           LV_names=NULL,
                           ordered = FALSE))


The names of one or more elements in LV_options is not valid.
The possibilities are: extraction and rotation

BIFACTOR output:

The number of factors =  4
The number of group factors =  3
The number of cases =  366 

bifactor_kind =  ESEM

LV_options group_keys:

R1 R2 R3 R4 C1 C2 C3 C4 A1 A2 A3 A4 
 1  1  1  1  2  2  2  2  3  3  3  3 


LV_options estimator =  MLR

LV_options ordered =  FALSE

lavaan model syntax:

efa("BI-ESEM")*G + 
efa("BI-ESEM")*LV_1 + 
efa("BI-ESEM")*LV_2 + 
efa("BI-ESEM")*LV_3 =~ 
 R1 + R2 + R3 + R4 + C1 + C2 + C3 + C4 + A1 + A2 + A3 + A4

Fit coefficients

 RMSR  = 0.02   GFI = 0.98   NFI  = 0.97   BIC = 10249.60   Chisq  = 68.85 
 RMSEA = 0.07   TLI = 0.94   NNFI = 0.94   AIC = 10038.85   df     = 24.00 
 SRMR  = 0.02   CFI = 0.98   IFI  = 0.98   MFI = 0.94       pvalue = 0.00

Robust fit coefficients

 CFI robust   = 0.98   TLI robust = 0.93   NNFI robust = 0.93 
 RMSEA robust = 0.07   RNI robust = 0.98   GFI robust  = 0.98


WARNING: The bifactor loadings do not seem consistent with a bifactor structure:
         There are items that have more than two loadings that are >= the
         min_loading value of 0.2


Bifactor Loadings

     General   Group 1   Group 2   Group 3
R1      0.44      0.60      0.10      0.10
R2      0.53      0.59      0.05      0.25
R3      0.54      0.52     -0.07     -0.09
R4      0.48      0.47     -0.05     -0.16
C1      0.59      0.01      0.67     -0.02
C2      0.57      0.02      0.57      0.07
C3      0.83     -0.07     -0.10     -0.09
C4      0.68      0.01      0.22     -0.01
A1      0.54      0.12      0.09      0.49
A2      0.59      0.04      0.02      0.18
A3      0.80     -0.15     -0.05      0.04
A4      0.60      0.10     -0.05      0.24

Factor correlations

           General    Group 1    Group 2    Group 3
General          1          0          0          0
Group 1          0          1          0          0
Group 2          0          0          1          0
Group 3          0          0          0          1


Bifactor Model Statistics

   Omega total = 0.92 

   Omega hierarchical = 0.8

   Cronbach alpha = 0.89

   Average relative parameter bias (ARPB) = 0.12

   Percent of uncontaminated correlations (PUC) = NA

   Root mean square of the residuals using the general & group factors = 0.03

   Root mean square of the residuals using only the general factor = 0.12


(min_loading = 0.2, which is important for PUC and for omega and ECV group factor statistics)


Bifactor Factor Statistics
                        General   Group 1   Group 2   Group 3
   Omega total             0.92      0.84      0.84      0.74
   Omega hierarchical      0.80      0.46      0.32      0.19
   ECV - SS                0.63      0.52      0.42      0.21
   ECV - SG                0.63      0.17      0.12      0.05
   ECV - GS                0.63      0.43      0.58      0.55
   Factor Determinacy      0.95      0.85      0.86      0.69
   coefficient H           0.90      0.65      0.58      0.35


Bifactor Item Statistics

      Communalities    Uniquenesses     IECV    Rel. Param. Bias
R1            0.575           0.425    0.337               0.302
R2            0.694           0.306    0.402               0.245
R3            0.584           0.416    0.508               0.117
R4            0.474           0.526    0.477               0.120
C1            0.803           0.197    0.436               0.085
C2            0.658           0.342    0.497               0.114
C3            0.704           0.296    0.969               0.134
C4            0.507           0.493    0.907               0.015
A1            0.551           0.449    0.524               0.125
A2            0.383           0.617    0.912               0.023
A3            0.670           0.330    0.958               0.133
A4            0.437           0.563    0.835               0.043



Eigenvalues and Proportions of Total Variance Explained

               Initial            Bifactor
            Eig  prop   cum   Eig    prop    cum 
Factor 1   5.42  0.45  0.45   4.46    0.37   0.37
Factor 2   1.39  0.12  0.57   1.26    0.11   0.48
Factor 3   0.98  0.08  0.65   0.87    0.07   0.55
Factor 4   0.83  0.07  0.72   0.45    0.04   0.59
Factor 5   0.63  0.05  0.77                      
Factor 6   0.60  0.05  0.82                      
Factor 7   0.45  0.04  0.86                      
Factor 8   0.43  0.04  0.90                      
Factor 9   0.39  0.03  0.93                      
Factor 10  0.33  0.03  0.96                      
Factor 11  0.27  0.02  0.98                      
Factor 12  0.26  0.02  1.00

BIFACTOR using ESEM, MPlus 4.7 data

MPlus_4.7_data <- read.table("https://www.statmodel.com/usersguide/chap4/ex4.7.dat", header=FALSE)
colnames(MPlus_4.7_data) <- paste0("y", 1:ncol(MPlus_4.7_data))

BIFACTOR(rawdata = MPlus_4.7_data,   
         Nfactors = 3,
         bifactor_kind = 'ESEM', 
         LV_options = list(group_keys = NULL,
                           estimator='ML',
                           rotation = 'bigeomin',
                           resid_correls=NULL,
                           LV_names=NULL,
                           ordered = FALSE))


BIFACTOR output:

The number of factors =  3
The number of group factors =  2
The number of cases =  5000 

bifactor_kind =  ESEM

LV_options group_keys = none provided

LV_options estimator =  ML

LV_options rotation =  bigeomin

LV_options ordered =  FALSE

lavaan model syntax:

efa("BI-ESEM")*G + 
efa("BI-ESEM")*LV_1 + 
efa("BI-ESEM")*LV_2 =~ 
 y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8 + y9 + y10

Fit coefficients

 RMSR  = 0.02   GFI = 1.00   NFI  = 1.00   BIC = 201183.13   Chisq  = 24.36 
 RMSEA = 0.01   TLI = 1.00   NNFI = 1.00   AIC = 200941.99   df     = 18.00 
 SRMR  = 0.01   CFI = 1.00   IFI  = 1.00   MFI = 1.00       pvalue = 0.14


Bifactor Loadings

      General   Group 1   Group 2
y1       0.47     -0.02      0.36
y2       0.45      0.00      0.12
y3       0.46     -0.01      0.26
y4       0.39      0.00      0.84
y5       0.38      0.22     -0.04
y6       0.30      0.43      0.00
y7       0.30      0.51      0.01
y8       0.32      0.57      0.00
y9       0.30      0.56     -0.03
y10      0.27      0.61      0.01

Factor correlations

           General    Group 1    Group 2
General          1       0.00       0.00
Group 1          0       1.00      -0.03
Group 2          0      -0.03       1.00


Bifactor Model Statistics

   Omega total = 0.79 

   Omega hierarchical = 0.44

   Cronbach alpha = 0.74

   Average relative parameter bias (ARPB) = 0.65

   Percent of uncontaminated correlations (PUC) = 0.6

   Root mean square of the residuals using the general & group factors = 0.01

   Root mean square of the residuals using only the general factor = 0.16


(min_loading = 0.2, which is important for PUC and for omega and ECV group factor statistics)


Bifactor Factor Statistics
                        General   Group 1   Group 2
   Omega total             0.79      0.75      0.72
   Omega hierarchical      0.44      0.53      0.40
   ECV - SS                0.36      0.72      0.61
   ECV - SG                0.36      0.40      0.24
   ECV - GS                0.36      0.28      0.39
   Factor Determinacy      0.72      0.80      0.87
   coefficient H           0.62      0.68      0.73


Bifactor Item Statistics

       Communalities    Uniquenesses     IECV    Rel. Param. Bias
y1             0.353           0.647    0.622               0.376
y2             0.215           0.785    0.928               0.361
y3             0.279           0.721    0.759               0.358
y4             0.862           0.138    0.175               0.296
y5             0.196           0.804    0.750               0.045
y6             0.278           0.722    0.331               0.736
y7             0.347           0.653    0.256               0.956
y8             0.432           0.568    0.234               1.022
y9             0.405           0.595    0.223               1.041
y10            0.448           0.552    0.168               1.318



Eigenvalues and Proportions of Total Variance Explained

               Initial            Bifactor
            Eig  prop   cum   Eig    prop    cum 
Factor 1   3.03  0.30  0.30   1.38    0.14   0.14
Factor 2   1.67  0.17  0.47   1.51    0.15   0.29
Factor 3   0.85  0.08  0.55   0.93    0.09   0.38
Factor 4   0.79  0.08  0.63                      
Factor 5   0.71  0.07  0.70                      
Factor 6   0.70  0.07  0.77                      
Factor 7   0.64  0.06  0.84                      
Factor 8   0.58  0.06  0.90                      
Factor 9   0.56  0.06  0.95                      
Factor 10  0.49  0.05  1.00                      
# compare with the MPlus output at  
# hhttps://www.statmodel.com/usersguide/chap4/ex4.7.html

Factorial_Invariance examples:

using LV_model, Holzinger & Swineford 1939 data

HS_1939_model <- 'visual  =~ test1 + test2 + test3
                  textual =~ test4 + test5 + test6
                  speed   =~ test7 + test8 + test9'

HS_1939_output <-
  Factorial_Invariance(data = data_HS_1939, group = 'school', LV_model = HS_1939_model)


Factorial_Invariance

The model:
visual  =~ test1 + test2 + test3
                  textual =~ test4 + test5 + test6
                  speed   =~ test7 + test8 + test9 


The levels of "group" and their Ns:
Grant-White     Pasteur 
        145         156 


Fit coefficients for each group separately:
               rmsea     srmr      cfi      tli         aic         bic
Pasteur        0.104    0.077    0.903    0.855    3936.617    4000.664
Grant-White    0.089    0.072    0.941    0.912    3511.778    3574.289


Invariance Model Fit Coefficients:
                  rmsea     srmr      cfi      tli         aic         bic
mod_Configural    0.097    0.068    0.923    0.885    7484.395    7706.822
mod_Metric        0.093    0.072    0.921    0.895    7480.587    7680.771
mod_Scalar        0.107    0.082    0.882    0.859    7508.647    7686.588
mod_Strict        0.104    0.088    0.873    0.867    7508.055    7652.632
mod_LV_vars       0.102    0.091    0.873    0.873    7504.782    7638.238
mod_LV_covars     0.100    0.095    0.873    0.878    7501.931    7624.266
mod_LV_means      0.113    0.119    0.832    0.845    7535.490    7646.703


Differences in Invariance Model Fit Coefficients:
                                 rmsea     srmr       cfi       tli       aic
mod_Metric - mod_Configural     -0.004    0.004    -0.002     0.009    -3.808
mod_Scalar - mod_Metric          0.015    0.011    -0.038    -0.036    28.059
mod_Strict - mod_Scalar         -0.003    0.006    -0.009     0.008    -0.591
mod_LV_vars - mod_Strict        -0.002    0.003     0.000     0.006    -3.273
mod_LV_covars - mod_LV_vars     -0.002    0.003     0.000     0.005    -2.851
mod_LV_means - mod_LV_covars     0.013    0.024    -0.041    -0.033    33.559
                                    bic
mod_Metric - mod_Configural     -26.050
mod_Scalar - mod_Metric           5.817
mod_Strict - mod_Scalar         -33.955
mod_LV_vars - mod_Strict        -14.395
mod_LV_covars - mod_LV_vars     -13.972
mod_LV_means - mod_LV_covars     22.437


Chi-Squared Difference Tests:
                  Df       AIC       BIC     Chisq    Chisq diff    RMSEA
mod_Configural    48    7484.4    7706.8    115.85                       
mod_Metric        54    7480.6    7680.8    124.04         8.192    0.049
mod_Scalar        60    7508.6    7686.6    164.10        40.059    0.194
mod_Strict        69    7508.1    7652.6    181.51        17.409    0.079
mod_LV_vars       72    7504.8    7638.2    184.24         2.727    0.000
mod_LV_covars     75    7501.9    7624.3    187.39         3.149    0.018
mod_LV_means      78    7535.5    7646.7    226.95        39.559    0.285
                  Df diff    Pr(>Chisq)       
mod_Configural                                
mod_Metric              6         0.224       
mod_Scalar              6        <2e-16    ***
mod_Strict              9         0.043    *  
mod_LV_vars             3         0.436       
mod_LV_covars           3         0.369       
mod_LV_means            3        <2e-16    ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Coefficients for ConfiguralStandardized loadings for each group (Configural):

$Pasteur
         visual    textul    speed
test1      0.89      0.00     0.00
test2      0.34      0.00     0.00
test3      0.51      0.00     0.00
test4      0.00      0.82     0.00
test5      0.00      0.86     0.00
test6      0.00      0.84     0.00
test7      0.00      0.00     0.55
test8      0.00      0.00     0.68
test9      0.00      0.00     0.55

$`Grant-White`
         visual    textul    speed
test1      0.68      0.00     0.00
test2      0.52      0.00     0.00
test3      0.69      0.00     0.00
test4      0.00      0.87     0.00
test5      0.00      0.83     0.00
test6      0.00      0.83     0.00
test7      0.00      0.00     0.66
test8      0.00      0.00     0.80
test9      0.00      0.00     0.70



WARNING: The latent trait means for the groups are apparently the same.
This occurs in a model where the intercepts & loadings were not constrained to be equal.
Consider using a "Scalar" model instead.
Latent variable group means (Configural):

           Pasteur    Grant-White
visual           0              0
textual          0              0
speed            0              0

Latent variable group standard deviations (Configural):

           Pasteur    Grant-White
visual        0.95           0.67
textual       0.89           0.92
speed         0.48           0.61

Latent variable group Ns (Configural):

           Pasteur    Grant-White
visual         156            145
textual        156            145
speed          156            145


Coefficients for ScalarStandardized loadings for each group (Scalar):

$Pasteur
         visual    textul    speed
test1      0.77      0.00     0.00
test2      0.41      0.00     0.00
test3      0.59      0.00     0.00
test4      0.00      0.81     0.00
test5      0.00      0.83     0.00
test6      0.00      0.86     0.00
test7      0.00      0.00     0.52
test8      0.00      0.00     0.66
test9      0.00      0.00     0.58

$`Grant-White`
         visual    textul    speed
test1      0.72      0.00     0.00
test2      0.44      0.00     0.00
test3      0.64      0.00     0.00
test4      0.00      0.85     0.00
test5      0.00      0.86     0.00
test6      0.00      0.80     0.00
test7      0.00      0.00     0.67
test8      0.00      0.00     0.77
test9      0.00      0.00     0.70

Latent variable group means (Scalar):

           Pasteur    Grant-White
visual           0          -0.15
textual          0           0.58
speed            0          -0.18

Latent variable group standard deviations (Scalar):

           Pasteur    Grant-White
visual        0.76           0.72
textual       0.88           0.88
speed         0.46           0.63

Latent variable group Ns (Scalar):

           Pasteur    Grant-White
visual         156            145
textual        156            145
speed          156            145
PLOT_Invariance(model_object = HS_1939_output,
                invar_model = 'Scalar',
                plot_types = c('loadings', 'LV_distribs')) 



Latent variable means:

           Pasteur    Grant-White
visual           0        -0.1477
textual          0         0.5764
speed            0        -0.1774
# for more detailed statistical comparisons of the latent variable means, try
# the GROUP.DIFFS function from the DFA.CANCOR package, as follows:
# install.packages(DFA.CANCOR); library(DFA.CANCOR)
# 
# LV_scores <- HS_1939_output$LV_scores
# 
# GROUP.DIFFS(data = LV_scores, GROUPS = 'school', DV = 'visual')
# 
# GROUP.DIFFS(data = LV_scores, GROUPS = 'school', DV = 'speed')

using LV_keys, Holzinger & Swineford 1939 data

HS_1939_output <-
Factorial_Invariance(data = data_HS_1939, group = 'school',
                     LV_keys = c(1,1,1, 2,2,2, 3,3,3),
                     LV_names = c('visual', 'textual', 'speed'))


Factorial_Invariance

The model:

visual =~  test1 + test2 + test3
textual =~  test4 + test5 + test6
speed =~  test7 + test8 + test9 


The levels of "group" and their Ns:
Grant-White     Pasteur 
        145         156 


Fit coefficients for each group separately:
               rmsea     srmr      cfi      tli         aic         bic
Pasteur        0.104    0.077    0.903    0.855    3936.617    4000.664
Grant-White    0.089    0.072    0.941    0.912    3511.778    3574.289


Invariance Model Fit Coefficients:
                  rmsea     srmr      cfi      tli         aic         bic
mod_Configural    0.097    0.068    0.923    0.885    7484.395    7706.822
mod_Metric        0.093    0.072    0.921    0.895    7480.587    7680.771
mod_Scalar        0.107    0.082    0.882    0.859    7508.647    7686.588
mod_Strict        0.104    0.088    0.873    0.867    7508.055    7652.632
mod_LV_vars       0.102    0.091    0.873    0.873    7504.782    7638.238
mod_LV_covars     0.100    0.095    0.873    0.878    7501.931    7624.266
mod_LV_means      0.113    0.119    0.832    0.845    7535.490    7646.703


Differences in Invariance Model Fit Coefficients:
                                 rmsea     srmr       cfi       tli       aic
mod_Metric - mod_Configural     -0.004    0.004    -0.002     0.009    -3.808
mod_Scalar - mod_Metric          0.015    0.011    -0.038    -0.036    28.059
mod_Strict - mod_Scalar         -0.003    0.006    -0.009     0.008    -0.591
mod_LV_vars - mod_Strict        -0.002    0.003     0.000     0.006    -3.273
mod_LV_covars - mod_LV_vars     -0.002    0.003     0.000     0.005    -2.851
mod_LV_means - mod_LV_covars     0.013    0.024    -0.041    -0.033    33.559
                                    bic
mod_Metric - mod_Configural     -26.050
mod_Scalar - mod_Metric           5.817
mod_Strict - mod_Scalar         -33.955
mod_LV_vars - mod_Strict        -14.395
mod_LV_covars - mod_LV_vars     -13.972
mod_LV_means - mod_LV_covars     22.437


Chi-Squared Difference Tests:
                  Df       AIC       BIC     Chisq    Chisq diff    RMSEA
mod_Configural    48    7484.4    7706.8    115.85                       
mod_Metric        54    7480.6    7680.8    124.04         8.192    0.049
mod_Scalar        60    7508.6    7686.6    164.10        40.059    0.194
mod_Strict        69    7508.1    7652.6    181.51        17.409    0.079
mod_LV_vars       72    7504.8    7638.2    184.24         2.727    0.000
mod_LV_covars     75    7501.9    7624.3    187.39         3.149    0.018
mod_LV_means      78    7535.5    7646.7    226.95        39.559    0.285
                  Df diff    Pr(>Chisq)       
mod_Configural                                
mod_Metric              6         0.224       
mod_Scalar              6        <2e-16    ***
mod_Strict              9         0.043    *  
mod_LV_vars             3         0.436       
mod_LV_covars           3         0.369       
mod_LV_means            3        <2e-16    ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Coefficients for ConfiguralStandardized loadings for each group (Configural):

$Pasteur
         visual    textul    speed
test1      0.89      0.00     0.00
test2      0.34      0.00     0.00
test3      0.51      0.00     0.00
test4      0.00      0.82     0.00
test5      0.00      0.86     0.00
test6      0.00      0.84     0.00
test7      0.00      0.00     0.55
test8      0.00      0.00     0.68
test9      0.00      0.00     0.55

$`Grant-White`
         visual    textul    speed
test1      0.68      0.00     0.00
test2      0.52      0.00     0.00
test3      0.69      0.00     0.00
test4      0.00      0.87     0.00
test5      0.00      0.83     0.00
test6      0.00      0.83     0.00
test7      0.00      0.00     0.66
test8      0.00      0.00     0.80
test9      0.00      0.00     0.70



WARNING: The latent trait means for the groups are apparently the same.
This occurs in a model where the intercepts & loadings were not constrained to be equal.
Consider using a "Scalar" model instead.
Latent variable group means (Configural):

           Pasteur    Grant-White
visual           0              0
textual          0              0
speed            0              0

Latent variable group standard deviations (Configural):

           Pasteur    Grant-White
visual        0.95           0.67
textual       0.89           0.92
speed         0.48           0.61

Latent variable group Ns (Configural):

           Pasteur    Grant-White
visual         156            145
textual        156            145
speed          156            145


Coefficients for ScalarStandardized loadings for each group (Scalar):

$Pasteur
         visual    textul    speed
test1      0.77      0.00     0.00
test2      0.41      0.00     0.00
test3      0.59      0.00     0.00
test4      0.00      0.81     0.00
test5      0.00      0.83     0.00
test6      0.00      0.86     0.00
test7      0.00      0.00     0.52
test8      0.00      0.00     0.66
test9      0.00      0.00     0.58

$`Grant-White`
         visual    textul    speed
test1      0.72      0.00     0.00
test2      0.44      0.00     0.00
test3      0.64      0.00     0.00
test4      0.00      0.85     0.00
test5      0.00      0.86     0.00
test6      0.00      0.80     0.00
test7      0.00      0.00     0.67
test8      0.00      0.00     0.77
test9      0.00      0.00     0.70

Latent variable group means (Scalar):

           Pasteur    Grant-White
visual           0          -0.15
textual          0           0.58
speed            0          -0.18

Latent variable group standard deviations (Scalar):

           Pasteur    Grant-White
visual        0.76           0.72
textual       0.88           0.88
speed         0.46           0.63

Latent variable group Ns (Scalar):

           Pasteur    Grant-White
visual         156            145
textual        156            145
speed          156            145
PLOT_Invariance(model_object = HS_1939_output, 
                invar_model = 'Configural',
                plot_types = c('loadings', 'ints_slopes'))

another way of using LV_keys, Holzinger & Swineford 1939 data

Factorial_Invariance(data = data_HS_1939, group = 'school',
                     LV_keys = c(test1 = 1, test2 = 1, test3 = 1,
                                 test4 = 2, test5 = 2, test6 = 2,
                                 test7 = 3, test8 = 3, test9 = 3) )


Factorial_Invariance

The model:

LV_1 =~  test1 + test2 + test3
LV_2 =~  test4 + test5 + test6
LV_3 =~  test7 + test8 + test9 


The levels of "group" and their Ns:
Grant-White     Pasteur 
        145         156 


Fit coefficients for each group separately:
               rmsea     srmr      cfi      tli         aic         bic
Pasteur        0.104    0.077    0.903    0.855    3936.617    4000.664
Grant-White    0.089    0.072    0.941    0.912    3511.778    3574.289


Invariance Model Fit Coefficients:
                  rmsea     srmr      cfi      tli         aic         bic
mod_Configural    0.097    0.068    0.923    0.885    7484.395    7706.822
mod_Metric        0.093    0.072    0.921    0.895    7480.587    7680.771
mod_Scalar        0.107    0.082    0.882    0.859    7508.647    7686.588
mod_Strict        0.104    0.088    0.873    0.867    7508.055    7652.632
mod_LV_vars       0.102    0.091    0.873    0.873    7504.782    7638.238
mod_LV_covars     0.100    0.095    0.873    0.878    7501.931    7624.266
mod_LV_means      0.113    0.119    0.832    0.845    7535.490    7646.703


Differences in Invariance Model Fit Coefficients:
                                 rmsea     srmr       cfi       tli       aic
mod_Metric - mod_Configural     -0.004    0.004    -0.002     0.009    -3.808
mod_Scalar - mod_Metric          0.015    0.011    -0.038    -0.036    28.059
mod_Strict - mod_Scalar         -0.003    0.006    -0.009     0.008    -0.591
mod_LV_vars - mod_Strict        -0.002    0.003     0.000     0.006    -3.273
mod_LV_covars - mod_LV_vars     -0.002    0.003     0.000     0.005    -2.851
mod_LV_means - mod_LV_covars     0.013    0.024    -0.041    -0.033    33.559
                                    bic
mod_Metric - mod_Configural     -26.050
mod_Scalar - mod_Metric           5.817
mod_Strict - mod_Scalar         -33.955
mod_LV_vars - mod_Strict        -14.395
mod_LV_covars - mod_LV_vars     -13.972
mod_LV_means - mod_LV_covars     22.437


Chi-Squared Difference Tests:
                  Df       AIC       BIC     Chisq    Chisq diff    RMSEA
mod_Configural    48    7484.4    7706.8    115.85                       
mod_Metric        54    7480.6    7680.8    124.04         8.192    0.049
mod_Scalar        60    7508.6    7686.6    164.10        40.059    0.194
mod_Strict        69    7508.1    7652.6    181.51        17.409    0.079
mod_LV_vars       72    7504.8    7638.2    184.24         2.727    0.000
mod_LV_covars     75    7501.9    7624.3    187.39         3.149    0.018
mod_LV_means      78    7535.5    7646.7    226.95        39.559    0.285
                  Df diff    Pr(>Chisq)       
mod_Configural                                
mod_Metric              6         0.224       
mod_Scalar              6        <2e-16    ***
mod_Strict              9         0.043    *  
mod_LV_vars             3         0.436       
mod_LV_covars           3         0.369       
mod_LV_means            3        <2e-16    ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Coefficients for ConfiguralStandardized loadings for each group (Configural):

$Pasteur
         LV_1    LV_2    LV_3
test1    0.89    0.00    0.00
test2    0.34    0.00    0.00
test3    0.51    0.00    0.00
test4    0.00    0.82    0.00
test5    0.00    0.86    0.00
test6    0.00    0.84    0.00
test7    0.00    0.00    0.55
test8    0.00    0.00    0.68
test9    0.00    0.00    0.55

$`Grant-White`
         LV_1    LV_2    LV_3
test1    0.68    0.00    0.00
test2    0.52    0.00    0.00
test3    0.69    0.00    0.00
test4    0.00    0.87    0.00
test5    0.00    0.83    0.00
test6    0.00    0.83    0.00
test7    0.00    0.00    0.66
test8    0.00    0.00    0.80
test9    0.00    0.00    0.70



WARNING: The latent trait means for the groups are apparently the same.
This occurs in a model where the intercepts & loadings were not constrained to be equal.
Consider using a "Scalar" model instead.
Latent variable group means (Configural):

        Pasteur    Grant-White
LV_1          0              0
LV_2          0              0
LV_3          0              0

Latent variable group standard deviations (Configural):

        Pasteur    Grant-White
LV_1       0.95           0.67
LV_2       0.89           0.92
LV_3       0.48           0.61

Latent variable group Ns (Configural):

        Pasteur    Grant-White
LV_1        156            145
LV_2        156            145
LV_3        156            145


Coefficients for ScalarStandardized loadings for each group (Scalar):

$Pasteur
         LV_1    LV_2    LV_3
test1    0.77    0.00    0.00
test2    0.41    0.00    0.00
test3    0.59    0.00    0.00
test4    0.00    0.81    0.00
test5    0.00    0.83    0.00
test6    0.00    0.86    0.00
test7    0.00    0.00    0.52
test8    0.00    0.00    0.66
test9    0.00    0.00    0.58

$`Grant-White`
         LV_1    LV_2    LV_3
test1    0.72    0.00    0.00
test2    0.44    0.00    0.00
test3    0.64    0.00    0.00
test4    0.00    0.85    0.00
test5    0.00    0.86    0.00
test6    0.00    0.80    0.00
test7    0.00    0.00    0.67
test8    0.00    0.00    0.77
test9    0.00    0.00    0.70

Latent variable group means (Scalar):

        Pasteur    Grant-White
LV_1          0          -0.15
LV_2          0           0.58
LV_3          0          -0.18

Latent variable group standard deviations (Scalar):

        Pasteur    Grant-White
LV_1       0.76           0.72
LV_2       0.88           0.88
LV_3       0.46           0.63

Latent variable group Ns (Scalar):

        Pasteur    Grant-White
LV_1        156            145
LV_2        156            145
LV_3        156            145

with a correlated residuals term, Brown 2015 data

# 2015 Brown - Confirmatory Factor Analysis for Applied Research  p 247, Table 7.11

# Tests of measurement invariance and population heterogeneity of DSM-IV
# major depressive disorder in men and women
data_Brown_2015 <- read.table("http://people.bu.edu/tabrown/Ch7/MDDALL.dat")
names(data_Brown_2015) <- c("sex", paste("mdd", 1:9, sep = ""))
data_Brown_2015$sex <- factor(data_Brown_2015$sex, levels = c(0, 1), labels = c("female", "male"))

# using LV_model
model.mdd <- '
MDD =~ mdd1 + mdd2 + mdd3 + mdd4 + mdd5 + mdd6 + mdd7 + mdd8 + mdd9
mdd1 ~~ mdd2
'
Factorial_Invariance(data = data_Brown_2015, group = 'sex', LV_model = model.mdd)


These variables have 10 or fewer levels/values:

mdd1 mdd2 mdd3 mdd4 mdd5 mdd6 mdd7 mdd8 mdd9 
   8    8    8    9    8    9    9    9    8 

Consider changing ordered to TRUE in such cases.



Factorial_Invariance

The model:

MDD =~ mdd1 + mdd2 + mdd3 + mdd4 + mdd5 + mdd6 + mdd7 + mdd8 + mdd9
mdd1 ~~ mdd2
 


The levels of "group" and their Ns:
female   male 
   375    375 


Fit coefficients for each group separately:
          rmsea     srmr      cfi      tli         aic         bic
female    0.053    0.044    0.959    0.944    13686.62    13761.23
male      0.045    0.041    0.967    0.955    13803.17    13877.79


Invariance Model Fit Coefficients:
                  rmsea     srmr      cfi      tli         aic         bic
mod_Configural    0.049    0.039    0.963    0.949    27525.80    27784.52
mod_Metric        0.044    0.042    0.966    0.960    27513.72    27735.49
mod_Scalar        0.043    0.045    0.963    0.961    27510.19    27695.00
mod_Strict        0.041    0.048    0.962    0.965    27501.91    27645.13
mod_LV_vars       0.040    0.050    0.962    0.965    27500.70    27639.30
mod_LV_means      0.041    0.051    0.962    0.965    27500.62    27634.60


Differences in Invariance Model Fit Coefficients:
                                rmsea     srmr       cfi      tli        aic
mod_Metric - mod_Configural    -0.005    0.003     0.003    0.011    -12.071
mod_Scalar - mod_Metric        -0.001    0.003    -0.004    0.001     -3.530
mod_Strict - mod_Scalar        -0.002    0.004    -0.001    0.004     -8.288
mod_LV_vars - mod_Strict        0.000    0.002     0.000    0.001     -1.207
mod_LV_means - mod_LV_vars      0.000    0.002    -0.001    0.000     -0.080
                                   bic
mod_Metric - mod_Configural    -49.032
mod_Scalar - mod_Metric        -40.491
mod_Strict - mod_Scalar        -49.869
mod_LV_vars - mod_Strict        -5.827
mod_LV_means - mod_LV_vars      -4.700


Chi-Squared Difference Tests:
                  Df      AIC      BIC      Chisq    Chisq diff    RMSEA
mod_Configural    52    27526    27784     98.911                       
mod_Metric        60    27514    27736    102.839         3.929    0.000
mod_Scalar        68    27510    27695    115.309        12.470    0.039
mod_Strict        77    27502    27645    125.021         9.712    0.015
mod_LV_vars       78    27501    27639    125.814         0.793    0.000
mod_LV_means      79    27501    27635    127.734         1.920    0.050
                  Df diff    Pr(>Chisq)
mod_Configural                         
mod_Metric              8         0.864
mod_Scalar              8         0.131
mod_Strict              9         0.374
mod_LV_vars             1         0.373
mod_LV_means            1         0.166


Coefficients for ConfiguralStandardized loadings for each group (Configural):

$female
         MDD
mdd1    0.73
mdd2    0.69
mdd3    0.44
mdd4    0.52
mdd5    0.48
mdd6    0.58
mdd7    0.37
mdd8    0.57
mdd9    0.34

$male
         MDD
mdd1    0.64
mdd2    0.63
mdd3    0.38
mdd4    0.54
mdd5    0.47
mdd6    0.55
mdd7    0.34
mdd8    0.48
mdd9    0.36



WARNING: The latent trait means for the groups are apparently the same.
This occurs in a model where the intercepts & loadings were not constrained to be equal.
Consider using a "Scalar" model instead.
Latent variable group means (Configural):

       female    male
MDD         0       0

Latent variable group standard deviations (Configural):

       female    male
MDD      1.11    0.88

Latent variable group Ns (Configural):

       female    male
MDD       375     375


Coefficients for ScalarStandardized loadings for each group (Scalar):

$female
         MDD
mdd1    0.70
mdd2    0.68
mdd3    0.43
mdd4    0.54
mdd5    0.48
mdd6    0.59
mdd7    0.39
mdd8    0.55
mdd9    0.36

$male
         MDD
mdd1    0.67
mdd2    0.63
mdd3    0.40
mdd4    0.51
mdd5    0.46
mdd6    0.53
mdd7    0.33
mdd8    0.50
mdd9    0.34

Latent variable group means (Scalar):

       female     male
MDD         0    -0.13

Latent variable group standard deviations (Scalar):

       female    male
MDD      1.04    0.94

Latent variable group Ns (Scalar):

       female    male
MDD       375     375
# using LV_keys
Brown_2015_output <-
  Factorial_Invariance(data = data_Brown_2015, group = 'sex', 
                     LV_keys = c(1,1,1,1,1,1,1,1,1),
                     LV_resid_correls = c('mdd1 ~~ mdd2'))


These variables have 10 or fewer levels/values:

mdd1 mdd2 mdd3 mdd4 mdd5 mdd6 mdd7 mdd8 mdd9 
   8    8    8    9    8    9    9    9    8 

Consider changing ordered to TRUE in such cases.



Factorial_Invariance

The model:

LV_1 =~  mdd1 + mdd2 + mdd3 + mdd4 + mdd5 + mdd6 + mdd7 + mdd8 + mdd9
mdd1 ~~ mdd2 


The levels of "group" and their Ns:
female   male 
   375    375 


Fit coefficients for each group separately:
          rmsea     srmr      cfi      tli         aic         bic
female    0.053    0.044    0.959    0.944    13686.62    13761.23
male      0.045    0.041    0.967    0.955    13803.17    13877.79


Invariance Model Fit Coefficients:
                  rmsea     srmr      cfi      tli         aic         bic
mod_Configural    0.049    0.039    0.963    0.949    27525.80    27784.52
mod_Metric        0.044    0.042    0.966    0.960    27513.72    27735.49
mod_Scalar        0.043    0.045    0.963    0.961    27510.19    27695.00
mod_Strict        0.041    0.048    0.962    0.965    27501.91    27645.13
mod_LV_vars       0.040    0.050    0.962    0.965    27500.70    27639.30
mod_LV_means      0.041    0.051    0.962    0.965    27500.62    27634.60


Differences in Invariance Model Fit Coefficients:
                                rmsea     srmr       cfi      tli        aic
mod_Metric - mod_Configural    -0.005    0.003     0.003    0.011    -12.071
mod_Scalar - mod_Metric        -0.001    0.003    -0.004    0.001     -3.530
mod_Strict - mod_Scalar        -0.002    0.004    -0.001    0.004     -8.288
mod_LV_vars - mod_Strict        0.000    0.002     0.000    0.001     -1.207
mod_LV_means - mod_LV_vars      0.000    0.002    -0.001    0.000     -0.080
                                   bic
mod_Metric - mod_Configural    -49.032
mod_Scalar - mod_Metric        -40.491
mod_Strict - mod_Scalar        -49.869
mod_LV_vars - mod_Strict        -5.827
mod_LV_means - mod_LV_vars      -4.700


Chi-Squared Difference Tests:
                  Df      AIC      BIC      Chisq    Chisq diff    RMSEA
mod_Configural    52    27526    27784     98.911                       
mod_Metric        60    27514    27736    102.839         3.929    0.000
mod_Scalar        68    27510    27695    115.309        12.470    0.039
mod_Strict        77    27502    27645    125.021         9.712    0.015
mod_LV_vars       78    27501    27639    125.814         0.793    0.000
mod_LV_means      79    27501    27635    127.734         1.920    0.050
                  Df diff    Pr(>Chisq)
mod_Configural                         
mod_Metric              8         0.864
mod_Scalar              8         0.131
mod_Strict              9         0.374
mod_LV_vars             1         0.373
mod_LV_means            1         0.166


Coefficients for ConfiguralStandardized loadings for each group (Configural):

$female
        LV_1
mdd1    0.73
mdd2    0.69
mdd3    0.44
mdd4    0.52
mdd5    0.48
mdd6    0.58
mdd7    0.37
mdd8    0.57
mdd9    0.34

$male
        LV_1
mdd1    0.64
mdd2    0.63
mdd3    0.38
mdd4    0.54
mdd5    0.47
mdd6    0.55
mdd7    0.34
mdd8    0.48
mdd9    0.36



WARNING: The latent trait means for the groups are apparently the same.
This occurs in a model where the intercepts & loadings were not constrained to be equal.
Consider using a "Scalar" model instead.
Latent variable group means (Configural):

        female    male
LV_1         0       0

Latent variable group standard deviations (Configural):

        female    male
LV_1      1.11    0.88

Latent variable group Ns (Configural):

        female    male
LV_1       375     375


Coefficients for ScalarStandardized loadings for each group (Scalar):

$female
        LV_1
mdd1    0.70
mdd2    0.68
mdd3    0.43
mdd4    0.54
mdd5    0.48
mdd6    0.59
mdd7    0.39
mdd8    0.55
mdd9    0.36

$male
        LV_1
mdd1    0.67
mdd2    0.63
mdd3    0.40
mdd4    0.51
mdd5    0.46
mdd6    0.53
mdd7    0.33
mdd8    0.50
mdd9    0.34

Latent variable group means (Scalar):

        female     male
LV_1         0    -0.13

Latent variable group standard deviations (Scalar):

        female    male
LV_1      1.04    0.94

Latent variable group Ns (Scalar):

        female    male
LV_1       375     375
PLOT_Invariance(model_object = Brown_2015_output,
                invar_model = 'Scalar',
                plot_types = c('loadings', 'LV_distribs')) 



Latent variable means:

        female       male
LV_1         0    -0.1308