Chapter 16:
Phil Wood
2025-08-10
Reasons for Drinking (Simulated) data
Load RFD Data
<a https://docs.google.com/document/d/1_cDpW3tD7eCP6Z_mFz4Z4u0vY5Nx5z4QK7D8sAwWHAM/edit?usp=sharing
Read in an example (simulated) data set for analysis
Two-Factor Unrotated Solution (Table 16.1)
To begin the presentation, this is simply an example of how one might do an exploratory factor analysis and the output which results. These two approaches yield the same values although in practice minor differences may emerge between the two approaches having to do with when the rotation is specified relative to extraction. In this analysis F1, F2, F3, and F4 from fa are represented as ML1, ML3, ML4, and ML2 respectively.
Using fa()
# Packages
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(psych)
library(GPArotation) ##
## Attaching package: 'GPArotation'## The following objects are masked from 'package:psych':
##
## equamax, varimin# Perform exploratory factor analysis with 2 factors, unrotated (maximum likelihood method)
fa_result <- factanal(RFD, factors = 2, rotation = "none")
# View the results
print(fa_result)##
## Call:
## factanal(x = RFD, factors = 2, rotation = "none")
##
## Uniquenesses:
## rridaw0 rridbw0 rridcw0 rriddw0 rridew0 rridsw0 rridtw0 rriduw0 rridvw0 rridww0
## 0.308 0.506 0.119 0.152 0.566 0.475 0.353 0.188 0.158 0.373
##
## Loadings:
## Factor1 Factor2
## rridaw0 0.744 -0.372
## rridbw0 0.664 -0.229
## rridcw0 0.822 -0.453
## rriddw0 0.808 -0.441
## rridew0 0.569 -0.333
## rridsw0 0.470 0.551
## rridtw0 0.475 0.649
## rriduw0 0.570 0.698
## rridvw0 0.584 0.708
## rridww0 0.523 0.594
##
## Factor1 Factor2
## SS loadings 4.034 2.767
## Proportion Var 0.403 0.277
## Cumulative Var 0.403 0.680
##
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 891.63 on 26 degrees of freedom.
## The p-value is 3.2e-171Using lavaan
Lavaan can also be used to do EFA. Inthis case, to extract four factors from the full measure
library(lavaan)## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.##
## Attaching package: 'lavaan'## The following object is masked from 'package:psych':
##
## cor2covmyvars <- c(
"rridaw0","rridbw0","rridcw0","rriddw0","rridew0","rridfw0","rridgw0","rridhw0",
"rridiw0","rridjw0","rridkw0","rridlw0","rridmw0","rridnw0","rridow0","rridpw0",
"rridqw0","rridrw0","rridsw0","rridtw0","rriduw0","rridvw0","rridww0"
)
RFD=AllRFDData[myvars]
#This is legacy code- original data had missing observations. Simulated data does not
CompleteRFD <- na.omit(RFD)
model <- paste0('
# 4 exploratory factors in the same block ("efa4")
efa("efa4")*F1 + efa("efa4")*F2 + efa("efa4")*F3 + efa("efa4")*F4 =~
', paste(myvars, collapse = " + "), '
')
fit_lav <- cfa(
model,
data = CompleteRFD,
std.lv = TRUE, # set factor variances = 1
missing = "fiml",
estimator = "ml",
rotation = "geomin", # <-- rotation name goes here
rotation.args = list(
geomin.epsilon = 0.01, # typical epsilon
orthogonal = FALSE # oblique rotation (factors can correlate)
)
)
summary(fit_lav, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 144
## Row rank of the constraints matrix 12
##
## Rotation method GEOMIN OBLIQUE
## Geomin epsilon 0.01
## Rotation algorithm (rstarts) GPA (30)
## Standardized metric TRUE
## Row weights None
##
## Number of observations 2422
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 2344.266
## Degrees of freedom 167
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 40456.892
## Degrees of freedom 253
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.946
## Tucker-Lewis Index (TLI) 0.918
##
## Robust Comparative Fit Index (CFI) 0.946
## Robust Tucker-Lewis Index (TLI) 0.918
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -54129.245
## Loglikelihood unrestricted model (H1) -52957.112
##
## Akaike (AIC) 108522.490
## Bayesian (BIC) 109287.080
## Sample-size adjusted Bayesian (SABIC) 108867.686
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.073
## 90 Percent confidence interval - lower 0.071
## 90 Percent confidence interval - upper 0.076
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 0.000
##
## Robust RMSEA 0.076
## 90 Percent confidence interval - lower 0.073
## 90 Percent confidence interval - upper 0.079
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 0.013
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.034
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~ efa4
## rridaw0 0.708 0.020 36.007 0.000 0.708 0.766
## rridbw0 0.625 0.022 28.931 0.000 0.625 0.681
## rridcw0 0.883 0.017 52.488 0.000 0.883 0.952
## rriddw0 0.858 0.018 47.849 0.000 0.858 0.893
## rridew0 0.479 0.023 20.560 0.000 0.479 0.543
## rridfw0 0.037 0.015 2.548 0.011 0.037 0.036
## rridgw0 -0.021 0.013 -1.652 0.098 -0.021 -0.021
## rridhw0 0.015 0.014 1.076 0.282 0.015 0.015
## rridiw0 0.145 0.026 5.570 0.000 0.145 0.140
## rridjw0 -0.005 0.013 -0.414 0.679 -0.005 -0.005
## rridkw0 0.086 0.021 4.036 0.000 0.086 0.084
## rridlw0 -0.023 0.030 -0.777 0.437 -0.023 -0.022
## rridmw0 -0.061 0.019 -3.214 0.001 -0.061 -0.091
## rridnw0 -0.100 0.030 -3.318 0.001 -0.100 -0.103
## rridow0 0.122 0.025 4.920 0.000 0.122 0.118
## rridpw0 -0.040 0.021 -1.901 0.057 -0.040 -0.035
## rridqw0 0.007 0.014 0.527 0.598 0.007 0.007
## rridrw0 0.270 0.021 12.955 0.000 0.270 0.288
## rridsw0 0.001 0.015 0.054 0.957 0.001 0.001
## rridtw0 0.002 0.011 0.218 0.827 0.002 0.004
## rriduw0 0.003 0.008 0.378 0.706 0.003 0.005
## rridvw0 -0.005 0.008 -0.575 0.566 -0.005 -0.007
## rridww0 0.050 0.015 3.257 0.001 0.050 0.064
## F2 =~ efa4
## rridaw0 0.025 0.012 2.187 0.029 0.025 0.027
## rridbw0 0.012 0.014 0.866 0.386 0.012 0.014
## rridcw0 -0.002 0.008 -0.184 0.854 -0.002 -0.002
## rriddw0 0.008 0.009 0.905 0.365 0.008 0.009
## rridew0 0.058 0.017 3.491 0.000 0.058 0.066
## rridfw0 0.917 0.018 50.076 0.000 0.917 0.884
## rridgw0 0.899 0.017 52.080 0.000 0.899 0.917
## rridhw0 0.837 0.018 46.142 0.000 0.837 0.824
## rridiw0 0.491 0.021 23.270 0.000 0.491 0.474
## rridjw0 0.893 0.017 51.715 0.000 0.893 0.894
## rridkw0 0.044 0.013 3.345 0.001 0.044 0.043
## rridlw0 0.013 0.024 0.541 0.588 0.013 0.012
## rridmw0 0.140 0.016 8.873 0.000 0.140 0.208
## rridnw0 0.065 0.024 2.709 0.007 0.065 0.067
## rridow0 0.014 0.014 1.009 0.313 0.014 0.013
## rridpw0 0.080 0.020 3.986 0.000 0.080 0.072
## rridqw0 -0.018 0.012 -1.565 0.118 -0.018 -0.018
## rridrw0 -0.056 0.012 -4.539 0.000 -0.056 -0.060
## rridsw0 -0.005 0.012 -0.418 0.676 -0.005 -0.007
## rridtw0 -0.008 0.009 -0.959 0.337 -0.008 -0.013
## rriduw0 -0.002 0.007 -0.259 0.796 -0.002 -0.003
## rridvw0 0.017 0.007 2.489 0.013 0.017 0.026
## rridww0 -0.003 0.011 -0.295 0.768 -0.003 -0.004
## F3 =~ efa4
## rridaw0 0.069 0.018 3.784 0.000 0.069 0.074
## rridbw0 -0.016 0.019 -0.870 0.384 -0.016 -0.018
## rridcw0 -0.012 0.011 -1.015 0.310 -0.012 -0.012
## rriddw0 0.026 0.013 1.947 0.051 0.026 0.027
## rridew0 0.115 0.024 4.725 0.000 0.115 0.130
## rridfw0 -0.029 0.014 -2.064 0.039 -0.029 -0.028
## rridgw0 -0.026 0.013 -1.955 0.051 -0.026 -0.027
## rridhw0 0.070 0.018 3.833 0.000 0.070 0.069
## rridiw0 0.140 0.028 5.005 0.000 0.140 0.135
## rridjw0 -0.010 0.013 -0.802 0.423 -0.010 -0.010
## rridkw0 0.795 0.024 33.162 0.000 0.795 0.772
## rridlw0 0.416 0.033 12.502 0.000 0.416 0.402
## rridmw0 0.107 0.021 5.033 0.000 0.107 0.159
## rridnw0 0.246 0.033 7.406 0.000 0.246 0.252
## rridow0 0.680 0.027 25.663 0.000 0.680 0.659
## rridpw0 0.772 0.028 27.684 0.000 0.772 0.690
## rridqw0 0.898 0.021 43.150 0.000 0.898 0.887
## rridrw0 0.606 0.022 27.269 0.000 0.606 0.645
## rridsw0 0.028 0.017 1.678 0.093 0.028 0.037
## rridtw0 -0.032 0.013 -2.512 0.012 -0.032 -0.048
## rriduw0 0.004 0.009 0.454 0.650 0.004 0.006
## rridvw0 0.011 0.008 1.295 0.195 0.011 0.017
## rridww0 -0.023 0.015 -1.590 0.112 -0.023 -0.030
## F4 =~ efa4
## rridaw0 0.018 0.010 1.813 0.070 0.018 0.019
## rridbw0 0.108 0.015 7.121 0.000 0.108 0.118
## rridcw0 -0.003 0.007 -0.400 0.689 -0.003 -0.003
## rriddw0 -0.001 0.008 -0.130 0.896 -0.001 -0.001
## rridew0 -0.037 0.014 -2.723 0.006 -0.037 -0.042
## rridfw0 -0.025 0.010 -2.500 0.012 -0.025 -0.024
## rridgw0 -0.004 0.009 -0.488 0.626 -0.004 -0.004
## rridhw0 -0.018 0.010 -1.860 0.063 -0.018 -0.018
## rridiw0 0.130 0.018 7.356 0.000 0.130 0.125
## rridjw0 0.020 0.009 2.233 0.026 0.020 0.020
## rridkw0 -0.059 0.012 -4.849 0.000 -0.059 -0.057
## rridlw0 -0.035 0.021 -1.656 0.098 -0.035 -0.033
## rridmw0 0.202 0.014 14.541 0.000 0.202 0.301
## rridnw0 0.142 0.021 6.672 0.000 0.142 0.145
## rridow0 0.082 0.015 5.424 0.000 0.082 0.079
## rridpw0 0.078 0.018 4.438 0.000 0.078 0.070
## rridqw0 -0.006 0.010 -0.652 0.514 -0.006 -0.006
## rridrw0 -0.022 0.008 -2.610 0.009 -0.022 -0.023
## rridsw0 0.540 0.014 38.593 0.000 0.540 0.720
## rridtw0 0.535 0.012 46.053 0.000 0.535 0.818
## rriduw0 0.575 0.010 55.294 0.000 0.575 0.901
## rridvw0 0.596 0.011 56.462 0.000 0.596 0.906
## rridww0 0.619 0.014 44.002 0.000 0.619 0.783
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.399 0.018 22.319 0.000 0.399 0.399
## F3 0.709 0.013 53.564 0.000 0.709 0.709
## F4 0.174 0.020 8.609 0.000 0.174 0.174
## F2 ~~
## F3 0.479 0.017 28.321 0.000 0.479 0.479
## F4 0.325 0.019 17.306 0.000 0.325 0.325
## F3 ~~
## F4 0.178 0.020 8.791 0.000 0.178 0.178
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw0 3.023 0.019 160.998 0.000 3.023 3.271
## .rridbw0 2.935 0.019 157.441 0.000 2.935 3.199
## .rridcw0 3.067 0.019 162.812 0.000 3.067 3.308
## .rriddw0 2.995 0.020 153.421 0.000 2.995 3.117
## .rridew0 3.122 0.018 174.332 0.000 3.122 3.542
## .rridfw0 1.904 0.021 90.315 0.000 1.904 1.835
## .rridgw0 1.761 0.020 88.422 0.000 1.761 1.797
## .rridhw0 1.878 0.021 90.961 0.000 1.878 1.848
## .rridiw0 1.935 0.021 91.941 0.000 1.935 1.868
## .rridjw0 1.801 0.020 88.784 0.000 1.801 1.804
## .rridkw0 2.881 0.021 137.706 0.000 2.881 2.798
## .rridlw0 2.692 0.021 128.059 0.000 2.692 2.602
## .rridmw0 1.389 0.014 101.651 0.000 1.389 2.065
## .rridnw0 1.839 0.020 92.487 0.000 1.839 1.879
## .rridow0 2.566 0.021 122.358 0.000 2.566 2.486
## .rridpw0 2.187 0.023 96.154 0.000 2.187 1.954
## .rridqw0 2.767 0.021 134.527 0.000 2.767 2.734
## .rridrw0 3.098 0.019 162.275 0.000 3.098 3.297
## .rridsw0 1.441 0.015 94.437 0.000 1.441 1.919
## .rridtw0 1.320 0.013 99.373 0.000 1.320 2.019
## .rriduw0 1.296 0.013 99.917 0.000 1.296 2.030
## .rridvw0 1.310 0.013 98.102 0.000 1.310 1.993
## .rridww0 1.469 0.016 91.442 0.000 1.469 1.858
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw0 0.258 0.009 30.128 0.000 0.258 0.302
## .rridbw0 0.424 0.013 32.506 0.000 0.424 0.503
## .rridcw0 0.097 0.006 16.854 0.000 0.097 0.113
## .rriddw0 0.148 0.006 23.283 0.000 0.148 0.161
## .rridew0 0.432 0.013 33.516 0.000 0.432 0.557
## .rridfw0 0.247 0.009 26.184 0.000 0.247 0.230
## .rridgw0 0.191 0.008 23.653 0.000 0.191 0.199
## .rridhw0 0.269 0.010 28.098 0.000 0.269 0.261
## .rridiw0 0.569 0.017 33.553 0.000 0.569 0.531
## .rridjw0 0.199 0.008 24.248 0.000 0.199 0.200
## .rridkw0 0.301 0.011 26.797 0.000 0.301 0.284
## .rridlw0 0.909 0.027 34.003 0.000 0.909 0.850
## .rridmw0 0.357 0.010 34.212 0.000 0.357 0.788
## .rridnw0 0.874 0.025 34.375 0.000 0.874 0.913
## .rridow0 0.429 0.014 30.457 0.000 0.429 0.403
## .rridpw0 0.604 0.020 30.938 0.000 0.604 0.482
## .rridqw0 0.226 0.011 21.107 0.000 0.226 0.221
## .rridrw0 0.257 0.009 29.042 0.000 0.257 0.292
## .rridsw0 0.267 0.009 31.418 0.000 0.267 0.474
## .rridtw0 0.148 0.005 28.844 0.000 0.148 0.347
## .rriduw0 0.076 0.003 23.728 0.000 0.076 0.186
## .rridvw0 0.069 0.003 21.173 0.000 0.069 0.159
## .rridww0 0.236 0.008 30.423 0.000 0.236 0.378
## F1 1.000 1.000 1.000
## F2 1.000 1.000 1.000
## F3 1.000 1.000 1.000
## F4 1.000 1.000 1.000lavInspect(fit_lav, "cor.lv") # factor correlations## F1 F2 F3 F4
## F1 1.000
## F2 0.399 1.000
## F3 0.709 0.479 1.000
## F4 0.174 0.325 0.178 1.000library(dplyr)
library(tidyr)
# Extract standardized factor loadings
loadings_df <- standardizedSolution(fit_lav) %>%
dplyr::filter(op == "=~") %>%
dplyr::select(lhs, rhs, est.std) %>% # pick columns
dplyr::rename( # rename after selecting
Factor = lhs,
Variable = rhs,
Loading = est.std
)
# Sort within each factor by absolute value of loading
loadings_sorted <- loadings_df %>%
arrange(Factor, desc(abs(Loading)))
# Make wide table: variables × factors
loadings_wide <- loadings_sorted %>%
pivot_wider(names_from = Factor, values_from = Loading)
library(knitr)
loadings_wide %>%
mutate(across(-Variable, ~ round(.x, 3))) %>%
kable()| Variable | F1 | F2 | F3 | F4 |
|---|---|---|---|---|
| rridcw0 | 0.952 | -0.002 | -0.012 | -0.003 |
| rriddw0 | 0.893 | 0.009 | 0.027 | -0.001 |
| rridaw0 | 0.766 | 0.027 | 0.074 | 0.019 |
| rridbw0 | 0.681 | 0.014 | -0.018 | 0.118 |
| rridew0 | 0.543 | 0.066 | 0.130 | -0.042 |
| rridrw0 | 0.288 | -0.060 | 0.645 | -0.023 |
| rridiw0 | 0.140 | 0.474 | 0.135 | 0.125 |
| rridow0 | 0.118 | 0.013 | 0.659 | 0.079 |
| rridnw0 | -0.103 | 0.067 | 0.252 | 0.145 |
| rridmw0 | -0.091 | 0.208 | 0.159 | 0.301 |
| rridkw0 | 0.084 | 0.043 | 0.772 | -0.057 |
| rridww0 | 0.064 | -0.004 | -0.030 | 0.783 |
| rridfw0 | 0.036 | 0.884 | -0.028 | -0.024 |
| rridpw0 | -0.035 | 0.072 | 0.690 | 0.070 |
| rridlw0 | -0.022 | 0.012 | 0.402 | -0.033 |
| rridgw0 | -0.021 | 0.917 | -0.027 | -0.004 |
| rridhw0 | 0.015 | 0.824 | 0.069 | -0.018 |
| rridqw0 | 0.007 | -0.018 | 0.887 | -0.006 |
| rridvw0 | -0.007 | 0.026 | 0.017 | 0.906 |
| rridjw0 | -0.005 | 0.894 | -0.010 | 0.020 |
| rriduw0 | 0.005 | -0.003 | 0.006 | 0.901 |
| rridtw0 | 0.004 | -0.013 | -0.048 | 0.818 |
| rridsw0 | 0.001 | -0.007 | 0.037 | 0.720 |
write.table(loadings_wide, "clipboard", sep="\t", row.names = FALSE)Vector Plot of Unrotated Loadings (Figure16.1)
# - Setup -
library(stats)
library(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha# 1) Read and subset
AllRFDData <- read.table("RFDSim.txt", sep = "\t", header = TRUE)
myvars <- c("rridaw0","rridbw0","rridcw0","rriddw0","rridew0",
"rridsw0","rridtw0","rriduw0","rridvw0","rridww0")
RFD <- AllRFDData[myvars]
# Optional: listwise deletion to match n.obs with the correlation matrix
RFD_complete <- na.omit(RFD)
# 2) Correlation matrix and ML EFA (unrotated)
R <- cor(RFD_complete, use = "pairwise.complete.obs")
fa_result <- factanal(covmat = R, n.obs = nrow(RFD_complete),
factors = 2, rotation = "none") # ML + unrotated
print(fa_result, digits = 3)##
## Call:
## factanal(factors = 2, covmat = R, n.obs = nrow(RFD_complete), rotation = "none")
##
## Uniquenesses:
## rridaw0 rridbw0 rridcw0 rriddw0 rridew0 rridsw0 rridtw0 rriduw0 rridvw0 rridww0
## 0.308 0.506 0.119 0.152 0.566 0.475 0.353 0.188 0.158 0.373
##
## Loadings:
## Factor1 Factor2
## rridaw0 0.744 -0.372
## rridbw0 0.664 -0.229
## rridcw0 0.822 -0.453
## rriddw0 0.808 -0.441
## rridew0 0.569 -0.333
## rridsw0 0.470 0.551
## rridtw0 0.475 0.649
## rriduw0 0.570 0.698
## rridvw0 0.584 0.708
## rridww0 0.523 0.594
##
## Factor1 Factor2
## SS loadings 4.034 2.767
## Proportion Var 0.403 0.277
## Cumulative Var 0.403 0.680
##
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 891.63 on 26 degrees of freedom.
## The p-value is 3.2e-171# 3) Extract loadings (variables x 2)
L <- as.matrix(fa_result$loadings)[, 1:2]
colnames(L) <- c("F1","F2")
# 4) Square vector plot WITHOUT labels (ggplot2)
dfL <- data.frame(var = rownames(L), F1 = L[,1], F2 = L[,2])
ggplot(dfL) +
# arrows from origin
geom_segment(aes(x = 0, y = 0, xend = F1, yend = F2),
arrow = arrow(length = unit(0.18, "cm"))) +
# reference axes
geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
geom_vline(xintercept = 0, linetype = 2, color = "grey60") +
# (optional) unit circle guide
annotate("path",
x = cos(seq(0, 2*pi, length.out = 400)),
y = sin(seq(0, 2*pi, length.out = 400)),
linetype = 3, color = "grey70") +
coord_equal(xlim = c(-1,1), ylim = c(-1,1)) + # square aspect & bounds
labs(x = "Factor 1 (unrotated ML)", y = "Factor 2 (unrotated ML)",
title = "Unrotated Two-Factor Vector Plot (No Labels)") +
theme_minimal(base_size = 12)
Example of Using Direction Cosines to Rotate Factor Loadings
library(dplyr)
library(ggplot2)
library(grid)
# Loadings were already calculated in previous chunk.
# 1) Making the translation matarix T of direction cosines to rotate factor loadings by 56° (clockwise)
angle_deg <- 56
theta <- angle_deg * pi / 180 #Making it radians
T <- matrix(c(cos(theta), -sin(theta),
sin(theta), cos(theta)),
nrow = 2, byrow = TRUE,
dimnames = list(c("F1","F2"), c("F1","F2")))
cat("\nDirection cosine rotation matrix for", angle_deg, "degrees (clockwise):\n")##
## Direction cosine rotation matrix for 56 degrees (clockwise):print(round(T, 6))## F1 F2
## F1 0.559193 -0.829038
## F2 0.829038 0.559193# 2) Rotated loadings "L_rot"
L_rot <- L %*% T
colnames(L_rot) <- c("F1_rot","F2_rot")
# 3) Build numeric table: original vs rotated (rounded)
library(dplyr)
tab <- data.frame(
Variable = rownames(L),
F1 = L[,1],
F2 = L[,2],
F1_rot = L_rot[,1],
F2_rot = L_rot[,2],
row.names = NULL
) %>%
mutate(across(-Variable, ~ round(.x, 3))) %>%
mutate(max_abs_rot = pmax(abs(F1_rot), abs(F2_rot))) %>%
arrange(desc(max_abs_rot)) %>%
dplyr::select(-max_abs_rot) # now explicitly using dplyr's select()
cat("\nOriginal and Rotated (56°) Loadings:\n")##
## Original and Rotated (56°) Loadings:print(tab, row.names = FALSE)## Variable F1 F2 F1_rot F2_rot
## rridcw0 0.822 -0.453 0.084 -0.935
## rriddw0 0.808 -0.441 0.087 -0.917
## rridvw0 0.584 0.708 0.914 -0.089
## rriduw0 0.570 0.698 0.897 -0.083
## rridaw0 0.744 -0.372 0.108 -0.825
## rridtw0 0.475 0.649 0.804 -0.031
## rridww0 0.523 0.594 0.785 -0.102
## rridsw0 0.470 0.551 0.720 -0.082
## rridbw0 0.664 -0.229 0.182 -0.679
## rridew0 0.569 -0.333 0.042 -0.658# 4) Vector plot: original (gray) vs rotated (red), square, no labels
df_orig <- data.frame(var = rownames(L), F1 = L[,1], F2 = L[,2], type = "Original")
df_rot <- data.frame(var = rownames(L), F1 = L_rot[,1], F2 = L_rot[,2], type = "Rotated (56°)")
df_all <- bind_rows(df_orig, df_rot)
unique(df_all$type)## [1] "Original" "Rotated (56°)"ggplot(df_all) +
geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
geom_vline(xintercept = 0, linetype = 2, color = "grey60") +
annotate("path",
x = cos(seq(0, 2*pi, length.out = 400)),
y = sin(seq(0, 2*pi, length.out = 400)),
linetype = 3, color = "grey75") +
geom_segment(aes(x = 0, y = 0, xend = F1, yend = F2, color = type),
arrow = arrow(length = unit(0.18, "cm")), linewidth = 0.7) +
scale_color_manual(values = c("Original" = "grey40", "Rotated (56°)" = "red")) +
coord_equal(xlim = c(-1, 1), ylim = c(-1, 1)) +
labs(x = "Factor 1", y = "Factor 2",
title = "Unrotated Loadings vs. 56° Rotation",
subtitle = "Rotation via 2×2 direction cosine matrix") +
theme_minimal(base_size = 12) +
theme(legend.title = element_blank())
Factor Analysis of an Oblique (“Correlated) Factors
The example uses the Enhancement and Social Reasons for drinking subscales (Figure 16.4)
# - Setup -
library(stats)
library(ggplot2)
# 1) Read and subset
AllRFDData <- read.table("RFDSim.txt", sep = "\t", header = TRUE)
myvars <- c("rridcw0","rriddw0","rridaw0","rridbw0","rridew0",
"rridqw0","rridkw0","rridpw0","rridow0","rridrw0")
RFD <- AllRFDData[myvars]
# Optional: Often researchers use listwise deletion to match n.obs with the correlation matrix
# In this example this isn't necessary as there is no missing data. As mentioned in the book,
# with large samples one may simply do pairwise deletion, but with moderate or small samples
# pairwise correlation matrices may be misleading or result in rank deficiencies.
RFD_complete <- na.omit(RFD)
# 2) Correlation matrix and ML EFA (unrotated)
R <- cor(RFD_complete, use = "pairwise.complete.obs")
fa_result <- factanal(covmat = R, n.obs = nrow(RFD_complete),
factors = 2, rotation = "none") # ML + unrotated
print(fa_result, digits = 3.,cutoff=0)##
## Call:
## factanal(factors = 2, covmat = R, n.obs = nrow(RFD_complete), rotation = "none")
##
## Uniquenesses:
## rridcw0 rriddw0 rridaw0 rridbw0 rridew0 rridqw0 rridkw0 rridpw0 rridow0 rridrw0
## 0.114 0.158 0.304 0.520 0.560 0.213 0.290 0.494 0.414 0.295
##
## Loadings:
## Factor1 Factor2
## rridcw0 0.905 -0.259
## rriddw0 0.889 -0.228
## rridaw0 0.819 -0.158
## rridbw0 0.668 -0.182
## rridew0 0.658 -0.089
## rridqw0 0.772 0.436
## rridkw0 0.762 0.360
## rridpw0 0.615 0.359
## rridow0 0.709 0.288
## rridrw0 0.807 0.233
##
## Factor1 Factor2
## SS loadings 5.867 0.771
## Proportion Var 0.587 0.077
## Cumulative Var 0.587 0.664
##
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 440.52 on 26 degrees of freedom.
## The p-value is 6.33e-77# 3) Extract loadings (variables x 2)
L <- as.matrix(fa_result$loadings)[, 1:2]
colnames(L) <- c("F1","F2")
# 4) Square vector plot WITHOUT labels (ggplot2)
dfL <- data.frame(var = rownames(L), F1 = L[,1], F2 = L[,2])
ggplot(dfL) +
# arrows from origin
geom_segment(aes(x = 0, y = 0, xend = F1, yend = F2),
arrow = arrow(length = unit(0.18, "cm"))) +
# reference axes
geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
geom_vline(xintercept = 0, linetype = 2, color = "grey60") +
# (optional) unit circle guide
annotate("path",
x = cos(seq(0, 2*pi, length.out = 400)),
y = sin(seq(0, 2*pi, length.out = 400)),
linetype = 3, color = "grey70") +
coord_equal(xlim = c(-1,1), ylim = c(-1,1)) + # square aspect & bounds
labs(x = "Factor 1 (unrotated ML)", y = "Factor 2 (unrotated ML)",
title = "Enhancement and Social Reasons (Figure 16.4)") +
theme_minimal(base_size = 12)
Oblique Rotation (Figure 16.5)
library(psych)
library(GPArotation)
# Again, Enhancement and Social Reasons Items
myvars <- c("rridcw0","rriddw0","rridaw0","rridbw0","rridew0",
"rridqw0","rridkw0","rridpw0","rridow0","rridrw0")
RFD_sub <- AllRFDData[myvars]
# 1. Unrotated maximum likelihood factor analysis (iterated)
fa_ml <- fa(RFD_sub, nfactors = 2, fm = "ml", rotate = "none")
# 2. Varimax rotation
fa_varimax <- fa(RFD_sub, nfactors = 2, fm = "ml", rotate = "varimax")
L_varimax <- fa_varimax$loadings
# 3. Promax rotation (m = 3)
fa_promax <- fa(RFD_sub, nfactors = 2, fm = "ml", rotate = "promax")
L_promax <- fa_promax$loadings
# 4. Target rotation to the Promax target (m = 3)
# Build target matrix U
A <- as.matrix(L_varimax)
U <- sign(A) * abs(A)^3 # promax target
# Target rotation
tgt <- GPArotation::targetQ(A, Target = U, normalize = FALSE)
L_target <- tgt$loadings
# 5. Combine into one table
rotations_df <- data.frame(
Variable = rownames(A),
Varimax_F1 = round(L_varimax[,1], 3),
Varimax_F2 = round(L_varimax[,2], 3),
Promax_F1 = round(L_promax[,1], 3),
Promax_F2 = round(L_promax[,2], 3),
Target_F1 = round(L_target[,1], 3),
Target_F2 = round(L_target[,2], 3)
)
require(knitr)
knitr::kable(rotations_df, caption = "Rotations: Varimax, Promax, Target Columns (F1–F2)")| Variable | Varimax_F1 | Varimax_F2 | Promax_F1 | Promax_F2 | Target_F1 | Target_F2 | |
|---|---|---|---|---|---|---|---|
| rridcw0 | rridcw0 | 0.867 | 0.367 | 0.960 | -0.026 | 0.870 | 0.122 |
| rriddw0 | rriddw0 | 0.834 | 0.382 | 0.908 | 0.013 | 0.828 | 0.149 |
| rridaw0 | rridaw0 | 0.736 | 0.392 | 0.769 | 0.085 | 0.713 | 0.194 |
| rridbw0 | rridbw0 | 0.634 | 0.278 | 0.697 | -0.006 | 0.633 | 0.100 |
| rridew0 | rridew0 | 0.567 | 0.345 | 0.567 | 0.124 | 0.535 | 0.199 |
| rridqw0 | rridqw0 | 0.326 | 0.825 | -0.043 | 0.918 | 0.101 | 0.828 |
| rridkw0 | rridkw0 | 0.366 | 0.759 | 0.050 | 0.805 | 0.169 | 0.739 |
| rridpw0 | rridpw0 | 0.253 | 0.665 | -0.049 | 0.747 | 0.070 | 0.672 |
| rridow0 | rridow0 | 0.370 | 0.670 | 0.108 | 0.682 | 0.202 | 0.636 |
| rridrw0 | rridrw0 | 0.481 | 0.688 | 0.248 | 0.640 | 0.324 | 0.619 |
Oblique Rotation Vector Plot
This plot will not resemble Figures 15.4 and 15.5 because the data are simulated. The regression lines will therefore look straight and the likert-format of the questions is not evident. Reversing the order of the questions (from rridew0 to rridaw0) will produce Figure 15.5
# - Packages -
# install.packages(c("psych","GPArotation","ggplot2"))
library(psych)
library(GPArotation)
library(ggplot2)
library(grid)
# Again, looking at Enhancement and Social Reasons
# AllRFDData <- read.table("RFDSim.txt", sep = "\t", header = TRUE)
myvars <- c("rridcw0","rriddw0","rridaw0","rridbw0","rridew0",
"rridqw0","rridkw0","rridpw0","rridow0","rridrw0")
X <- na.omit(AllRFDData[, myvars])
# - Settings -
nf <- 2 # number of factors for the plot; code works for nf >= 2, plot shows first two
# - Unrotated ML EFA on the correlation matrix -
R <- cor(X, use = "pairwise.complete.obs")
fa0 <- factanal(covmat = R, n.obs = nrow(X), factors = nf, rotation = "none")
# Unrotated loadings (variables x nf)
A <- as.matrix(fa0$loadings)[, 1:nf, drop = FALSE]
rownames(A) <- colnames(R)
# - Oblique GEOMIN rotation (GeominQ) -
# Note: normalize = TRUE mimics common practice (Kaiser normalization before rotation),
# matching what psych::fa(..., rotate="geominQ") would typically do.
geo <- GPArotation::geominQ(A, normalize = TRUE) # oblique rotation
# Extract rotated results
L_rot <- as.matrix(geo$loadings) # Pattern loadings (Λ)
Phi <- geo$Phi # Factor correlation matrix (Φ)
Tmat <- geo$Th # Transformation (rotation) matrix (T)
# Structure matrix (optional): Σ = Λ Φ
Struct <- L_rot %*% Phi
# - Print tables (rounded) -
cat("\nRotated pattern loadings (Geomin, oblique):\n")##
## Rotated pattern loadings (Geomin, oblique):print(round(L_rot, 3))## Factor1 Factor2
## rridcw0 0.969 -0.038
## rriddw0 0.916 0.002
## rridaw0 0.776 0.076
## rridbw0 0.704 -0.015
## rridew0 0.573 0.117
## rridqw0 -0.039 0.916
## rridkw0 0.055 0.801
## rridpw0 -0.046 0.745
## rridow0 0.112 0.679
## rridrw0 0.253 0.635cat("\nFactor correlations (Phi):\n")##
## Factor correlations (Phi):print(round(Phi, 3))## Factor1 Factor2
## Factor1 1.000 0.741
## Factor2 0.741 1.000cat("\nTransformation (rotation) matrix T:\n")##
## Transformation (rotation) matrix T:print(round(Tmat, 3))## [,1] [,2]
## [1,] 0.968 0.885
## [2,] -0.249 0.466# - Plot of rotated solution (first two factors) -
if (ncol(L_rot) < 2) stop("Need at least 2 factors to make the 2D vector plot.")
dfL <- data.frame(var = rownames(L_rot), F1 = L_rot[,1], F2 = L_rot[,2])
library(ggplot2)
library(grid) # for unit(), arrow()
# - Angle (from factor correlation or fixed) -
theta_deg <- 42.18
theta <- theta_deg * pi/180
# - Transform oblique coordinates (F1,F2) -> Cartesian (x,y) -
df_plot <- dfL
df_plot$x <- with(df_plot, F1 + F2 * cos(theta))
df_plot$y <- with(df_plot, F2 * sin(theta))
# - Oblique axes (length a bit beyond the plotting region) -
ax_len <- 1.1
axis_F1 <- data.frame(x1 = -ax_len, y1 = 0, x2 = ax_len, y2 = 0)
axis_F2 <- data.frame(
x1 = -ax_len * cos(theta), y1 = -ax_len * sin(theta),
x2 = ax_len * cos(theta), y2 = ax_len * sin(theta)
)
# - "Unit circle" in oblique coords becomes an ellipse in (x,y) -
t <- seq(0, 2*pi, length.out = 400)
ell <- data.frame(
x = cos(t) + sin(t) * cos(theta),
y = sin(t) * sin(theta)
)
ggplot(df_plot) +
# ellipse = image of the unit circle under the oblique mapping
geom_path(data = ell, aes(x, y), linetype = 3, color = "grey70") +
# oblique axes (both directions)
geom_segment(data = axis_F1, aes(x = x1, y = y1, xend = x2, yend = y2),
linewidth = 0.7,
arrow = arrow(length = unit(0.16, "cm"), ends = "both")) +
geom_segment(data = axis_F2, aes(x = x1, y = y1, xend = x2, yend = y2),
linewidth = 0.7,
arrow = arrow(length = unit(0.16, "cm"), ends = "both")) +
# your vectors, now plotted IN the oblique frame
geom_segment(aes(x = 0, y = 0, xend = x, yend = y),
arrow = arrow(length = unit(0.18, "cm")), linewidth = 0.7) +
geom_text(aes(x, y, label = var), vjust = -0.6, size = 3) +
coord_equal(xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1)) +
labs(
x = "F1 axis",
y = "F2 axis",
title = sprintf("Oblique Factor Space (axes separated by %.2f°)", theta_deg),
subtitle = "Enhancement & Social vectors in oblique space\nEllipse is the image of the unit circle"
) +
theme_minimal(base_size = 12)
Tucker’s Factor Congruence Measure (Table 16.5)
# - Packages - (already installed but if you're running this module by itself)
# install.packages(c("psych","GPArotation"))
library(psych)
library(GPArotation)
# - Data & variables -
# AllRFDData <- read.table("RFDSim.txt", sep="\t", header=TRUE) # if not already loaded
myvars <- c("rridcw0","rriddw0","rridaw0","rridbw0","rridew0",
"rridqw0","rridkw0","rridpw0","rridow0","rridrw0")
X <- na.omit(AllRFDData[ , myvars])
# - Fit 2-factor ML EFA with three rotations (pattern loadings) -
nf <- 2
fa_varimax <- fa(X, nfactors = nf, fm = "ml", rotate = "varimax") # orthogonal
fa_promax <- fa(X, nfactors = nf, fm = "ml", rotate = "promax") # oblique
fa_geomin <- fa(X, nfactors = nf, fm = "ml", rotate = "geominQ") # oblique
# - Extract loadings
L_vmx <- as.matrix(fa_varimax$loadings)
L_pmx <- as.matrix(fa_promax$loadings)
L_geo <- as.matrix(fa_geomin$loadings)
# - Rather than go throught the work of figuring out which factor goes with which
# - it makes sense to try to match up the most representative factors
# - Helper function to align columns/signs of B to A (2-factor case) using congruence -
align_to <- function(A, B) {
# Try both permutations
perms <- list(c(1,2), c(2,1))
best <- NULL; best_score <- -Inf
for (p in perms) {
Bp <- B[, p, drop=FALSE]
# Flip signs to maximize diagonal congruence
C <- factor.congruence(A, Bp) # 2x2 matrix of φ
sgn <- sign(diag(C)); sgn[sgn==0] <- 1
Bps <- sweep(Bp, 2, sgn, `*`)
C2 <- factor.congruence(A, Bps)
score <- sum(abs(diag(C2)))
if (score > best_score) { best <- Bps; best_score <- score }
}
best
}
# Align Promax and Geomin to Varimax
L_pmx_a <- align_to(L_vmx, L_pmx)
L_geo_a <- align_to(L_vmx, L_geo)
# - Compute Tucker's congruence (φ) matrices -
phi_vmx_pmx <- factor.congruence(L_vmx, L_pmx_a)
phi_vmx_geo <- factor.congruence(L_vmx, L_geo_a)
# Also Promax vs Geomin (align geomin to promax as well)
L_geo_to_pmx <- align_to(L_pmx, L_geo)
phi_pmx_geo <- factor.congruence(L_pmx, L_geo_to_pmx)
# - Print results -
cat("\nTucker's congruence (φ) Varimax vs Promax (aligned to Varimax):\n")##
## Tucker's congruence (φ) Varimax vs Promax (aligned to Varimax):print(round(phi_vmx_pmx, 3))## ML1 ML2
## ML1 0.93 0.46
## ML2 0.50 0.91cat("\nTucker's congruence (φ) Varimax vs Geomin (aligned to Varimax):\n")##
## Tucker's congruence (φ) Varimax vs Geomin (aligned to Varimax):print(round(phi_vmx_geo, 3))## ML1 ML2
## ML1 0.93 0.46
## ML2 0.50 0.91cat("\nTucker's congruence (φ) Promax vs Geomin (aligned to Promax):\n")##
## Tucker's congruence (φ) Promax vs Geomin (aligned to Promax):print(round(phi_pmx_geo, 3))## ML1 ML2
## ML1 1.0 0.1
## ML2 0.1 1.0# Optional: diagonals only (per-factor matches)
cat("\nDiagonal φ (best-matching factors):\n")##
## Diagonal φ (best-matching factors):diag_tab <- data.frame(
Pair = c("Vmx–Pmx F1","Vmx–Pmx F2","Vmx–Geo F1","Vmx–Geo F2","Pmx–Geo F1","Pmx–Geo F2"),
Phi = round(c(diag(phi_vmx_pmx), diag(phi_vmx_geo), diag(phi_pmx_geo)), 3)
)
print(diag_tab, row.names = FALSE)## Pair Phi
## Vmx–Pmx F1 0.93
## Vmx–Pmx F2 0.91
## Vmx–Geo F1 0.93
## Vmx–Geo F2 0.91
## Pmx–Geo F1 1.00
## Pmx–Geo F2 1.00# Examine factor correlations for the two oblique rotations
cat("\nFactor correlations (Phi) — Promax:\n"); print(round(fa_promax$Phi, 3))##
## Factor correlations (Phi) — Promax:## ML1 ML2
## ML1 1.000 0.737
## ML2 0.737 1.000cat("\nFactor correlations (Phi) — Geomin:\n"); print(round(fa_geomin$Phi, 3))##
## Factor correlations (Phi) — Geomin:## ML1 ML2
## ML1 1.000 0.738
## ML2 0.738 1.000Bifactor Rotations (Table 16.6) and Congruence Coefficients for Bifactor Rotations (Table 16.7)
# - Packages -
# install.packages(c("GPArotation","psych")) # if needed
library(psych) # for quick data handling / correlations
library(GPArotation) # rotations incl. bigeominQ and CF family
# - Data -
AllRFDData <- read.table("RFDSim.txt", sep = "\t", header = TRUE)
myvars <- c(
"rridaw0","rridbw0","rridcw0","rriddw0","rridew0","rridfw0","rridgw0",
"rridhw0","rridiw0","rridjw0","rridkw0","rridlw0","rridmw0","rridnw0",
"rridow0","rridpw0","rridqw0","rridrw0","rridsw0","rridtw0","rriduw0",
"rridvw0","rridww0"
)
X <- na.omit(AllRFDData[ , myvars]) # drop rows with any missing on these vars
R <- cor(X) # correlation matrix for ML EFA
# - Factor extraction (unrotated) -
# Choose total number of factors (1 general + S specifics). Adjust as needed.
nf <- 4 # example: 1 general + 3 specifics. Change to what your design warrants.
fa0 <- factanal(covmat = R, factors = nf, rotation = "none", n.obs = nrow(X))
L <- loadings(fa0) # unrotated loading matrix (class "loadings"); coerce to matrix if needed
# Optional: Kaiser normalization before rotation often helps GPA convergence
norm <- TRUE
rs <- 200 # random starts to mitigate local minima (esp. for geomin/bigeomin)
# You can't, at least as of this writing, do a bifactor CF-Quartimax rotation in R.
# The code here uses cfQ(..., kappa = 0) which is CF-Quartimax (oblique), not the BI-CF-Quartimax variant.
# If you have the Mplus installed on your computer, you can use MplusAutomation to submit the code and retrieve
# the rotated factor loadings back into R as shown in the RMarkdown chunk below this one.
# - BI-GEOMIN (oblique) -
# GPArotation 2025+ provides bigeominQ/bigeominT
rot_bigeomin <- bigeominQ(L, normalize = norm, randomStarts = rs)
# - CF-QUARTIMAX (oblique, κ = 0) -
# This is CF family with kappa = 0 (aka CF-quartimax). NOTE: this is NOT the bifactor version.
rot_cf_quartimax <- cfQ(L, kappa = 0, normalize = norm, randomStarts = rs)
# - Summaries -
cat("\n=== BI-GEOMIN (oblique) loadings ===\n")##
## === BI-GEOMIN (oblique) loadings ===print(rot_bigeomin, digits = 3, sortLoadings = TRUE)## Oblique rotation method Bi-Geomin converged at lowest minimum.
## Of 200 random starts 100% converged, 86% at the same lowest minimum.
## Random starts converged to 2 different local minima
## Loadings at lowest minimum:
## Factor1 Factor2 Factor3 Factor4
## rridaw0 0.665 0.00695 -0.000344 0.50658
## rridbw0 0.531 0.10884 0.005049 0.45579
## rridcw0 0.694 -0.00989 -0.009654 0.63560
## rriddw0 0.696 -0.01019 -0.007887 0.59419
## rridew0 0.561 -0.05439 0.022612 0.35458
## rridfw0 0.551 -0.02762 0.689889 0.01747
## rridgw0 0.537 -0.00862 0.715820 -0.02105
## rridhw0 0.594 -0.02798 0.626464 -0.00115
## rridiw0 0.571 0.10674 0.341230 0.08205
## rridjw0 0.558 0.01399 0.694980 -0.01083
## rridkw0 0.830 -0.10797 -0.101968 0.01346
## rridlw0 0.376 -0.05903 -0.059871 -0.03677
## rridmw0 0.304 0.27798 0.132329 -0.07006
## rridnw0 0.253 0.12306 0.007599 -0.08219
## rridow0 0.764 0.03114 -0.106493 0.04321
## rridpw0 0.715 0.02055 -0.064711 -0.06126
## rridqw0 0.861 -0.06555 -0.168651 -0.04287
## rridrw0 0.801 -0.06733 -0.161506 0.15740
## rridsw0 0.235 0.68989 -0.017775 0.00124
## rridtw0 0.177 0.79032 -0.008575 0.00791
## rriduw0 0.261 0.86641 -0.010964 0.00606
## rridvw0 0.282 0.87041 0.009434 -0.00278
## rridww0 0.235 0.75451 -0.005403 0.04672
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 7.339 3.333 2.095 1.402
## Proportion Var 0.319 0.145 0.091 0.061
## Cumulative Var 0.319 0.464 0.555 0.616
##
## Phi:
## Factor1 Factor2 Factor3 Factor4
## Factor1 1.00e+00 -1.98e-05 1.38e-05 -2.15e-06
## Factor2 -1.98e-05 1.00e+00 2.15e-01 -4.36e-02
## Factor3 1.38e-05 2.15e-01 1.00e+00 -1.06e-01
## Factor4 -2.15e-06 -4.36e-02 -1.06e-01 1.00e+00cat("\nFactor correlations (Phi), BI-GEOMIN:\n")##
## Factor correlations (Phi), BI-GEOMIN:print(rot_bigeomin$Phi, digits = 3)## Factor1 Factor2 Factor3 Factor4
## Factor1 1.00e+00 1.38e-05 -1.98e-05 2.15e-06
## Factor2 1.38e-05 1.00e+00 2.15e-01 1.06e-01
## Factor3 -1.98e-05 2.15e-01 1.00e+00 4.36e-02
## Factor4 2.15e-06 1.06e-01 4.36e-02 1.00e+00cat("\n=== CF-QUARTIMAX (oblique, kappa=0) loadings ===\n")##
## === CF-QUARTIMAX (oblique, kappa=0) loadings ===print(rot_cf_quartimax, digits = 3, sortLoadings = TRUE)## Oblique rotation method Crawford-Ferguson Quartimax/Quartimin converged at lowest minimum.
## Of 200 random starts 100% converged, 100% at the same lowest minimum.
## Loadings at lowest minimum:
## Factor1 Factor2 Factor3 Factor4
## rridaw0 0.74952 0.02771 0.04820 0.08958
## rridbw0 0.66338 0.13023 0.03186 -0.00359
## rridcw0 0.92272 0.01296 0.02636 0.01112
## rriddw0 0.86886 0.01203 0.03431 0.04795
## rridew0 0.53667 -0.04130 0.08052 0.13901
## rridfw0 0.04621 -0.04421 0.89843 -0.04537
## rridgw0 -0.00820 -0.02627 0.93003 -0.04648
## rridhw0 0.03228 -0.04155 0.83460 0.04955
## rridiw0 0.15639 0.11049 0.47961 0.12136
## rridjw0 0.00906 -0.00135 0.90675 -0.02999
## rridkw0 0.13716 -0.09187 0.03052 0.75215
## rridlw0 0.00732 -0.05224 0.00412 0.39006
## rridmw0 -0.06518 0.28942 0.19994 0.14073
## rridnw0 -0.07615 0.13171 0.05659 0.23734
## rridow0 0.16532 0.05231 0.00170 0.64062
## rridpw0 0.01951 0.03719 0.05613 0.66580
## rridqw0 0.07204 -0.04530 -0.03663 0.86200
## rridrw0 0.32493 -0.04574 -0.06513 0.63465
## rridsw0 0.02296 0.72386 -0.02043 0.02062
## rridtw0 0.02161 0.82705 -0.02665 -0.06433
## rriduw0 0.02932 0.90803 -0.01900 -0.01337
## rridvw0 0.01900 0.91169 0.00931 -0.00419
## rridww0 0.08041 0.79135 -0.01576 -0.04407
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 3.658 3.639 3.589 3.282
## Proportion Var 0.159 0.158 0.156 0.143
## Cumulative Var 0.159 0.317 0.473 0.616
##
## Phi:
## Factor1 Factor2 Factor3 Factor4
## Factor1 1.000 0.139 0.360 0.647
## Factor2 0.139 1.000 0.364 0.237
## Factor3 0.360 0.364 1.000 0.513
## Factor4 0.647 0.237 0.513 1.000cat("\nFactor correlations (Phi), CF-QUARTIMAX:\n")##
## Factor correlations (Phi), CF-QUARTIMAX:print(rot_cf_quartimax$Phi, digits = 3)## Factor1 Factor2 Factor3 Factor4
## Factor1 1.000 -0.647 -0.139 -0.360
## Factor2 -0.647 1.000 0.237 0.513
## Factor3 -0.139 0.237 1.000 0.364
## Factor4 -0.360 0.513 0.364 1.000# Matrices you might want to keep explicitly:
L_bigeomin <- rot_bigeomin$loadings
Phi_bigeomin <- rot_bigeomin$Phi
T_bigeomin <- rot_bigeomin$Th # transformation (rotation) matrix
L_cfquart <- rot_cf_quartimax$loadings
Phi_cfquart <- rot_cf_quartimax$Phi
T_cfquart <- rot_cf_quartimax$Th
# - Quick comparison table of largest loadings by factor (optional) -
# shows, for each variable, its largest absolute loading and which factor that is
largest_loading <- function(Lmat) {
M <- as.matrix(Lmat)
idx <- apply(abs(M), 1, which.max)
val <- M[cbind(seq_len(nrow(M)), idx)]
data.frame(variable = rownames(M),
max_loading = round(val, 3),
max_factor = paste0("F", idx),
row.names = NULL)
}
cat("\n=== Largest loading per variable: BI-GEOMIN ===\n")##
## === Largest loading per variable: BI-GEOMIN ===print(largest_loading(L_bigeomin))## variable max_loading max_factor
## 1 rridaw0 -0.665 F1
## 2 rridbw0 -0.531 F1
## 3 rridcw0 -0.694 F1
## 4 rriddw0 -0.696 F1
## 5 rridew0 -0.561 F1
## 6 rridfw0 -0.690 F2
## 7 rridgw0 -0.716 F2
## 8 rridhw0 -0.626 F2
## 9 rridiw0 -0.571 F1
## 10 rridjw0 -0.695 F2
## 11 rridkw0 -0.830 F1
## 12 rridlw0 -0.376 F1
## 13 rridmw0 -0.304 F1
## 14 rridnw0 -0.253 F1
## 15 rridow0 -0.764 F1
## 16 rridpw0 -0.715 F1
## 17 rridqw0 -0.861 F1
## 18 rridrw0 -0.801 F1
## 19 rridsw0 -0.690 F3
## 20 rridtw0 -0.790 F3
## 21 rriduw0 -0.866 F3
## 22 rridvw0 -0.870 F3
## 23 rridww0 -0.755 F3cat("\n=== Largest loading per variable: CF-QUARTIMAX ===\n")##
## === Largest loading per variable: CF-QUARTIMAX ===print(largest_loading(L_cfquart))## variable max_loading max_factor
## 1 rridaw0 -0.750 F1
## 2 rridbw0 -0.663 F1
## 3 rridcw0 -0.923 F1
## 4 rriddw0 -0.869 F1
## 5 rridew0 -0.537 F1
## 6 rridfw0 0.898 F4
## 7 rridgw0 0.930 F4
## 8 rridhw0 0.835 F4
## 9 rridiw0 0.480 F4
## 10 rridjw0 0.907 F4
## 11 rridkw0 0.752 F2
## 12 rridlw0 0.390 F2
## 13 rridmw0 0.289 F3
## 14 rridnw0 0.237 F2
## 15 rridow0 0.641 F2
## 16 rridpw0 0.666 F2
## 17 rridqw0 0.862 F2
## 18 rridrw0 0.635 F2
## 19 rridsw0 0.724 F3
## 20 rridtw0 0.827 F3
## 21 rriduw0 0.908 F3
## 22 rridvw0 0.912 F3
## 23 rridww0 0.791 F3# - Congruence (Tucker's φ) between the two rotated solutions -
# Ensure we have numeric matrices with same variables in same order
Lb <- as.matrix(L_bigeomin)
Lc <- as.matrix(L_cfquart)
stopifnot(nrow(Lb) == nrow(Lc), ncol(Lb) == ncol(Lc))
stopifnot(identical(rownames(Lb), rownames(Lc)))
# Function to compute Tucker's congruence for two vectors
tucker_phi_vec <- function(x, y) sum(x * y) / sqrt(sum(x^2) * sum(y^2))
# Matrix of φ across all factor pairs (BI-GEOMIN columns vs CF-Quartimax columns)
tucker_matrix <- function(L1, L2) {
p <- ncol(L1)
M <- matrix(NA_real_, nrow = p, ncol = p,
dimnames = list(paste0("BIgeo_F", 1:p), paste0("CFQ_F", 1:p)))
for (i in seq_len(p)) {
for (j in seq_len(p)) M[i, j] <- tucker_phi_vec(L1[, i], L2[, j])
}
M
}
Phi_cross <- tucker_matrix(Lb, Lc)
cat("\n=== Tucker's congruence (φ) matrix: BI-GEOMIN vs CF-Quartimax ===\n")##
## === Tucker's congruence (φ) matrix: BI-GEOMIN vs CF-Quartimax ===print(round(Phi_cross, 3))## CFQ_F1 CFQ_F2 CFQ_F3 CFQ_F4
## BIgeo_F1 0.652 -0.716 -0.211 -0.484
## BIgeo_F2 -0.001 0.184 0.002 -0.980
## BIgeo_F3 0.036 0.059 -0.999 -0.005
## BIgeo_F4 -0.988 0.076 0.042 0.044# - Best one-to-one factor matching (maximize |φ|) via Hungarian algorithm -
if (!requireNamespace("clue", quietly = TRUE)) install.packages("clue")
library(clue)## Warning: package 'clue' was built under R version 4.5.1cost <- 1 - abs(Phi_cross) # maximize |φ| <=> minimize (1 - |φ|)
perm <- solve_LSAP(cost) # column index in CF-Quartimax matched to each BI-GEOMIN factor
matched_phi <- Phi_cross[cbind(seq_len(ncol(Lb)), perm)]
# Optional: flip signs in CF-Quartimax to align with BI-GEOMIN (purely for reporting)
sign_flips <- sign(matched_phi); sign_flips[sign_flips == 0] <- 1
Lc_aligned <- Lc[, perm, drop = FALSE] %*% diag(sign_flips, ncol(Lc))
colnames(Lc_aligned) <- paste0("CFQ_F", perm)
# Recompute φ after sign alignment (these will be |φ|)
Phi_aligned <- tucker_matrix(Lb, Lc_aligned)
cat("\n=== Best matches (BI-GEOMIN factor -> CF-Quartimax factor) ===\n")##
## === Best matches (BI-GEOMIN factor -> CF-Quartimax factor) ===match_report <- data.frame(
BIgeo = paste0("F", seq_len(ncol(Lb))),
CFQ = paste0("F", as.integer(perm)),
phi = round(matched_phi, 3),
abs_phi = round(abs(matched_phi), 3),
sign = ifelse(matched_phi >= 0, "+", "-"),
row.names = NULL
)
print(match_report)## BIgeo CFQ phi abs_phi sign
## 1 F1 F2 -0.716 0.716 -
## 2 F2 F4 -0.980 0.980 -
## 3 F3 F3 -0.999 0.999 -
## 4 F4 F1 -0.988 0.988 -cat("\nSummary of matched |φ|:\n")##
## Summary of matched |φ|:summary_abs <- c(
min = round(min(abs(matched_phi)), 3),
q1 = round(quantile(abs(matched_phi), .25), 3),
median= round(median(abs(matched_phi)), 3),
mean = round(mean(abs(matched_phi)), 3),
q3 = round(quantile(abs(matched_phi), .75), 3),
max = round(max(abs(matched_phi)), 3)
)
print(summary_abs)## min q1.25% median mean q3.75% max
## 0.716 0.914 0.984 0.921 0.991 0.999# (Optional) Keep aligned matrices if you’ll compare patterns side-by-side later
L_cfquart_aligned <- Lc_aligned
# - Side-by-side loadings table: BI-GEOMIN vs CF-Quartimax -
# Starting from objects created earlier:
# L_bigeomin, L_cfquart, Phi_cross, perm, L_cfquart_aligned (optional)
Lb <- as.matrix(L_bigeomin)
Lc <- as.matrix(L_cfquart)
# If you didn't run the congruence matching code yet, do a quick alignment now:
if (!exists("Phi_cross") || !exists("perm") || !exists("L_cfquart_aligned")) {
tucker_phi_vec <- function(x, y) sum(x * y) / sqrt(sum(x^2) * sum(y^2))
tucker_matrix <- function(L1, L2) {
p <- ncol(L1)
M <- matrix(NA_real_, nrow = p, ncol = p)
for (i in seq_len(p)) for (j in seq_len(p)) M[i, j] <- tucker_phi_vec(L1[, i], L2[, j])
M
}
Phi_cross <- tucker_matrix(Lb, Lc)
if (!requireNamespace("clue", quietly = TRUE)) install.packages("clue")
library(clue)
cost <- 1 - abs(Phi_cross)
perm <- solve_LSAP(cost)
sign_flips <- sign(Phi_cross[cbind(seq_len(ncol(Lb)), perm)])
sign_flips[sign_flips == 0] <- 1
L_cfquart_aligned <- Lc[, perm, drop = FALSE] %*% diag(sign_flips, ncol(Lc))
colnames(L_cfquart_aligned) <- paste0("CFQ_F", as.integer(perm))
} else {
L_cfquart_aligned <- as.matrix(L_cfquart_aligned)
}
# Build interleaved, side-by-side table
p <- ncol(Lb)
var_names <- rownames(Lb)
if (is.null(var_names)) var_names <- rownames(L_cfquart_aligned)
if (is.null(var_names)) var_names <- make.names(seq_len(nrow(Lb)))
# Round for display
Lb_r <- round(Lb, 3)
Lc_r <- round(L_cfquart_aligned, 3)
# Interleave: BIgeo_F1, CFQ_F1, BIgeo_F2, CFQ_F2, ...
col_list <- list()
for (k in seq_len(p)) {
col_list[[length(col_list) + 1]] <- Lb_r[, k, drop = FALSE]
colnames(col_list[[length(col_list)]]) <- sprintf("BIgeo_F%g", k)
col_list[[length(col_list) + 1]] <- Lc_r[, k, drop = FALSE]
colnames(col_list[[length(col_list)]]) <- sprintf("CFQ_F%g", k)
}
SideBySide <- do.call(cbind, col_list)
SideBySide <- data.frame(Variable = var_names, SideBySide, row.names = NULL, check.names = FALSE)
# Optional helpers: primary factor (by absolute loading) under BI-GEOMIN and CFQ
primary_idx <- function(M) apply(abs(M), 1, which.max)
SideBySide$Primary_BIgeo <- paste0("F", primary_idx(Lb))
SideBySide$Primary_CFQ <- paste0("F", primary_idx(L_cfquart_aligned))
# Print a neat view
print(SideBySide, row.names = FALSE)## Variable BIgeo_F1 CFQ_F1 BIgeo_F2 CFQ_F2 BIgeo_F3 CFQ_F3 BIgeo_F4 CFQ_F4
## rridaw0 -0.665 -0.090 0.000 -0.048 -0.007 -0.028 0.507 0.750
## rridbw0 -0.531 0.004 -0.005 -0.032 -0.109 -0.130 0.456 0.663
## rridcw0 -0.694 -0.011 0.010 -0.026 0.010 -0.013 0.636 0.923
## rriddw0 -0.696 -0.048 0.008 -0.034 0.010 -0.012 0.594 0.869
## rridew0 -0.561 -0.139 -0.023 -0.081 0.054 0.041 0.355 0.537
## rridfw0 -0.551 0.045 -0.690 -0.898 0.028 0.044 0.017 0.046
## rridgw0 -0.537 0.046 -0.716 -0.930 0.009 0.026 -0.021 -0.008
## rridhw0 -0.594 -0.050 -0.626 -0.835 0.028 0.042 -0.001 0.032
## rridiw0 -0.571 -0.121 -0.341 -0.480 -0.107 -0.110 0.082 0.156
## rridjw0 -0.558 0.030 -0.695 -0.907 -0.014 0.001 -0.011 0.009
## rridkw0 -0.830 -0.752 0.102 -0.031 0.108 0.092 0.013 0.137
## rridlw0 -0.376 -0.390 0.060 -0.004 0.059 0.052 -0.037 0.007
## rridmw0 -0.304 -0.141 -0.132 -0.200 -0.278 -0.289 -0.070 -0.065
## rridnw0 -0.253 -0.237 -0.008 -0.057 -0.123 -0.132 -0.082 -0.076
## rridow0 -0.764 -0.641 0.106 -0.002 -0.031 -0.052 0.043 0.165
## rridpw0 -0.715 -0.666 0.065 -0.056 -0.021 -0.037 -0.061 0.020
## rridqw0 -0.861 -0.862 0.169 0.037 0.066 0.045 -0.043 0.072
## rridrw0 -0.801 -0.635 0.162 0.065 0.067 0.046 0.157 0.325
## rridsw0 -0.235 -0.021 0.018 0.020 -0.690 -0.724 0.001 0.023
## rridtw0 -0.177 0.064 0.009 0.027 -0.790 -0.827 0.008 0.022
## rriduw0 -0.261 0.013 0.011 0.019 -0.866 -0.908 0.006 0.029
## rridvw0 -0.282 0.004 -0.009 -0.009 -0.870 -0.912 -0.003 0.019
## rridww0 -0.235 0.044 0.005 0.016 -0.755 -0.791 0.047 0.080
## Primary_BIgeo Primary_CFQ
## F1 F4
## F1 F4
## F1 F4
## F1 F4
## F1 F4
## F2 F2
## F2 F2
## F2 F2
## F1 F2
## F2 F2
## F1 F1
## F1 F1
## F1 F3
## F1 F1
## F1 F1
## F1 F1
## F1 F1
## F1 F1
## F3 F3
## F3 F3
## F3 F3
## F3 F3
## F3 F3# (Optional) Write to CSV for inspection in Excel
# write.csv(SideBySide, file = "RFDSideBySide_Loadings_BIgeo_vs_CFQ.csv", row.names = FALSE)using MplusAutomation to get Bi-CF-Quartimax rotation
You will need to have Mplus installed on the computer where you run this code and will need to have the MplusAutomation package installed.
library(MplusAutomation)## Warning: package 'MplusAutomation' was built under R version 4.5.1## Version: 1.1.1
## We work hard to write this free software. Please help us get credit by citing:
##
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
##
## -- see citation("MplusAutomation").# your data + vars
myvars <- c(
"rridaw0","rridbw0","rridcw0","rriddw0","rridew0","rridfw0","rridgw0",
"rridhw0","rridiw0","rridjw0","rridkw0","rridlw0","rridmw0","rridnw0",
"rridow0","rridpw0","rridqw0","rridrw0","rridsw0","rridtw0","rriduw0",
"rridvw0","rridww0")
dat <- na.omit(AllRFDData[ , myvars])
nf <- 4 # factors
dat <- na.omit(AllRFDData[ , myvars])
# install.packages("MplusAutomation") # if needed
library(MplusAutomation)
# 1) Your analysis data frame
# (replace `dat` with whatever you already built)
# dat <- na.omit(AllRFDData[ , myvars])
# 2) Build the Mplus input *without* DATA: and without NAMES/USEVARIABLES
# You can still add other VARIABLE: statements (e.g., MISSING, CATEGORICAL)
mobj <- mplusObject(
TITLE = "EFA with BI-CF-QUARTIMAX (OBLIQUE);",
# VARIABLE section: include *only* things other than NAMES/USEVARIABLES
VARIABLE = "MISSING ARE ALL (-9999);",
ANALYSIS = paste0(
"ESTIMATOR = ML;\n",
"ROTATION = BI-CF-QUARTIMAX (OBLIQUE);\n",
"TYPE = EFA 4 4;\n" # set factor count here
),
OUTPUT = "SAMPSTAT STANDARDIZED TECH4;",
# >>> This is the key <<<
# Let MplusAutomation write the .dat file and the NAMES/USEVARIABLES lines
rdata = dat,
usevariables = colnames(dat)
)
# 3) Write and run. Do NOT include writeData='never' here since we want it to write the data.
# If Mplus.exe isn’t on PATH, set Mplus_command to its full path.
res <- mplusModeler(
mobj,
modelout = "efa_bicf.inp",
run = 1L,
writeData = "always" # force writing the .dat & adding NAMES/USEVARIABLES
# , Mplus_command = "C:/Program Files/Mplus/Mplus.exe"
)## The file(s)
## 'efa_bicf_9e2ddb076b5b5d995514e3758fc1bbb7.dat'
## currently exist(s) and will be overwritten# 4) Read results (optional)
out <- readModels("efa_bicf.out")
# Try common locations
Lraw <- out$parameters$efa$stdyx$loadings %||%
out$parameters$efa$standardized$loadings %||%
out$parameters$efa$unstandardized$loadings
`%||%` <- function(x, y) if (!is.null(x)) x else y
# Fallback parser keyed on "F1 F2 ..." header
# Robust fixed-width parser for rotated EFA loadings in Mplus .out
parse_mplus_efa_loadings <- function(path, nf = NULL, vars = NULL, take = c("best","last","first")) {
take <- match.arg(take)
lines <- readLines(path, warn = FALSE)
# find all lines that look like factor headers: "F1 F2 ... Fk"
hdr_idx <- which(grepl("^\\s*F\\d+(\\s+F\\d+)+\\s*$", lines))
if (!length(hdr_idx)) stop("No factor header lines found (e.g., 'F1 F2 ...').")
parse_block <- function(s) {
header <- lines[s]
# start positions of each "F#" token in header
hstarts <- as.integer(gregexpr("F\\d+", header)[[1]])
if (any(hstarts < 0)) return(NULL)
# number of factors from header if not supplied
k <- length(hstarts)
if (!is.null(nf) && nf != k) {
# if nf specified and doesn't match, keep going but respect header k
k <- nf
hstarts <- hstarts[seq_len(k)]
}
# compute end positions for each column region
hends <- c(hstarts[-1] - 1L, nchar(header))
# variable name field ends just before first factor column
name_end <- max(1L, hstarts[1] - 1L)
# walk lines below header and slice fixed-width fields
i <- s + 1L
out <- list(); r <- 0L
while (i <= length(lines)) {
ln <- lines[i]
# stop when we hit a blank line or another header/table
if (grepl("^\\s*$", ln)) { if (r > 0) break; i <- i + 1L; next }
if (grepl("^\\s*F\\d+(\\s+F\\d+)+\\s*$", ln)) break
# grab name & k numeric fields by fixed positions
vname <- trimws(substr(ln, 1L, name_end))
# skip non-data rows
if (vname == "" || (!is.null(vars) && !(vname %in% vars))) { i <- i + 1L; next }
vals <- numeric(k)
ok <- TRUE
for (j in seq_len(k)) {
seg <- trimws(substr(ln, hstarts[j], hends[j]))
vals[j] <- suppressWarnings(as.numeric(seg))
if (is.na(vals[j])) vals[j] <- 0 # treat blanks/periods as zero loading
}
if (ok) {
r <- r + 1L
out[[r]] <- data.frame(Variable = vname, t(vals), check.names = FALSE)
}
i <- i + 1L
}
if (r == 0L) return(NULL)
df <- do.call(rbind, out)
names(df) <- c("Variable", paste0("F", seq_len(k)))
# If a var list provided, keep that order and subset to it
if (!is.null(vars)) {
df <- df[df$Variable %in% vars, , drop = FALSE]
df <- df[order(match(df$Variable, vars)), , drop = FALSE]
}
df
}
# try all headers, keep the best block (most rows)
blocks <- lapply(hdr_idx, parse_block)
sizes <- vapply(blocks, function(x) if (is.null(x)) 0L else nrow(x), 1L)
if (max(sizes) == 0L) {
stop("Found header(s) but failed to parse any loading rows. ",
"Check that your .out actually printed rotated loadings.")
}
pick <- switch(take,
best = which.max(sizes),
last = tail(which(sizes > 0L), 1),
first = head(which(sizes > 0L), 1)
)
blocks[[pick]]
}
loadings_wide %>%
mutate(across(-Variable, ~ round(.x, 2))) %>%
print(n = nrow(.))## # A tibble: 23 × 5
## Variable F1 F2 F3 F4
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 rridcw0 0.95 0 -0.01 0
## 2 rriddw0 0.89 0.01 0.03 0
## 3 rridaw0 0.77 0.03 0.07 0.02
## 4 rridbw0 0.68 0.01 -0.02 0.12
## 5 rridew0 0.54 0.07 0.13 -0.04
## 6 rridrw0 0.29 -0.06 0.64 -0.02
## 7 rridiw0 0.14 0.47 0.13 0.13
## 8 rridow0 0.12 0.01 0.66 0.08
## 9 rridnw0 -0.1 0.07 0.25 0.15
## 10 rridmw0 -0.09 0.21 0.16 0.3
## 11 rridkw0 0.08 0.04 0.77 -0.06
## 12 rridww0 0.06 0 -0.03 0.78
## 13 rridfw0 0.04 0.88 -0.03 -0.02
## 14 rridpw0 -0.04 0.07 0.69 0.07
## 15 rridlw0 -0.02 0.01 0.4 -0.03
## 16 rridgw0 -0.02 0.92 -0.03 0
## 17 rridhw0 0.01 0.82 0.07 -0.02
## 18 rridqw0 0.01 -0.02 0.89 -0.01
## 19 rridvw0 -0.01 0.03 0.02 0.91
## 20 rridjw0 -0.01 0.89 -0.01 0.02
## 21 rriduw0 0 0 0.01 0.9
## 22 rridtw0 0 -0.01 -0.05 0.82
## 23 rridsw0 0 -0.01 0.04 0.72Target Rotation (Table 16.8)
R Code
# install.packages("lavaan") # if needed
# install.packages("GPArotation") # if needed
library(lavaan)
myvars <- c(
"rridaw0","rridbw0","rridcw0","rriddw0","rridew0","rridfw0","rridgw0",
"rridhw0","rridiw0","rridjw0","rridkw0","rridlw0","rridmw0","rridnw0",
"rridow0","rridpw0","rridqw0","rridrw0","rridsw0","rridtw0","rriduw0",
"rridvw0","rridww0"
)
dat <- na.omit(AllRFDData[ , myvars])
# 4-factor EFA block
efa_block <- paste(
sprintf('efa("efa4")*F%d =~ %s', 1, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 2, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 3, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 4, paste(myvars, collapse = " + ")),
sep = "\n"
)
# Target matrix: only anchor items constrained to be ~0 on non-target factors
Tmat <- matrix(NA_real_, nrow = length(myvars), ncol = 4,
dimnames = list(myvars, paste0("F", 1:4)))
Tmat["rridaw0", c("F2","F3","F4")] <- 0 # F1 anchor
Tmat["rridgw0", c("F1","F3","F4")] <- 0 # F2 anchor
Tmat["rridqw0", c("F1","F2","F4")] <- 0 # F3 anchor
Tmat["rridvw0", c("F1","F2","F3")] <- 0 # F4 anchor
# sanity check
stopifnot(identical(rownames(Tmat), colnames(dat)))
TargetModel <- cfa(
model = efa_block,
data = dat,
std.lv = TRUE,
estimator = "ml",
rotation = "target",
rotation.args = list(
target = Tmat, # <-- lower-case 'target'
orthogonal = FALSE # oblique target rotation
)
)
summary(TargetModel, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 121
## Row rank of the constraints matrix 12
##
## Rotation method PST OBLIQUE
## Rotation algorithm (rstarts) GPA (30)
## Standardized metric TRUE
## Row weights None
##
## Number of observations 2422
##
## Model Test User Model:
##
## Test statistic 2344.266
## Degrees of freedom 167
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 40456.892
## Degrees of freedom 253
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.946
## Tucker-Lewis Index (TLI) 0.918
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -54129.245
## Loglikelihood unrestricted model (H1) -52957.112
##
## Akaike (AIC) 108476.490
## Bayesian (BIC) 109107.856
## Sample-size adjusted Bayesian (SABIC) 108761.538
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.073
## 90 Percent confidence interval - lower 0.071
## 90 Percent confidence interval - upper 0.076
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.035
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~ efa4
## rridaw0 0.772 0.015 50.224 0.000 0.772 0.836
## rridbw0 0.683 0.028 24.178 0.000 0.683 0.745
## rridcw0 0.963 0.026 36.657 0.000 0.963 1.038
## rriddw0 0.936 0.026 35.559 0.000 0.936 0.974
## rridew0 0.522 0.026 19.791 0.000 0.522 0.593
## rridfw0 0.064 0.024 2.691 0.007 0.064 0.061
## rridgw0 0.000 0.000 0.000
## rridhw0 0.037 0.023 1.575 0.115 0.037 0.036
## rridiw0 0.170 0.029 5.965 0.000 0.170 0.164
## rridjw0 0.017 0.022 0.777 0.437 0.017 0.017
## rridkw0 0.088 0.026 3.359 0.001 0.088 0.085
## rridlw0 -0.029 0.035 -0.817 0.414 -0.029 -0.028
## rridmw0 -0.062 0.022 -2.874 0.004 -0.062 -0.093
## rridnw0 -0.109 0.034 -3.223 0.001 -0.109 -0.111
## rridow0 0.128 0.027 4.666 0.000 0.128 0.124
## rridpw0 -0.047 0.032 -1.455 0.146 -0.047 -0.042
## rridqw0 -0.000 NA -0.000 -0.000
## rridrw0 0.288 0.023 12.640 0.000 0.288 0.306
## rridsw0 0.005 0.020 0.241 0.810 0.005 0.006
## rridtw0 0.007 0.016 0.429 0.668 0.007 0.011
## rriduw0 0.008 0.013 0.600 0.549 0.008 0.012
## rridvw0 0.000 0.000 77.842 0.000 0.000 0.000
## rridww0 0.060 0.020 3.023 0.003 0.060 0.076
## F2 =~ efa4
## rridaw0 0.000 NA 0.000 0.000
## rridbw0 -0.014 0.020 -0.696 0.486 -0.014 -0.015
## rridcw0 -0.034 0.018 -1.888 0.059 -0.034 -0.036
## rriddw0 -0.022 0.018 -1.222 0.222 -0.022 -0.023
## rridew0 0.043 0.019 2.284 0.022 0.043 0.048
## rridfw0 0.894 0.019 46.732 0.000 0.894 0.861
## rridgw0 0.877 0.016 55.719 0.000 0.877 0.895
## rridhw0 0.818 0.019 44.077 0.000 0.818 0.805
## rridiw0 0.473 0.021 22.862 0.000 0.473 0.457
## rridjw0 0.870 0.018 48.331 0.000 0.870 0.872
## rridkw0 0.057 0.018 3.269 0.001 0.057 0.056
## rridlw0 0.023 0.025 0.922 0.356 0.023 0.022
## rridmw0 0.135 0.015 8.786 0.000 0.135 0.201
## rridnw0 0.068 0.024 2.855 0.004 0.068 0.070
## rridow0 0.020 0.019 1.074 0.283 0.020 0.020
## rridpw0 0.093 0.022 4.192 0.000 0.093 0.083
## rridqw0 -0.000 NA -0.000 -0.000
## rridrw0 -0.052 0.016 -3.283 0.001 -0.052 -0.055
## rridsw0 -0.020 0.014 -1.378 0.168 -0.020 -0.026
## rridtw0 -0.024 0.011 -2.123 0.034 -0.024 -0.037
## rriduw0 -0.018 0.009 -1.919 0.055 -0.018 -0.028
## rridvw0 0.000 0.000 77.857 0.000 0.000 0.000
## rridww0 -0.023 0.014 -1.641 0.101 -0.023 -0.029
## F3 =~ efa4
## rridaw0 0.000 0.000 0.000
## rridbw0 -0.079 0.029 -2.694 0.007 -0.079 -0.086
## rridcw0 -0.097 0.027 -3.609 0.000 -0.097 -0.105
## rriddw0 -0.057 0.027 -2.102 0.036 -0.057 -0.059
## rridew0 0.070 0.027 2.570 0.010 0.070 0.080
## rridfw0 -0.008 0.025 -0.317 0.751 -0.008 -0.008
## rridgw0 -0.000 -0.000 -0.000
## rridhw0 0.091 0.024 3.740 0.000 0.091 0.089
## rridiw0 0.136 0.029 4.607 0.000 0.136 0.131
## rridjw0 0.014 0.023 0.605 0.545 0.014 0.014
## rridkw0 0.786 0.028 28.258 0.000 0.786 0.764
## rridlw0 0.418 0.036 11.506 0.000 0.418 0.404
## rridmw0 0.112 0.022 4.983 0.000 0.112 0.166
## rridnw0 0.254 0.035 7.267 0.000 0.254 0.260
## rridow0 0.665 0.029 23.042 0.000 0.665 0.644
## rridpw0 0.774 0.034 22.783 0.000 0.774 0.691
## rridqw0 0.894 0.017 53.352 0.000 0.894 0.883
## rridrw0 0.576 0.024 23.961 0.000 0.576 0.613
## rridsw0 0.017 0.021 0.807 0.419 0.017 0.022
## rridtw0 -0.043 0.017 -2.569 0.010 -0.043 -0.065
## rriduw0 -0.008 0.014 -0.573 0.566 -0.008 -0.012
## rridvw0 0.000 0.000 77.842 0.000 0.000 0.000
## rridww0 -0.041 0.020 -1.980 0.048 -0.041 -0.051
## F4 =~ efa4
## rridaw0 0.000 NA 0.000 0.000
## rridbw0 0.093 0.018 5.152 0.000 0.093 0.101
## rridcw0 -0.026 0.016 -1.619 0.105 -0.026 -0.028
## rriddw0 -0.023 0.017 -1.405 0.160 -0.023 -0.024
## rridew0 -0.049 0.017 -2.911 0.004 -0.049 -0.056
## rridfw0 -0.022 0.015 -1.462 0.144 -0.022 -0.021
## rridgw0 -0.000 NA -0.000 -0.000
## rridhw0 -0.015 0.015 -0.986 0.324 -0.015 -0.015
## rridiw0 0.130 0.018 7.124 0.000 0.130 0.126
## rridjw0 0.024 0.014 1.725 0.084 0.024 0.025
## rridkw0 -0.056 0.016 -3.540 0.000 -0.056 -0.054
## rridlw0 -0.031 0.022 -1.403 0.161 -0.031 -0.030
## rridmw0 0.208 0.014 14.739 0.000 0.208 0.309
## rridnw0 0.148 0.022 6.859 0.000 0.148 0.152
## rridow0 0.085 0.017 4.979 0.000 0.085 0.082
## rridpw0 0.086 0.020 4.294 0.000 0.086 0.077
## rridqw0 -0.000 NA -0.000 -0.000
## rridrw0 -0.025 0.014 -1.740 0.082 -0.025 -0.026
## rridsw0 0.547 0.014 38.186 0.000 0.547 0.728
## rridtw0 0.541 0.012 45.241 0.000 0.541 0.827
## rriduw0 0.582 0.011 53.889 0.000 0.582 0.911
## rridvw0 0.603 0.010 58.075 0.000 0.603 0.917
## rridww0 0.624 0.015 43.026 0.000 0.624 0.790
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.412 0.023 17.769 0.000 0.412 0.412
## F3 0.760 0.016 48.771 0.000 0.760 0.760
## F4 0.225 0.025 8.933 0.000 0.225 0.225
## F2 ~~
## F3 0.425 0.022 19.362 0.000 0.425 0.425
## F4 0.348 0.022 15.697 0.000 0.348 0.348
## F3 ~~
## F4 0.190 0.024 7.855 0.000 0.190 0.190
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw0 0.258 0.008 30.547 0.000 0.258 0.302
## .rridbw0 0.424 0.013 32.741 0.000 0.424 0.503
## .rridcw0 0.097 0.006 16.326 0.000 0.097 0.113
## .rriddw0 0.148 0.007 22.485 0.000 0.148 0.161
## .rridew0 0.432 0.013 33.614 0.000 0.432 0.557
## .rridfw0 0.247 0.009 26.394 0.000 0.247 0.230
## .rridgw0 0.191 0.008 24.294 0.000 0.191 0.199
## .rridhw0 0.269 0.009 28.411 0.000 0.269 0.261
## .rridiw0 0.569 0.017 33.573 0.000 0.569 0.531
## .rridjw0 0.199 0.008 24.908 0.000 0.199 0.200
## .rridkw0 0.301 0.011 26.712 0.000 0.301 0.284
## .rridlw0 0.909 0.027 34.086 0.000 0.909 0.850
## .rridmw0 0.357 0.010 34.250 0.000 0.357 0.788
## .rridnw0 0.874 0.025 34.438 0.000 0.874 0.913
## .rridow0 0.429 0.014 30.747 0.000 0.429 0.403
## .rridpw0 0.604 0.019 31.014 0.000 0.604 0.482
## .rridqw0 0.226 0.011 21.182 0.000 0.226 0.221
## .rridrw0 0.257 0.009 29.118 0.000 0.257 0.292
## .rridsw0 0.267 0.008 32.105 0.000 0.267 0.474
## .rridtw0 0.148 0.005 29.966 0.000 0.148 0.347
## .rriduw0 0.076 0.003 23.462 0.000 0.076 0.186
## .rridvw0 0.069 0.003 21.567 0.000 0.069 0.159
## .rridww0 0.236 0.008 30.751 0.000 0.236 0.378
## F1 1.000 1.000 1.000
## F2 1.000 1.000 1.000
## F3 1.000 1.000 1.000
## F4 1.000 1.000 1.000lavInspect(TargetModel, "cor.lv") # factor correlations (Phi)## F1 F2 F3 F4
## F1 1.000
## F2 0.412 1.000
## F3 0.760 0.425 1.000
## F4 0.225 0.348 0.190 1.000# Tidy standardized loadings
loadings_df <- standardizedSolution(TargetModel)
loadings_df <- subset(loadings_df, op == "=~", select = c("lhs","rhs","est.std"))
names(loadings_df) <- c("Factor","Variable","Loading")
loadings_df <- loadings_df[order(loadings_df$Factor, -abs(loadings_df$Loading)), ]
library(dplyr)
library(tidyr)
# Wide table of standardized loadings (F1..F4 as columns)
loadings_wide <- standardizedSolution(TargetModel) %>%
filter(op == "=~") %>% # keep loadings only
transmute(Variable = rhs, Factor = lhs, Loading = est.std) %>%
mutate(Factor = factor(Factor, levels = paste0("F", 1:4))) %>%
arrange(match(Variable, myvars), Factor) %>%
pivot_wider(names_from = Factor, values_from = Loading) %>%
arrange(match(Variable, myvars))
# Round to two decimals for display
loadings_wide_print <- loadings_wide %>%
mutate(across(starts_with("F"), ~ round(.x, 2)))
# Print all rows
print(loadings_wide_print, row.names = FALSE)## # A tibble: 23 × 5
## Variable F1 F2 F3 F4
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 rridaw0 0.84 0 0 0
## 2 rridbw0 0.74 -0.02 -0.09 0.1
## 3 rridcw0 1.04 -0.04 -0.1 -0.03
## 4 rriddw0 0.97 -0.02 -0.06 -0.02
## 5 rridew0 0.59 0.05 0.08 -0.06
## 6 rridfw0 0.06 0.86 -0.01 -0.02
## 7 rridgw0 0 0.9 0 0
## 8 rridhw0 0.04 0.81 0.09 -0.01
## 9 rridiw0 0.16 0.46 0.13 0.13
## 10 rridjw0 0.02 0.87 0.01 0.02
## # ℹ 13 more rows# install.packages(c("knitr","kableExtra"))
library(dplyr)
library(knitr)
library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.5.1##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows# if you have the numeric version `loadings_wide`, start from it:
tbl <- loadings_wide %>%
mutate(across(starts_with("F"), ~ sprintf("%.2f", .x))) # fixed 2 decimals
# (optional) bold |loading| >= .30
fmt <- if (knitr::is_latex_output()) "latex" else "html"
tbl <- tbl %>%
mutate(across(starts_with("F"),
~ ifelse(abs(as.numeric(.x)) >= .30,
kableExtra::cell_spec(.x, format = fmt, bold = TRUE),
.x)
))
kbl(tbl,
booktabs = TRUE,
align = c("l","r","r","r","r"),
caption = "Table 16.8\nStandardized Factor Loadings from a 4-Factor Target Rotation",
escape = FALSE # allow bold from cell_spec
) %>%
kable_classic(full_width = FALSE, html_font = "Times New Roman") %>%
add_header_above(c(" " = 1, "Factors" = 4)) %>%
footnote(
general = "Note. Values are standardized loadings (est.std). Bold indicates |loading| ≥ .30.",
threeparttable = TRUE
)|
Factors
|
||||
|---|---|---|---|---|
| Variable | F1 | F2 | F3 | F4 |
| rridaw0 | 0.84 | 0.00 | 0.00 | 0.00 |
| rridbw0 | 0.74 | -0.02 | -0.09 | 0.10 |
| rridcw0 | 1.04 | -0.04 | -0.10 | -0.03 |
| rriddw0 | 0.97 | -0.02 | -0.06 | -0.02 |
| rridew0 | 0.59 | 0.05 | 0.08 | -0.06 |
| rridfw0 | 0.06 | 0.86 | -0.01 | -0.02 |
| rridgw0 | 0.00 | 0.90 | -0.00 | -0.00 |
| rridhw0 | 0.04 | 0.81 | 0.09 | -0.01 |
| rridiw0 | 0.16 | 0.46 | 0.13 | 0.13 |
| rridjw0 | 0.02 | 0.87 | 0.01 | 0.02 |
| rridkw0 | 0.09 | 0.06 | 0.76 | -0.05 |
| rridlw0 | -0.03 | 0.02 | 0.40 | -0.03 |
| rridmw0 | -0.09 | 0.20 | 0.17 | 0.31 |
| rridnw0 | -0.11 | 0.07 | 0.26 | 0.15 |
| rridow0 | 0.12 | 0.02 | 0.64 | 0.08 |
| rridpw0 | -0.04 | 0.08 | 0.69 | 0.08 |
| rridqw0 | -0.00 | -0.00 | 0.88 | -0.00 |
| rridrw0 | 0.31 | -0.05 | 0.61 | -0.03 |
| rridsw0 | 0.01 | -0.03 | 0.02 | 0.73 |
| rridtw0 | 0.01 | -0.04 | -0.07 | 0.83 |
| rriduw0 | 0.01 | -0.03 | -0.01 | 0.91 |
| rridvw0 | 0.00 | 0.00 | 0.00 | 0.92 |
| rridww0 | 0.08 | -0.03 | -0.05 | 0.79 |
| Note: | ||||
| Note. Values are standardized loadings (est.std). Bold indicates |loading| ≥ .30. | ||||
require(lavaangui)## Loading required package: lavaangui## Warning: package 'lavaangui' was built under R version 4.5.1## This is lavaangui 0.2.5
## lavaangui is BETA software! Please report any bugs at https://github.com/karchjd/lavaangui/issues#plot_lavaan(fit)Multi-Group Invariance with ESEM Model C (Intercepts and Loadings Invariant)
Download Simulated Item-level data
<a https://docs.google.com/document/d/1-64WYwrcii6y0T7yxxcTYqrLXzXiUUjKyh7L6E8sMXc/edit?usp=sharing
Unstandardized factor loadings, factor means and factor variances by group andstandardized factor loadings
Remember that, because the factor variances differ by group, the standardized loadings are not going to be identical across group, only the unstandardized loadings will be.
# install.packages("lavaan") # if needed
library(lavaan)
# --- data ---
RFDItems <- read.table(
"RFDItems.txt",
sep="\t", header=TRUE)
RFDItems$sex <- factor(RFDItems$sex)
myvars <- c(
"rridaw1","rridbw1","rridcw1","rriddw1","rridew1","rridfw1","rridgw1",
"rridhw1","rridiw1","rridjw1","rridkw1","rridlw1","rridmw1","rridnw1",
"rridow1","rridpw1","rridqw1","rridrw1","rridsw1","rridtw1","rriduw1",
"rridvw1","rridww1"
)
# --- ESEM model with a 4-factor EFA block in each group ---
efa_block <- paste(
sprintf('efa("efa4")*F%d =~ %s', 1, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 2, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 3, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 4, paste(myvars, collapse = " + ")),
sep = "\n"
)
# --- Minimal shared target to anchor factor orientation across groups ---
# (Zeros on non-target columns for four anchor items; all else NA = don't care)
T_shared <- matrix(NA_real_, nrow = length(myvars), ncol = 4,
dimnames = list(myvars, paste0("F", 1:4)))
T_shared["rridaw1", c("F2","F3","F4")] <- 0 # anchor for F1
T_shared["rridgw1", c("F1","F3","F4")] <- 0 # anchor for F2
T_shared["rridqw1", c("F1","F2","F4")] <- 0 # anchor for F3
T_shared["rridvw1", c("F1","F2","F3")] <- 0 # anchor for F4
# --- Fit: MG-ESEM with loadings+intercepts invariant; other things variant ---
fit_mg_esemC <- cfa(
model = efa_block,
data = RFDItems,
group = "sex",
meanstructure = TRUE, # needed for intercept/mean invariance
estimator = "ml",
# rotation: target (oblique), same target in both groups
rotation = "target",
rotation.args = list(target = T_shared, orthogonal = FALSE),
# equality constraints ACROSS groups:
group.equal = c("loadings", "intercepts")
# (Note: we are NOT constraining lv.variances, lv.covariances, or residuals,
# so factor means/variances/covariances and error variances are free by group.)
)
# --- Inspect results ---
summary(fit_mg_esemC, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 208 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 296
## Number of equality constraints 115
## Row rank of the constraints matrix 127
##
## Rotation method PST OBLIQUE
## Rotation algorithm (rstarts) GPA (30)
## Standardized metric TRUE
## Row weights None
##
## Number of observations per group:
## 1 865
## 2 1294
##
## Model Test User Model:
##
## Test statistic 2812.116
## Degrees of freedom 429
## P-value (Chi-square) 0.000
## Test statistic for each group:
## 1 1363.111
## 2 1449.004
##
## Model Test Baseline Model:
##
## Test statistic 37163.568
## Degrees of freedom 506
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.935
## Tucker-Lewis Index (TLI) 0.923
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -47308.523
## Loglikelihood unrestricted model (H1) -45902.465
##
## Akaike (AIC) 94955.046
## Bayesian (BIC) 95914.527
## Sample-size adjusted Bayesian (SABIC) 95377.592
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.072
## 90 Percent confidence interval - lower 0.069
## 90 Percent confidence interval - upper 0.074
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.044
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
##
## Group 1 [1]:
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~ efa4
## rridaw1 (.p1.) 0.709 0.021 34.405 0.000 0.709 0.819
## rridbw1 (.p2.) 0.643 0.027 23.599 0.000 0.643 0.771
## rridcw1 (.p3.) 0.859 0.029 29.750 0.000 0.859 1.008
## rriddw1 (.p4.) 0.865 0.029 29.644 0.000 0.865 0.990
## rridew1 (.p5.) 0.518 0.026 19.990 0.000 0.518 0.646
## rridfw1 (.p6.) 0.104 0.021 4.837 0.000 0.104 0.101
## rridgw1 (.p7.) 0.000 0.000 0.000
## rridhw1 (.p8.) 0.064 0.021 2.999 0.003 0.064 0.064
## rridiw1 (.p9.) 0.284 0.028 10.157 0.000 0.284 0.273
## rridjw1 (.10.) 0.033 0.021 1.586 0.113 0.033 0.033
## rridkw1 (.11.) 0.192 0.024 7.905 0.000 0.192 0.202
## rridlw1 (.12.) -0.002 0.031 -0.079 0.937 -0.002 -0.002
## rridmw1 (.13.) -0.151 0.019 -7.913 0.000 -0.151 -0.192
## rridnw1 (.14.) -0.134 0.033 -4.091 0.000 -0.134 -0.128
## rridow1 (.15.) 0.047 0.027 1.759 0.079 0.047 0.047
## rridpw1 (.16.) -0.025 0.030 -0.815 0.415 -0.025 -0.022
## rridqw1 (.17.) -0.000 -0.000 -0.000
## rridrw1 (.18.) 0.241 0.022 10.892 0.000 0.241 0.285
## rridsw1 (.19.) -0.031 0.018 -1.679 0.093 -0.031 -0.035
## rridtw1 (.20.) -0.030 0.015 -2.011 0.044 -0.030 -0.037
## rriduw1 (.21.) -0.017 0.013 -1.333 0.182 -0.017 -0.021
## rridvw1 (.22.) 0.000 0.000 0.000
## rridww1 (.23.) 0.064 0.020 3.282 0.001 0.064 0.071
## F2 =~ efa4
## rridaw1 0.000 0.000 50.175 0.000 0.000 0.000
## rridbw1 (.25.) 0.001 0.019 0.077 0.939 0.001 0.002
## rridcw1 (.26.) -0.051 0.018 -2.875 0.004 -0.051 -0.060
## rriddw1 (.27.) -0.041 0.018 -2.258 0.024 -0.041 -0.047
## rridew1 (.28.) 0.003 0.019 0.169 0.866 0.003 0.004
## rridfw1 (.29.) 0.901 0.027 33.653 0.000 0.901 0.877
## rridgw1 (.30.) 0.908 0.025 36.936 0.000 0.908 0.914
## rridhw1 (.31.) 0.844 0.026 32.879 0.000 0.844 0.843
## rridiw1 (.32.) 0.452 0.024 18.903 0.000 0.452 0.436
## rridjw1 (.33.) 0.898 0.026 33.944 0.000 0.898 0.886
## rridkw1 (.34.) 0.110 0.018 5.970 0.000 0.110 0.115
## rridlw1 (.35.) 0.064 0.024 2.619 0.009 0.064 0.064
## rridmw1 (.36.) 0.086 0.015 5.675 0.000 0.086 0.109
## rridnw1 (.37.) -0.027 0.026 -1.066 0.287 -0.027 -0.026
## rridow1 (.38.) 0.039 0.020 1.946 0.052 0.039 0.040
## rridpw1 (.39.) 0.092 0.023 4.001 0.000 0.092 0.084
## rridqw1 (.40.) -0.000 -0.000 -0.000
## rridrw1 (.41.) -0.062 0.016 -3.813 0.000 -0.062 -0.073
## rridsw1 (.42.) -0.003 0.015 -0.198 0.843 -0.003 -0.003
## rridtw1 (.43.) -0.003 0.012 -0.268 0.788 -0.003 -0.004
## rriduw1 (.44.) -0.009 0.010 -0.878 0.380 -0.009 -0.011
## rridvw1 (.45.) -0.000 -0.000 -0.000
## rridww1 (.46.) -0.027 0.015 -1.741 0.082 -0.027 -0.030
## F3 =~ efa4
## rridaw1 -0.000 -0.000 -0.000
## rridbw1 -0.041 0.025 -1.649 0.099 -0.041 -0.049
## rridcw1 (.49.) -0.059 0.024 -2.489 0.013 -0.059 -0.069
## rriddw1 (.50.) -0.038 0.024 -1.573 0.116 -0.038 -0.043
## rridew1 (.51.) 0.049 0.025 2.014 0.044 0.049 0.062
## rridfw1 (.52.) -0.034 0.022 -1.527 0.127 -0.034 -0.033
## rridgw1 (.53.) -0.000 -0.000 -0.000
## rridhw1 (.54.) 0.057 0.022 2.587 0.010 0.057 0.057
## rridiw1 (.55.) 0.091 0.028 3.222 0.001 0.091 0.088
## rridjw1 (.56.) 0.006 0.021 0.301 0.763 0.006 0.006
## rridkw1 (.57.) 0.572 0.028 20.547 0.000 0.572 0.601
## rridlw1 (.58.) 0.325 0.033 9.862 0.000 0.325 0.324
## rridmw1 (.59.) 0.154 0.020 7.718 0.000 0.154 0.196
## rridnw1 (.60.) 0.285 0.034 8.276 0.000 0.285 0.272
## rridow1 (.61.) 0.651 0.031 20.790 0.000 0.651 0.652
## rridpw1 (.62.) 0.745 0.036 20.967 0.000 0.745 0.676
## rridqw1 (.63.) 0.811 0.023 34.626 0.000 0.811 0.896
## rridrw1 (.64.) 0.537 0.025 21.442 0.000 0.537 0.636
## rridsw1 (.65.) 0.048 0.019 2.506 0.012 0.048 0.054
## rridtw1 (.66.) -0.015 0.015 -0.984 0.325 -0.015 -0.019
## rriduw1 (.67.) 0.011 0.013 0.840 0.401 0.011 0.014
## rridvw1 (.68.) -0.000 -0.000 -0.000
## rridww1 (.69.) -0.032 0.020 -1.610 0.107 -0.032 -0.036
## F4 =~ efa4
## rridaw1 0.000 0.000 50.176 0.000 0.000 0.000
## rridbw1 0.012 0.020 0.601 0.548 0.012 0.014
## rridcw1 -0.084 0.019 -4.415 0.000 -0.084 -0.099
## rriddw1 (.73.) -0.083 0.019 -4.281 0.000 -0.083 -0.095
## rridew1 (.74.) -0.102 0.020 -5.102 0.000 -0.102 -0.127
## rridfw1 (.75.) -0.058 0.018 -3.202 0.001 -0.058 -0.056
## rridgw1 (.76.) -0.000 -0.000 -0.000
## rridhw1 (.77.) -0.031 0.018 -1.726 0.084 -0.031 -0.031
## rridiw1 (.78.) 0.071 0.023 3.065 0.002 0.071 0.068
## rridjw1 (.79.) -0.032 0.018 -1.809 0.070 -0.032 -0.031
## rridkw1 (.80.) -0.120 0.019 -6.191 0.000 -0.120 -0.127
## rridlw1 (.81.) -0.115 0.026 -4.389 0.000 -0.115 -0.115
## rridmw1 (.82.) 0.342 0.018 18.739 0.000 0.342 0.435
## rridnw1 (.83.) 0.226 0.028 8.044 0.000 0.226 0.215
## rridow1 (.84.) 0.090 0.022 4.162 0.000 0.090 0.090
## rridpw1 (.85.) 0.116 0.024 4.738 0.000 0.116 0.105
## rridqw1 (.86.) -0.000 -0.000 -0.000
## rridrw1 (.87.) -0.082 0.017 -4.800 0.000 -0.082 -0.097
## rridsw1 (.88.) 0.696 0.022 31.065 0.000 0.696 0.792
## rridtw1 (.89.) 0.701 0.021 34.005 0.000 0.701 0.862
## rriduw1 (.90.) 0.761 0.021 36.435 0.000 0.761 0.925
## rridvw1 (.91.) 0.769 0.020 37.672 0.000 0.769 0.930
## rridww1 (.92.) 0.784 0.024 32.945 0.000 0.784 0.871
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.326 0.036 8.972 0.000 0.326 0.326
## F3 0.723 0.023 31.130 0.000 0.723 0.723
## F4 0.217 0.038 5.654 0.000 0.217 0.217
## F2 ~~
## F3 0.245 0.037 6.600 0.000 0.245 0.245
## F4 0.458 0.031 14.966 0.000 0.458 0.458
## F3 ~~
## F4 0.111 0.038 2.906 0.004 0.111 0.111
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw1 (.126) 3.006 0.027 112.000 0.000 3.006 3.472
## .rridbw1 (.127) 3.031 0.025 123.234 0.000 3.031 3.636
## .rridcw1 (.128) 3.125 0.028 110.955 0.000 3.125 3.669
## .rriddw1 (.129) 3.062 0.029 105.717 0.000 3.062 3.504
## .rridew1 (.130) 3.156 0.023 137.211 0.000 3.156 3.932
## .rridfw1 (.131) 2.039 0.033 61.604 0.000 2.039 1.985
## .rridgw1 (.132) 1.929 0.033 58.938 0.000 1.929 1.941
## .rridhw1 (.133) 2.059 0.032 64.122 0.000 2.059 2.057
## .rridiw1 (.134) 2.331 0.029 79.817 0.000 2.331 2.246
## .rridjw1 (.135) 1.957 0.033 59.784 0.000 1.957 1.932
## .rridkw1 (.136) 3.046 0.029 104.643 0.000 3.046 3.198
## .rridlw1 (.137) 2.736 0.023 116.675 0.000 2.736 2.730
## .rridmw1 (.138) 1.344 0.018 74.127 0.000 1.344 1.712
## .rridnw1 (.139) 1.901 0.024 80.003 0.000 1.901 1.812
## .rridow1 (.140) 2.664 0.029 91.063 0.000 2.664 2.671
## .rridpw1 (.141) 2.434 0.032 75.379 0.000 2.434 2.208
## .rridqw1 (.142) 2.976 0.030 99.285 0.000 2.976 3.286
## .rridrw1 (.143) 3.199 0.027 118.616 0.000 3.199 3.788
## .rridsw1 (.144) 1.576 0.026 59.821 0.000 1.576 1.794
## .rridtw1 (.145) 1.437 0.025 56.744 0.000 1.437 1.767
## .rriduw1 (.146) 1.452 0.027 54.104 0.000 1.452 1.766
## .rridvw1 (.147) 1.481 0.027 54.226 0.000 1.481 1.790
## .rridww1 (.148) 1.640 0.029 56.550 0.000 1.640 1.823
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw1 0.247 0.013 18.673 0.000 0.247 0.330
## .rridbw1 0.314 0.016 19.432 0.000 0.314 0.451
## .rridcw1 0.102 0.008 12.821 0.000 0.102 0.140
## .rriddw1 0.103 0.008 12.782 0.000 0.103 0.134
## .rridew1 0.349 0.017 19.976 0.000 0.349 0.541
## .rridfw1 0.237 0.015 16.185 0.000 0.237 0.224
## .rridgw1 0.163 0.012 14.040 0.000 0.163 0.165
## .rridhw1 0.242 0.014 16.895 0.000 0.242 0.242
## .rridiw1 0.598 0.030 20.171 0.000 0.598 0.555
## .rridjw1 0.222 0.014 15.964 0.000 0.222 0.216
## .rridkw1 0.350 0.019 18.490 0.000 0.350 0.385
## .rridlw1 0.888 0.043 20.524 0.000 0.888 0.884
## .rridmw1 0.466 0.023 20.397 0.000 0.466 0.755
## .rridnw1 1.012 0.049 20.568 0.000 1.012 0.919
## .rridow1 0.483 0.026 18.825 0.000 0.483 0.486
## .rridpw1 0.604 0.032 18.731 0.000 0.604 0.497
## .rridqw1 0.162 0.014 11.541 0.000 0.162 0.197
## .rridrw1 0.209 0.012 17.149 0.000 0.209 0.293
## .rridsw1 0.290 0.015 19.062 0.000 0.290 0.376
## .rridtw1 0.182 0.010 17.995 0.000 0.182 0.275
## .rriduw1 0.108 0.007 15.098 0.000 0.108 0.159
## .rridvw1 0.093 0.007 14.022 0.000 0.093 0.136
## .rridww1 0.195 0.011 17.376 0.000 0.195 0.242
## F1 1.000 1.000 1.000
## F2 1.000 1.000 1.000
## F3 1.000 1.000 1.000
## F4 1.000 1.000 1.000
##
##
## Group 2 [2]:
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~ efa4
## rridaw1 (.p1.) 0.709 0.021 34.405 0.000 0.746 0.827
## rridbw1 (.p2.) 0.643 0.027 23.599 0.000 0.677 0.793
## rridcw1 (.p3.) 0.859 0.029 29.750 0.000 0.904 1.033
## rriddw1 (.p4.) 0.865 0.029 29.644 0.000 0.911 0.997
## rridew1 (.p5.) 0.518 0.026 19.990 0.000 0.546 0.650
## rridfw1 (.p6.) 0.104 0.021 4.837 0.000 0.109 0.106
## rridgw1 (.p7.) 0.000 0.000 0.000
## rridhw1 (.p8.) 0.064 0.021 2.999 0.003 0.068 0.066
## rridiw1 (.p9.) 0.284 0.028 10.157 0.000 0.299 0.288
## rridjw1 (.10.) 0.033 0.021 1.586 0.113 0.035 0.034
## rridkw1 (.11.) 0.192 0.024 7.905 0.000 0.203 0.210
## rridlw1 (.12.) -0.002 0.031 -0.079 0.937 -0.003 -0.003
## rridmw1 (.13.) -0.151 0.019 -7.913 0.000 -0.159 -0.277
## rridnw1 (.14.) -0.134 0.033 -4.091 0.000 -0.141 -0.146
## rridow1 (.15.) 0.047 0.027 1.759 0.079 0.050 0.051
## rridpw1 (.16.) -0.025 0.030 -0.815 0.415 -0.026 -0.024
## rridqw1 (.17.) -0.000 -0.000 -0.000
## rridrw1 (.18.) 0.241 0.022 10.892 0.000 0.254 0.282
## rridsw1 (.19.) -0.031 0.018 -1.679 0.093 -0.033 -0.046
## rridtw1 (.20.) -0.030 0.015 -2.011 0.044 -0.031 -0.050
## rriduw1 (.21.) -0.017 0.013 -1.333 0.182 -0.018 -0.029
## rridvw1 (.22.) -0.000 -0.000 -0.000
## rridww1 (.23.) 0.064 0.020 3.282 0.001 0.068 0.083
## F2 =~ efa4
## rridaw1 0.000 0.000 3.662 0.000 0.000 0.000
## rridbw1 (.25.) 0.001 0.019 0.077 0.939 0.001 0.002
## rridcw1 (.26.) -0.051 0.018 -2.875 0.004 -0.051 -0.058
## rriddw1 (.27.) -0.041 0.018 -2.258 0.024 -0.041 -0.045
## rridew1 (.28.) 0.003 0.019 0.169 0.866 0.003 0.004
## rridfw1 (.29.) 0.901 0.027 33.653 0.000 0.891 0.867
## rridgw1 (.30.) 0.908 0.025 36.936 0.000 0.898 0.903
## rridhw1 (.31.) 0.844 0.026 32.879 0.000 0.835 0.820
## rridiw1 (.32.) 0.452 0.024 18.903 0.000 0.447 0.430
## rridjw1 (.33.) 0.898 0.026 33.944 0.000 0.888 0.881
## rridkw1 (.34.) 0.110 0.018 5.970 0.000 0.109 0.112
## rridlw1 (.35.) 0.064 0.024 2.619 0.009 0.063 0.067
## rridmw1 (.36.) 0.086 0.015 5.675 0.000 0.085 0.148
## rridnw1 (.37.) -0.027 0.026 -1.066 0.287 -0.027 -0.028
## rridow1 (.38.) 0.039 0.020 1.946 0.052 0.039 0.040
## rridpw1 (.39.) 0.092 0.023 4.001 0.000 0.091 0.084
## rridqw1 (.40.) -0.000 -0.000 -0.000
## rridrw1 (.41.) -0.062 0.016 -3.813 0.000 -0.061 -0.068
## rridsw1 (.42.) -0.003 0.015 -0.198 0.843 -0.003 -0.004
## rridtw1 (.43.) -0.003 0.012 -0.268 0.788 -0.003 -0.005
## rriduw1 (.44.) -0.009 0.010 -0.878 0.380 -0.009 -0.014
## rridvw1 (.45.) -0.000 -0.000 -0.000
## rridww1 (.46.) -0.027 0.015 -1.741 0.082 -0.026 -0.032
## F3 =~ efa4
## rridaw1 -0.000 -0.000 -0.000
## rridbw1 -0.041 0.025 -1.649 0.099 -0.041 -0.048
## rridcw1 (.49.) -0.059 0.024 -2.489 0.013 -0.059 -0.067
## rriddw1 (.50.) -0.038 0.024 -1.573 0.116 -0.038 -0.042
## rridew1 (.51.) 0.049 0.025 2.014 0.044 0.049 0.059
## rridfw1 (.52.) -0.034 0.022 -1.527 0.127 -0.034 -0.033
## rridgw1 (.53.) 0.000 0.000 0.000
## rridhw1 (.54.) 0.057 0.022 2.587 0.010 0.057 0.056
## rridiw1 (.55.) 0.091 0.028 3.222 0.001 0.091 0.088
## rridjw1 (.56.) 0.006 0.021 0.301 0.763 0.006 0.006
## rridkw1 (.57.) 0.572 0.028 20.547 0.000 0.573 0.594
## rridlw1 (.58.) 0.325 0.033 9.862 0.000 0.325 0.342
## rridmw1 (.59.) 0.154 0.020 7.718 0.000 0.154 0.269
## rridnw1 (.60.) 0.285 0.034 8.276 0.000 0.285 0.294
## rridow1 (.61.) 0.651 0.031 20.790 0.000 0.651 0.665
## rridpw1 (.62.) 0.745 0.036 20.967 0.000 0.746 0.685
## rridqw1 (.63.) 0.811 0.023 34.626 0.000 0.813 0.854
## rridrw1 (.64.) 0.537 0.025 21.442 0.000 0.538 0.598
## rridsw1 (.65.) 0.048 0.019 2.506 0.012 0.048 0.067
## rridtw1 (.66.) -0.015 0.015 -0.984 0.325 -0.015 -0.024
## rriduw1 (.67.) 0.011 0.013 0.840 0.401 0.011 0.018
## rridvw1 (.68.) -0.000 -0.000 -0.000
## rridww1 (.69.) -0.032 0.020 -1.610 0.107 -0.032 -0.040
## F4 =~ efa4
## rridaw1 0.000 0.000 3.662 0.000 0.000 0.000
## rridbw1 0.012 0.020 0.601 0.548 0.009 0.011
## rridcw1 -0.084 0.019 -4.415 0.000 -0.065 -0.074
## rriddw1 (.73.) -0.083 0.019 -4.281 0.000 -0.064 -0.071
## rridew1 (.74.) -0.102 0.020 -5.102 0.000 -0.079 -0.094
## rridfw1 (.75.) -0.058 0.018 -3.202 0.001 -0.045 -0.043
## rridgw1 (.76.) 0.000 0.000 0.000
## rridhw1 (.77.) -0.031 0.018 -1.726 0.084 -0.024 -0.024
## rridiw1 (.78.) 0.071 0.023 3.065 0.002 0.055 0.053
## rridjw1 (.79.) -0.032 0.018 -1.809 0.070 -0.025 -0.024
## rridkw1 (.80.) -0.120 0.019 -6.191 0.000 -0.093 -0.097
## rridlw1 (.81.) -0.115 0.026 -4.389 0.000 -0.089 -0.094
## rridmw1 (.82.) 0.342 0.018 18.739 0.000 0.265 0.463
## rridnw1 (.83.) 0.226 0.028 8.044 0.000 0.175 0.181
## rridow1 (.84.) 0.090 0.022 4.162 0.000 0.070 0.071
## rridpw1 (.85.) 0.116 0.024 4.738 0.000 0.090 0.083
## rridqw1 (.86.) -0.000 -0.000 -0.000
## rridrw1 (.87.) -0.082 0.017 -4.800 0.000 -0.063 -0.070
## rridsw1 (.88.) 0.696 0.022 31.065 0.000 0.540 0.759
## rridtw1 (.89.) 0.701 0.021 34.005 0.000 0.543 0.864
## rriduw1 (.90.) 0.761 0.021 36.435 0.000 0.590 0.934
## rridvw1 (.91.) 0.769 0.020 37.672 0.000 0.596 0.905
## rridww1 (.92.) 0.784 0.024 32.945 0.000 0.608 0.742
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.417 0.040 10.462 0.000 0.400 0.400
## F3 0.727 0.052 14.004 0.000 0.689 0.689
## F4 0.204 0.029 7.160 0.000 0.250 0.250
## F2 ~~
## F3 0.330 0.037 9.032 0.000 0.333 0.333
## F4 0.285 0.028 10.214 0.000 0.371 0.371
## F3 ~~
## F4 0.133 0.027 5.013 0.000 0.171 0.171
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw1 (.126) 3.006 0.027 112.000 0.000 3.006 3.331
## .rridbw1 (.127) 3.031 0.025 123.234 0.000 3.031 3.551
## .rridcw1 (.128) 3.125 0.028 110.955 0.000 3.125 3.572
## .rriddw1 (.129) 3.062 0.029 105.717 0.000 3.062 3.352
## .rridew1 (.130) 3.156 0.023 137.211 0.000 3.156 3.757
## .rridfw1 (.131) 2.039 0.033 61.604 0.000 2.039 1.984
## .rridgw1 (.132) 1.929 0.033 58.938 0.000 1.929 1.938
## .rridhw1 (.133) 2.059 0.032 64.122 0.000 2.059 2.023
## .rridiw1 (.134) 2.331 0.029 79.817 0.000 2.331 2.243
## .rridjw1 (.135) 1.957 0.033 59.784 0.000 1.957 1.942
## .rridkw1 (.136) 3.046 0.029 104.643 0.000 3.046 3.155
## .rridlw1 (.137) 2.736 0.023 116.675 0.000 2.736 2.878
## .rridmw1 (.138) 1.344 0.018 74.127 0.000 1.344 2.349
## .rridnw1 (.139) 1.901 0.024 80.003 0.000 1.901 1.959
## .rridow1 (.140) 2.664 0.029 91.063 0.000 2.664 2.719
## .rridpw1 (.141) 2.434 0.032 75.379 0.000 2.434 2.233
## .rridqw1 (.142) 2.976 0.030 99.285 0.000 2.976 3.128
## .rridrw1 (.143) 3.199 0.027 118.616 0.000 3.199 3.559
## .rridsw1 (.144) 1.576 0.026 59.821 0.000 1.576 2.217
## .rridtw1 (.145) 1.437 0.025 56.744 0.000 1.437 2.285
## .rriduw1 (.146) 1.452 0.027 54.104 0.000 1.452 2.298
## .rridvw1 (.147) 1.481 0.027 54.226 0.000 1.481 2.247
## .rridww1 (.148) 1.640 0.029 56.550 0.000 1.640 2.001
## F1 -0.127 0.046 -2.736 0.006 -0.120 -0.120
## F2 -0.055 0.045 -1.208 0.227 -0.055 -0.055
## F3 -0.227 0.047 -4.807 0.000 -0.226 -0.226
## F4 -0.242 0.042 -5.787 0.000 -0.312 -0.312
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw1 0.257 0.011 23.171 0.000 0.257 0.316
## .rridbw1 0.303 0.013 23.861 0.000 0.303 0.416
## .rridcw1 0.072 0.006 12.636 0.000 0.072 0.094
## .rriddw1 0.100 0.006 15.635 0.000 0.100 0.120
## .rridew1 0.383 0.016 24.609 0.000 0.383 0.543
## .rridfw1 0.225 0.012 19.204 0.000 0.225 0.213
## .rridgw1 0.183 0.010 17.762 0.000 0.183 0.185
## .rridhw1 0.264 0.013 20.876 0.000 0.264 0.255
## .rridiw1 0.579 0.024 24.634 0.000 0.579 0.536
## .rridjw1 0.213 0.011 18.979 0.000 0.213 0.209
## .rridkw1 0.358 0.016 22.019 0.000 0.358 0.384
## .rridlw1 0.788 0.032 24.957 0.000 0.788 0.872
## .rridmw1 0.227 0.009 24.521 0.000 0.227 0.694
## .rridnw1 0.865 0.035 25.058 0.000 0.865 0.919
## .rridow1 0.444 0.020 22.087 0.000 0.444 0.463
## .rridpw1 0.570 0.026 22.028 0.000 0.570 0.479
## .rridqw1 0.244 0.015 16.040 0.000 0.244 0.270
## .rridrw1 0.310 0.014 21.880 0.000 0.310 0.383
## .rridsw1 0.214 0.009 23.313 0.000 0.214 0.423
## .rridtw1 0.111 0.005 21.220 0.000 0.111 0.280
## .rriduw1 0.058 0.004 15.823 0.000 0.058 0.145
## .rridvw1 0.078 0.004 18.065 0.000 0.078 0.181
## .rridww1 0.297 0.013 23.460 0.000 0.297 0.443
## F1 1.109 0.072 15.323 0.000 1.000 1.000
## F2 0.979 0.065 15.066 0.000 1.000 1.000
## F3 1.003 0.071 14.196 0.000 1.000 1.000
## F4 0.601 0.040 15.223 0.000 1.000 1.000# per-group factor correlations (Phi)
lavInspect(fit_mg_esemC, "cor.lv")## $`1`
## F1 F2 F3 F4
## F1 1.000
## F2 0.326 1.000
## F3 0.723 0.245 1.000
## F4 0.217 0.458 0.111 1.000
##
## $`2`
## F1 F2 F3 F4
## F1 1.000
## F2 0.400 1.000
## F3 0.689 0.333 1.000
## F4 0.250 0.371 0.171 1.000# tidy loadings by group (standardized)
ld <- standardizedSolution(fit_mg_esemC)
ld <- subset(ld, op == "=~", select = c("group","lhs","rhs","est.std"))
names(ld) <- c("Group","Factor","Variable","Loading")
ld <- ld[order(ld$Group, ld$Factor, -abs(ld$Loading)), ]
print(ld, row.names = FALSE)## Group Factor Variable Loading
## 1 F1 rridcw1 1.008
## 1 F1 rriddw1 0.990
## 1 F1 rridaw1 0.819
## 1 F1 rridbw1 0.771
## 1 F1 rridew1 0.646
## 1 F1 rridrw1 0.285
## 1 F1 rridiw1 0.273
## 1 F1 rridkw1 0.202
## 1 F1 rridmw1 -0.192
## 1 F1 rridnw1 -0.128
## 1 F1 rridfw1 0.101
## 1 F1 rridww1 0.071
## 1 F1 rridhw1 0.064
## 1 F1 rridow1 0.047
## 1 F1 rridtw1 -0.037
## 1 F1 rridsw1 -0.035
## 1 F1 rridjw1 0.033
## 1 F1 rridpw1 -0.022
## 1 F1 rriduw1 -0.021
## 1 F1 rridlw1 -0.002
## 1 F1 rridvw1 0.000
## 1 F1 rridgw1 0.000
## 1 F1 rridqw1 0.000
## 1 F2 rridgw1 0.914
## 1 F2 rridjw1 0.886
## 1 F2 rridfw1 0.877
## 1 F2 rridhw1 0.843
## 1 F2 rridiw1 0.436
## 1 F2 rridkw1 0.115
## 1 F2 rridmw1 0.109
## 1 F2 rridpw1 0.084
## 1 F2 rridrw1 -0.073
## 1 F2 rridlw1 0.064
## 1 F2 rridcw1 -0.060
## 1 F2 rriddw1 -0.047
## 1 F2 rridow1 0.040
## 1 F2 rridww1 -0.030
## 1 F2 rridnw1 -0.026
## 1 F2 rriduw1 -0.011
## 1 F2 rridew1 0.004
## 1 F2 rridtw1 -0.004
## 1 F2 rridsw1 -0.003
## 1 F2 rridbw1 0.002
## 1 F2 rridaw1 0.000
## 1 F2 rridvw1 0.000
## 1 F2 rridqw1 0.000
## 1 F3 rridqw1 0.896
## 1 F3 rridpw1 0.676
## 1 F3 rridow1 0.652
## 1 F3 rridrw1 0.636
## 1 F3 rridkw1 0.601
## 1 F3 rridlw1 0.324
## 1 F3 rridnw1 0.272
## 1 F3 rridmw1 0.196
## 1 F3 rridiw1 0.088
## 1 F3 rridcw1 -0.069
## 1 F3 rridew1 0.062
## 1 F3 rridhw1 0.057
## 1 F3 rridsw1 0.054
## 1 F3 rridbw1 -0.049
## 1 F3 rriddw1 -0.043
## 1 F3 rridww1 -0.036
## 1 F3 rridfw1 -0.033
## 1 F3 rridtw1 -0.019
## 1 F3 rriduw1 0.014
## 1 F3 rridjw1 0.006
## 1 F3 rridvw1 0.000
## 1 F3 rridgw1 0.000
## 1 F3 rridaw1 0.000
## 1 F4 rridvw1 0.930
## 1 F4 rriduw1 0.925
## 1 F4 rridww1 0.871
## 1 F4 rridtw1 0.862
## 1 F4 rridsw1 0.792
## 1 F4 rridmw1 0.435
## 1 F4 rridnw1 0.215
## 1 F4 rridew1 -0.127
## 1 F4 rridkw1 -0.127
## 1 F4 rridlw1 -0.115
## 1 F4 rridpw1 0.105
## 1 F4 rridcw1 -0.099
## 1 F4 rridrw1 -0.097
## 1 F4 rriddw1 -0.095
## 1 F4 rridow1 0.090
## 1 F4 rridiw1 0.068
## 1 F4 rridfw1 -0.056
## 1 F4 rridjw1 -0.031
## 1 F4 rridhw1 -0.031
## 1 F4 rridbw1 0.014
## 1 F4 rridaw1 0.000
## 1 F4 rridqw1 0.000
## 1 F4 rridgw1 0.000
## 2 F1 rridcw1 1.033
## 2 F1 rriddw1 0.997
## 2 F1 rridaw1 0.827
## 2 F1 rridbw1 0.793
## 2 F1 rridew1 0.650
## 2 F1 rridiw1 0.288
## 2 F1 rridrw1 0.282
## 2 F1 rridmw1 -0.277
## 2 F1 rridkw1 0.210
## 2 F1 rridnw1 -0.146
## 2 F1 rridfw1 0.106
## 2 F1 rridww1 0.083
## 2 F1 rridhw1 0.066
## 2 F1 rridow1 0.051
## 2 F1 rridtw1 -0.050
## 2 F1 rridsw1 -0.046
## 2 F1 rridjw1 0.034
## 2 F1 rriduw1 -0.029
## 2 F1 rridpw1 -0.024
## 2 F1 rridlw1 -0.003
## 2 F1 rridgw1 0.000
## 2 F1 rridvw1 0.000
## 2 F1 rridqw1 0.000
## 2 F2 rridgw1 0.903
## 2 F2 rridjw1 0.881
## 2 F2 rridfw1 0.867
## 2 F2 rridhw1 0.820
## 2 F2 rridiw1 0.430
## 2 F2 rridmw1 0.148
## 2 F2 rridkw1 0.112
## 2 F2 rridpw1 0.084
## 2 F2 rridrw1 -0.068
## 2 F2 rridlw1 0.067
## 2 F2 rridcw1 -0.058
## 2 F2 rriddw1 -0.045
## 2 F2 rridow1 0.040
## 2 F2 rridww1 -0.032
## 2 F2 rridnw1 -0.028
## 2 F2 rriduw1 -0.014
## 2 F2 rridtw1 -0.005
## 2 F2 rridsw1 -0.004
## 2 F2 rridew1 0.004
## 2 F2 rridbw1 0.002
## 2 F2 rridqw1 0.000
## 2 F2 rridvw1 0.000
## 2 F2 rridaw1 0.000
## 2 F3 rridqw1 0.854
## 2 F3 rridpw1 0.685
## 2 F3 rridow1 0.665
## 2 F3 rridrw1 0.598
## 2 F3 rridkw1 0.594
## 2 F3 rridlw1 0.342
## 2 F3 rridnw1 0.294
## 2 F3 rridmw1 0.269
## 2 F3 rridiw1 0.088
## 2 F3 rridsw1 0.067
## 2 F3 rridcw1 -0.067
## 2 F3 rridew1 0.059
## 2 F3 rridhw1 0.056
## 2 F3 rridbw1 -0.048
## 2 F3 rriddw1 -0.042
## 2 F3 rridww1 -0.040
## 2 F3 rridfw1 -0.033
## 2 F3 rridtw1 -0.024
## 2 F3 rriduw1 0.018
## 2 F3 rridjw1 0.006
## 2 F3 rridgw1 0.000
## 2 F3 rridvw1 0.000
## 2 F3 rridaw1 0.000
## 2 F4 rriduw1 0.934
## 2 F4 rridvw1 0.905
## 2 F4 rridtw1 0.864
## 2 F4 rridsw1 0.759
## 2 F4 rridww1 0.742
## 2 F4 rridmw1 0.463
## 2 F4 rridnw1 0.181
## 2 F4 rridkw1 -0.097
## 2 F4 rridew1 -0.094
## 2 F4 rridlw1 -0.094
## 2 F4 rridpw1 0.083
## 2 F4 rridcw1 -0.074
## 2 F4 rridow1 0.071
## 2 F4 rriddw1 -0.071
## 2 F4 rridrw1 -0.070
## 2 F4 rridiw1 0.053
## 2 F4 rridfw1 -0.043
## 2 F4 rridjw1 -0.024
## 2 F4 rridhw1 -0.024
## 2 F4 rridbw1 0.011
## 2 F4 rridqw1 0.000
## 2 F4 rridgw1 0.000
## 2 F4 rridaw1 0.000library(dplyr)
library(tidyr)
# Group labels (e.g., "male","female")
gl <- lavInspect(fit_mg_esemC, "group.label")
# Long loadings by group
ld <- standardizedSolution(fit_mg_esemC) %>%
filter(op == "=~") %>%
transmute(
Variable = rhs,
Factor = factor(lhs, levels = paste0("F", 1:4)), # F1..F4 in order
Group = factor(group, levels = seq_along(gl), labels = gl),
Loading = est.std
)
# Wide table: columns like F1_male, F1_female, F2_male, ...
loadings_wide_all <- ld %>%
arrange(match(Variable, myvars), Group, Factor) %>%
pivot_wider(
names_from = c(Factor, Group),
values_from = Loading,
names_glue = "{Factor}_{Group}"
) %>%
arrange(match(Variable, myvars))
# Round for display (keep the numeric version if you’ll compute further)
loadings_wide_print <- loadings_wide_all %>%
mutate(across(-Variable, ~ round(.x, 2)))
# Print
print(loadings_wide_print, row.names = FALSE)## # A tibble: 23 × 9
## Variable F1_1 F2_1 F3_1 F4_1 F1_2 F2_2 F3_2 F4_2
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 rridaw1 0.82 0 0 0 0.83 0 0 0
## 2 rridbw1 0.77 0 -0.05 0.01 0.79 0 -0.05 0.01
## 3 rridcw1 1.01 -0.06 -0.07 -0.1 1.03 -0.06 -0.07 -0.07
## 4 rriddw1 0.99 -0.05 -0.04 -0.1 1 -0.04 -0.04 -0.07
## 5 rridew1 0.65 0 0.06 -0.13 0.65 0 0.06 -0.09
## 6 rridfw1 0.1 0.88 -0.03 -0.06 0.11 0.87 -0.03 -0.04
## 7 rridgw1 0 0.91 0 0 0 0.9 0 0
## 8 rridhw1 0.06 0.84 0.06 -0.03 0.07 0.82 0.06 -0.02
## 9 rridiw1 0.27 0.44 0.09 0.07 0.29 0.43 0.09 0.05
## 10 rridjw1 0.03 0.89 0.01 -0.03 0.03 0.88 0.01 -0.02
## # ℹ 13 more rows# install.packages(c("dplyr","tidyr","knitr","kableExtra")) # if needed
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
# Label groups explicitly: 1 = Male, 2 = Female
grp_labels <- c("1" = "Male", "2" = "Female")
factors <- paste0("F", 1:4)
# Long loadings with explicit group labels
ld <- standardizedSolution(fit_mg_esemC) %>%
filter(op == "=~") %>%
transmute(
Variable = rhs,
Factor = factor(lhs, levels = factors),
Group = recode(as.character(group), !!!grp_labels),
Loading = est.std
)
# Wide per group
male <- ld %>% filter(Group == "Male") %>%
select(Variable, Factor, Loading) %>%
pivot_wider(names_from = Factor, values_from = Loading) %>%
arrange(match(Variable, myvars))
female <- ld %>% filter(Group == "Female") %>%
select(Variable, Factor, Loading) %>%
pivot_wider(names_from = Factor, values_from = Loading) %>%
arrange(match(Variable, myvars))
# Suffix columns and join
names(male)[-1] <- paste0(names(male)[-1], "_Male")
names(female)[-1] <- paste0(names(female)[-1], "_Female")
loadings_apa <- left_join(male, female, by = "Variable")
# Order columns: Variable, F1..F4 (Male), F1..F4 (Female)
loadings_apa <- loadings_apa %>%
select(Variable,
F1_Male, F2_Male, F3_Male, F4_Male,
F1_Female, F2_Female, F3_Female, F4_Female)
# Format: two decimals; bold |loading| >= .30 (APA-ish)
fmt <- if (knitr::is_latex_output()) "latex" else "html"
thresh <- 0.30
loadings_apa_fmt <- loadings_apa %>%
mutate(across(-Variable, ~ {
x <- .
lbl <- sprintf("%.2f", x)
ifelse(abs(x) >= thresh,
kableExtra::cell_spec(lbl, format = fmt, bold = TRUE),
lbl)
}))
# Render APA-style table
# Notice these will not match across groups because factor variances vary
# across groups.
kbl(loadings_apa_fmt,
booktabs = TRUE,
align = c("l", rep("r", 8)),
caption = "Sup. Table 16.1\nStandardized Factor Loadings by Group (will not match across groups-factor variances vary)",
escape = FALSE) %>%
add_header_above(c(" " = 1, "Male (Group 1)" = 4, "Female (Group 2)" = 4)) %>%
kable_classic(full_width = FALSE, html_font = "Times New Roman") %>%
footnote(
general = "Note. Values are standardized loadings (est.std). Bold indicates |loading| ≥ .30.",
threeparttable = TRUE
)| Variable | F1_Male | F2_Male | F3_Male | F4_Male | F1_Female | F2_Female | F3_Female | F4_Female |
|---|---|---|---|---|---|---|---|---|
| rridaw1 | 0.82 | 0.00 | -0.00 | 0.00 | 0.83 | 0.00 | -0.00 | 0.00 |
| rridbw1 | 0.77 | 0.00 | -0.05 | 0.01 | 0.79 | 0.00 | -0.05 | 0.01 |
| rridcw1 | 1.01 | -0.06 | -0.07 | -0.10 | 1.03 | -0.06 | -0.07 | -0.07 |
| rriddw1 | 0.99 | -0.05 | -0.04 | -0.10 | 1.00 | -0.04 | -0.04 | -0.07 |
| rridew1 | 0.65 | 0.00 | 0.06 | -0.13 | 0.65 | 0.00 | 0.06 | -0.09 |
| rridfw1 | 0.10 | 0.88 | -0.03 | -0.06 | 0.11 | 0.87 | -0.03 | -0.04 |
| rridgw1 | 0.00 | 0.91 | -0.00 | -0.00 | 0.00 | 0.90 | 0.00 | 0.00 |
| rridhw1 | 0.06 | 0.84 | 0.06 | -0.03 | 0.07 | 0.82 | 0.06 | -0.02 |
| rridiw1 | 0.27 | 0.44 | 0.09 | 0.07 | 0.29 | 0.43 | 0.09 | 0.05 |
| rridjw1 | 0.03 | 0.89 | 0.01 | -0.03 | 0.03 | 0.88 | 0.01 | -0.02 |
| rridkw1 | 0.20 | 0.12 | 0.60 | -0.13 | 0.21 | 0.11 | 0.59 | -0.10 |
| rridlw1 | -0.00 | 0.06 | 0.32 | -0.11 | -0.00 | 0.07 | 0.34 | -0.09 |
| rridmw1 | -0.19 | 0.11 | 0.20 | 0.44 | -0.28 | 0.15 | 0.27 | 0.46 |
| rridnw1 | -0.13 | -0.03 | 0.27 | 0.22 | -0.15 | -0.03 | 0.29 | 0.18 |
| rridow1 | 0.05 | 0.04 | 0.65 | 0.09 | 0.05 | 0.04 | 0.66 | 0.07 |
| rridpw1 | -0.02 | 0.08 | 0.68 | 0.11 | -0.02 | 0.08 | 0.68 | 0.08 |
| rridqw1 | -0.00 | -0.00 | 0.90 | -0.00 | -0.00 | -0.00 | 0.85 | -0.00 |
| rridrw1 | 0.29 | -0.07 | 0.64 | -0.10 | 0.28 | -0.07 | 0.60 | -0.07 |
| rridsw1 | -0.04 | -0.00 | 0.05 | 0.79 | -0.05 | -0.00 | 0.07 | 0.76 |
| rridtw1 | -0.04 | -0.00 | -0.02 | 0.86 | -0.05 | -0.00 | -0.02 | 0.86 |
| rriduw1 | -0.02 | -0.01 | 0.01 | 0.92 | -0.03 | -0.01 | 0.02 | 0.93 |
| rridvw1 | 0.00 | -0.00 | -0.00 | 0.93 | -0.00 | -0.00 | -0.00 | 0.91 |
| rridww1 | 0.07 | -0.03 | -0.04 | 0.87 | 0.08 | -0.03 | -0.04 | 0.74 |
| Note: | ||||||||
| Note. Values are standardized loadings (est.std). Bold indicates |loading| ≥ .30. |
# install.packages(c("dplyr","tidyr","knitr","kableExtra")) # if needed
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
# -------------------------------
# 0) Setup
# -------------------------------
grp_labels <- c("1" = "Male", "2" = "Female") # map lavaan group indices -> labels
factors <- paste0("F", 1:4)
# Unstandardized parameter table from lavaan
pe <- parameterEstimates(fit_mg_esemC, standardized = FALSE)
# -------------------------------
# 1) Unstandardized loadings by group (APA-style table)
# -------------------------------
ld_u <- pe %>%
filter(op == "=~", lhs %in% factors, rhs %in% myvars) %>%
transmute(
Variable = rhs,
Factor = factor(lhs, levels = factors),
Group = recode(as.character(group), !!!grp_labels),
Loading = est
)
# Build wide tables per group then join
male_u <- ld_u %>%
filter(Group == "Male") %>%
select(Variable, Factor, Loading) %>%
pivot_wider(names_from = Factor, values_from = Loading) %>%
arrange(match(Variable, myvars))
female_u <- ld_u %>%
filter(Group == "Female") %>%
select(Variable, Factor, Loading) %>%
pivot_wider(names_from = Factor, values_from = Loading) %>%
arrange(match(Variable, myvars))
names(male_u)[-1] <- paste0(names(male_u)[-1], "_Male")
names(female_u)[-1] <- paste0(names(female_u)[-1], "_Female")
loadings_unstd <- left_join(male_u, female_u, by = "Variable") %>%
select(Variable,
F1_Male, F2_Male, F3_Male, F4_Male,
F1_Female, F2_Female, F3_Female, F4_Female)
# Format to two decimals (APA-ish)
loadings_unstd_fmt <- loadings_unstd %>%
mutate(across(-Variable, ~ sprintf("%.2f", .)))
# Render APA table
kbl(loadings_unstd_fmt,
booktabs = TRUE,
align = c("l", rep("r", 8)),
caption = "Table 16.10\nUnstandardized Factor Loadings by Group (4-Factor Target-Rotated ESEM)",
escape = FALSE) %>%
add_header_above(c(" " = 1, "Male (Group 1)" = 4, "Female (Group 2)" = 4)) %>%
kable_classic(full_width = FALSE, html_font = "Times New Roman") %>%
footnote(
general = "Note. Entries are unstandardized loadings (est).",
threeparttable = TRUE
)| Variable | F1_Male | F2_Male | F3_Male | F4_Male | F1_Female | F2_Female | F3_Female | F4_Female |
|---|---|---|---|---|---|---|---|---|
| rridaw1 | 0.71 | 0.00 | -0.00 | 0.00 | 0.71 | 0.00 | -0.00 | 0.00 |
| rridbw1 | 0.64 | 0.00 | -0.04 | 0.01 | 0.64 | 0.00 | -0.04 | 0.01 |
| rridcw1 | 0.86 | -0.05 | -0.06 | -0.08 | 0.86 | -0.05 | -0.06 | -0.08 |
| rriddw1 | 0.87 | -0.04 | -0.04 | -0.08 | 0.87 | -0.04 | -0.04 | -0.08 |
| rridew1 | 0.52 | 0.00 | 0.05 | -0.10 | 0.52 | 0.00 | 0.05 | -0.10 |
| rridfw1 | 0.10 | 0.90 | -0.03 | -0.06 | 0.10 | 0.90 | -0.03 | -0.06 |
| rridgw1 | 0.00 | 0.91 | -0.00 | -0.00 | 0.00 | 0.91 | 0.00 | 0.00 |
| rridhw1 | 0.06 | 0.84 | 0.06 | -0.03 | 0.06 | 0.84 | 0.06 | -0.03 |
| rridiw1 | 0.28 | 0.45 | 0.09 | 0.07 | 0.28 | 0.45 | 0.09 | 0.07 |
| rridjw1 | 0.03 | 0.90 | 0.01 | -0.03 | 0.03 | 0.90 | 0.01 | -0.03 |
| rridkw1 | 0.19 | 0.11 | 0.57 | -0.12 | 0.19 | 0.11 | 0.57 | -0.12 |
| rridlw1 | -0.00 | 0.06 | 0.32 | -0.12 | -0.00 | 0.06 | 0.32 | -0.12 |
| rridmw1 | -0.15 | 0.09 | 0.15 | 0.34 | -0.15 | 0.09 | 0.15 | 0.34 |
| rridnw1 | -0.13 | -0.03 | 0.28 | 0.23 | -0.13 | -0.03 | 0.28 | 0.23 |
| rridow1 | 0.05 | 0.04 | 0.65 | 0.09 | 0.05 | 0.04 | 0.65 | 0.09 |
| rridpw1 | -0.02 | 0.09 | 0.75 | 0.12 | -0.02 | 0.09 | 0.75 | 0.12 |
| rridqw1 | -0.00 | -0.00 | 0.81 | -0.00 | -0.00 | -0.00 | 0.81 | -0.00 |
| rridrw1 | 0.24 | -0.06 | 0.54 | -0.08 | 0.24 | -0.06 | 0.54 | -0.08 |
| rridsw1 | -0.03 | -0.00 | 0.05 | 0.70 | -0.03 | -0.00 | 0.05 | 0.70 |
| rridtw1 | -0.03 | -0.00 | -0.02 | 0.70 | -0.03 | -0.00 | -0.02 | 0.70 |
| rriduw1 | -0.02 | -0.01 | 0.01 | 0.76 | -0.02 | -0.01 | 0.01 | 0.76 |
| rridvw1 | 0.00 | -0.00 | -0.00 | 0.77 | -0.00 | -0.00 | -0.00 | 0.77 |
| rridww1 | 0.06 | -0.03 | -0.03 | 0.78 | 0.06 | -0.03 | -0.03 | 0.78 |
| Note: | ||||||||
| Note. Entries are unstandardized loadings (est). |
# -------------------------------
# 2) Factor variances and factor means (compact APA table)
# -------------------------------
# Factor variances (lhs == rhs, op '~~')
var_u <- pe %>%
filter(op == "~~", lhs == rhs, lhs %in% factors) %>%
transmute(Factor = lhs,
Group = recode(as.character(group), !!!grp_labels),
Variance = est) %>%
distinct() %>%
pivot_wider(names_from = Group, values_from = Variance) %>%
arrange(Factor)
# Factor means (op '~1')
mean_u <- pe %>%
filter(op == "~1", lhs %in% factors) %>%
transmute(Factor = lhs,
Group = recode(as.character(group), !!!grp_labels),
Mean = est) %>%
distinct() %>%
pivot_wider(names_from = Group, values_from = Mean) %>%
arrange(Factor)
# Merge variances + means into one APA-friendly table
vm <- var_u %>%
rename(Var_Male = Male, Var_Female = Female) %>%
left_join(mean_u %>% rename(Mean_Male = Male, Mean_Female = Female),
by = "Factor")
vm_fmt <- vm %>%
mutate(across(-Factor, ~ sprintf("%.2f", .)))
kbl(vm_fmt,
booktabs = TRUE,
align = c("l","r","r","r","r"),
caption = "Table 16.10 Part II \nLatent Factor Variances and Means by Group",
col.names = c("Factor", "Variance (Male)", "Variance (Female)",
"Mean (Male)", "Mean (Female)")) %>%
kable_classic(full_width = FALSE, html_font = "Times New Roman") %>%
footnote(
general = "Note. Entries are unstandardized estimates (est). Group 1 is Male; Group 2 is Female.",
threeparttable = TRUE
)| Factor | Variance (Male) | Variance (Female) | Mean (Male) | Mean (Female) |
|---|---|---|---|---|
| F1 | 1.00 | 1.11 | 0.00 | -0.13 |
| F2 | 1.00 | 0.98 | 0.00 | -0.05 |
| F3 | 1.00 | 1.00 | 0.00 | -0.23 |
| F4 | 1.00 | 0.60 | 0.00 | -0.24 |
| Note: | ||||
| Note. Entries are unstandardized estimates (est). Group 1 is Male; Group 2 is Female. |
Multi-Group Invariance with ESEM Models A-G (Table 16.9)
(These will not match values reported in the book at all because the data are simulated and the raw data contained missing data.)
The code, however, shows how you would set this up.
It is not easy to fit Model E using lavaan’s defaults and, given that it’s unlikely for this type of data it was omitted.
# install.packages("lavaan") # if needed
library(lavaan)
# --- data ---
RFDItems <- read.table(
"RFDItems.txt",
sep="\t", header=TRUE)
RFDItems$sex <- factor(RFDItems$sex)
myvars <- c(
"rridaw1","rridbw1","rridcw1","rriddw1","rridew1","rridfw1","rridgw1",
"rridhw1","rridiw1","rridjw1","rridkw1","rridlw1","rridmw1","rridnw1",
"rridow1","rridpw1","rridqw1","rridrw1","rridsw1","rridtw1","rriduw1",
"rridvw1","rridww1"
)
# --- ESEM model with a 4-factor EFA block in each group ---
efa_block <- paste(
sprintf('efa("efa4")*F%d =~ %s', 1, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 2, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 3, paste(myvars, collapse = " + ")),
sprintf('efa("efa4")*F%d =~ %s', 4, paste(myvars, collapse = " + ")),
sep = "\n"
)
# --- Minimal shared target to anchor factor orientation across groups ---
# (Zeros on non-target columns for four anchor items; all else NA = don't care)
T_shared <- matrix(NA_real_, nrow = length(myvars), ncol = 4,
dimnames = list(myvars, paste0("F", 1:4)))
T_shared["rridaw1", c("F2","F3","F4")] <- 0 # anchor for F1
T_shared["rridgw1", c("F1","F3","F4")] <- 0 # anchor for F2
T_shared["rridqw1", c("F1","F2","F4")] <- 0 # anchor for F3
T_shared["rridvw1", c("F1","F2","F3")] <- 0 # anchor for F4
# --- Fit: MG-ESEM with loadings+intercepts invariant; other things variant ---
fit_mg_esemA <- cfa(
model = efa_block,
data = RFDItems,
group = "sex",
meanstructure = TRUE, # needed for intercept/mean invariance
estimator = "ml",
# rotation: target (oblique), same target in both groups
rotation = "target",
rotation.args = list(target = T_shared, orthogonal = FALSE),
)
fit_mg_esemB <- cfa(
model = efa_block,
data = RFDItems,
group = "sex",
meanstructure = TRUE, # needed for intercept/mean invariance
estimator = "ml",
# rotation: target (oblique), same target in both groups
rotation = "target",
rotation.args = list(target = T_shared, orthogonal = FALSE),
# equality constraints ACROSS groups:
group.equal = c("loadings")
# (Note: we are NOT constraining lv.variances, lv.covariances, or residuals,
# so factor means/variances/covariances and error variances are free by group.)
)
fit_mg_esemC <- cfa(
model = efa_block,
data = RFDItems,
group = "sex",
meanstructure = TRUE,
estimator = "ml",
rotation = "target",
rotation.args = list(target = T_shared, orthogonal = FALSE),
group.equal = c("loadings", "intercepts")
)
fit_mg_esemD <- cfa(
model = efa_block,
data = RFDItems,
group = "sex",
meanstructure = TRUE,
estimator = "ml",
rotation = "target",
rotation.args = list(target = T_shared, orthogonal = FALSE),
group.equal = c("loadings", "intercepts", "lv.covariances", "lv.variances")
)
# F: loadings + intercepts + residuals + lv (co)variances equal across groups
fit_mg_esemF <- cfa(
model = efa_block,
data = RFDItems,
group = "sex",
meanstructure = TRUE,
estimator = "ml",
rotation = "target",
rotation.args = list(target = T_shared, orthogonal = FALSE),
group.equal = c("loadings", "intercepts", "lv.covariances", "residuals")
)
# G: same as F, and ALSO constrain LATENT means across groups
fit_mg_esemG <- cfa(
model = efa_block,
data = RFDItems,
group = "sex",
meanstructure = TRUE,
estimator = "ml",
rotation = "target",
rotation.args = list(target = T_shared, orthogonal = FALSE),
group.equal = c("loadings", "intercepts", "lv.covariances", "lv.variances",
"residuals", "means") # <-- use lv.means (not "means")
)
# install.packages(c("dplyr","knitr","kableExtra")) # if needed
library(dplyr)
library(knitr)
library(kableExtra)
library(lavaan)
# Helper: extract fit measures and format for APA-style table
make_fit_table <- function(named_fits) {
rows <- lapply(names(named_fits), function(nm) {
fit <- named_fits[[nm]]
ok <- inherits(fit, "lavaan") && lavInspect(fit, "converged")
if (!ok) {
return(data.frame(Model = nm,
`χ² (df)` = "—",
BIC = "—",
TLI = "—",
RMSEA = "—",
AIC = "—",
stringsAsFactors = FALSE))
}
fm <- fitMeasures(fit, c("chisq","df","aic","bic","tli",
"rmsea","rmsea.ci.lower","rmsea.ci.upper"))
data.frame(
Model = nm,
`χ² (df)` = sprintf("%.2f (%d)", fm["chisq"], as.integer(fm["df"])),
BIC = formatC(fm["bic"], digits = 1, format = "f"),
TLI = sprintf("%.3f", fm["tli"]),
# RMSEA with CI underneath (HTML line break so it prints nicely)
RMSEA = sprintf("%.3f<br/><span style='font-size:85%%'>(%.3f, %.3f)</span>",
fm["rmsea"], fm["rmsea.ci.lower"], fm["rmsea.ci.upper"]),
AIC = formatC(fm["aic"], digits = 1, format = "f"),
stringsAsFactors = FALSE
)
})
do.call(rbind, rows)
}
# Build the table from your fits
tbl_fit <- make_fit_table(list(
"A (Configural)" = fit_mg_esemA,
"B (Loadings equal)" = fit_mg_esemB,
"C (Loadings + Intercepts equal)" = fit_mg_esemC,
"D (+ LV vars & covs equal)" = fit_mg_esemD,
"F (+ Residuals equal)" = fit_mg_esemF,
"G (+ LV means equal)" = fit_mg_esemG
))
# Render APA-style table (HTML/Word)
kbl(tbl_fit,
booktabs = TRUE,
align = c("l","r","r","r","r","r"),
caption = "Table 16.9\nModel Fit Indices for Multi-Group ESEM",
escape = FALSE # allow the RMSEA line break/HTML span
) %>%
kable_classic(full_width = FALSE, html_font = "Times New Roman") %>%
footnote(
general = "Note. RMSEA values include the 90% confidence interval on the line below.",
threeparttable = TRUE
)| Model | χ…df. | BIC | TLI | RMSEA | AIC | |
|---|---|---|---|---|---|---|
| bic | A (Configural) | 2444.32 (334) | 96276.1 | 0.913 |
0.077 (0.074, 0.079) |
94777.3 |
| bic1 | B (Loadings equal) | 2598.84 (410) | 95847.1 | 0.926 |
0.070 (0.068, 0.073) |
94779.8 |
| bic2 | C (Loadings + Intercepts equal) | 2812.12 (429) | 95914.5 | 0.923 |
0.072 (0.069, 0.074) |
94955.0 |
| bic3 | D (+ LV vars & covs equal) | 2890.75 (439) | 95916.4 | 0.923 |
0.072 (0.069, 0.074) |
95013.7 |
| bic4 | F (+ Residuals equal) | 3171.37 (458) | 96051.1 | 0.918 |
0.074 (0.072, 0.077) |
95256.3 |
| bic5 | G (+ LV means equal) | 3307.14 (466) | 96125.5 | 0.916 |
0.075 (0.073, 0.078) |
95376.1 |
| Note: | ||||||
| Note. RMSEA values include the 90% confidence interval on the line below. |