## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Warning: package 'lavaan' was built under R version 4.5.3
## This is lavaan 0.6-21
## lavaan is FREE software! Please report any bugs.
## Warning: package 'semTools' was built under R version 4.5.3
##
## ###############################################################################
## This is semTools 0.5-8
## All users of R (or SEM) are invited to submit functions or ideas for functions.
## ###############################################################################
##
## Attaching package: 'semTools'
##
## The following object is masked from 'package:readr':
##
## clipboard
## Warning: package 'MVN' was built under R version 4.5.3
##
## Attaching package: 'psych'
##
## The following object is masked from 'package:MVN':
##
## mardia
##
## The following objects are masked from 'package:semTools':
##
## reliability, skew
##
## The following object is masked from 'package:lavaan':
##
## cor2cov
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:psych':
##
## logit
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
## Warning: package 'semPlot' was built under R version 4.5.3
df_sem <- read.csv("Data POS-ITS 2.csv")[, -(1:8)]
ordinal_items <- c(
paste0("POS", 1:8),
paste0("JS", 1:5),
paste0("WLB", 1:3),
paste0("ITS", 1:5)
)
head(df_sem)read.csv("Data POS-ITS 2.csv")[, (1:8)] |>
rename(
`Position` = Q1,
`Tenure (Years)` = Q2,
`Companies Worked For` = Q3,
`Age` = Q4,
`Gender` = Q5,
`Education` = Q6,
`Monthly Income` = Q7,
`Marital Status` = Q8
) |>
mutate(
Position = factor(Position, levels = 1:5,
labels = c("Worker", "Sub-leader", "Leader", "Exec. Staff", "Other")),
`Tenure (Years)` = factor(`Tenure (Years)`, levels = 1:5,
labels = c("<2 yrs", "2-4 yrs", "5-7 yrs", "8-10 yrs", ">10 yrs")),
`Companies Worked For` = factor(`Companies Worked For`, levels = 1:5,
labels = c("<2", "2-4", "5-7", "8-10", ">10")),
Age = factor(Age, levels = 1:4,
labels = c("18-27", "28-37", "37-46", ">47")),
Gender = factor(Gender, levels = 0:1,
labels = c("Male", "Female")),
Education = factor(Education, levels = 1:3,
labels = c("High School", "College/Uni", "Master+")),
`Monthly Income` = factor(`Monthly Income`, levels = 1:4,
labels = c("<$250", "$250-400", "$401-600", ">$600")),
`Marital Status` = factor(`Marital Status`, levels = 0:1,
labels = c("Single", "Married"))
) |>
pivot_longer(cols = everything(), names_to = "Characteristic", values_to = "Category") |>
ggplot(aes(y = Category)) +
geom_bar(fill = "#3498db", color = "black", alpha = 0.8) +
facet_wrap(~ Characteristic, scales = "free", ncol = 2) +
theme_minimal() +
labs(
title = "Respondent Characteristics Profile",
x = "Number of Respondents",
y = NULL
) +
theme(
strip.text = element_text(face = "bold", size = 10),
plot.title = element_text(face = "bold", size = 14, hjust = 0.5)
)Although the material not include the “no outliers” assumption, SEM relies on covariance matrix which could be distort by significant multivariate outliers. Therefore, dropping outliers using mahalanobis with Chi-Squared cut-off is needed.
## Before removal: 604
m_dist <- mahalanobis(df_sem,
colMeans(df_sem),
cov(df_sem))
cutoff <- qchisq(0.95, df = ncol(df_sem))
keep <- m_dist < cutoff
df_sem <- df_sem[keep, ]
cat("After removal:", nrow(df_sem), "\n")## After removal: 552
data.frame(
Index = 1:length(m_dist),
Distance = m_dist,
IsOutlier = m_dist >= cutoff
) |>
ggplot(aes(x = Index, y = Distance, color = IsOutlier)) +
geom_point(alpha = 0.5) +
geom_hline(yintercept = cutoff, linetype = "dashed", color = "red") +
scale_color_manual(values = c("FALSE" = "#2c3e50", "TRUE" = "#e74c3c"),
labels = c("Normal", "Outlier")) +
theme_minimal() +
labs(title = "Mahalanobis Distance Outlier Detection",
subtitle = paste("Cutoff =", round(cutoff, 2)),
x = "Observation Index", y = "Mahalanobis Distance")Standard method in SEM is Maximum Likelihood (ML) estimation. ML was mathematically designed assuming that your data follows a perfect multivariate normal distribution. Therefore, checking the is needed the multivariate normality is necessarily.
##
## ── Multivariate Normality Test Results ─────────────────────────────────────────
## Test Statistic p.value Method MVN
## 1 Henze-Zirkler 11.645 <0.001 asymptotic ✗ Not normal
## Test Statistic p.value Method MVN
## 1 Henze-Zirkler 11.645 <0.001 asymptotic ✗ Not normal
##
## ── Multivariate Normality Test Results ─────────────────────────────────────────
## Test Statistic p.value Method MVN
## 1 Henze-Zirkler 11.645 <0.001 asymptotic ✗ Not normal
## Test Variable Statistic p.value Normality
## 1 Anderson-Darling POS1 117.435 <0.001 ✗ Not normal
## 2 Anderson-Darling POS2 46.197 <0.001 ✗ Not normal
## 3 Anderson-Darling POS3 76.890 <0.001 ✗ Not normal
## 4 Anderson-Darling POS4 72.607 <0.001 ✗ Not normal
## 5 Anderson-Darling POS5 73.902 <0.001 ✗ Not normal
## 6 Anderson-Darling POS6 64.302 <0.001 ✗ Not normal
## 7 Anderson-Darling POS7 64.615 <0.001 ✗ Not normal
## 8 Anderson-Darling POS8 77.440 <0.001 ✗ Not normal
## 9 Anderson-Darling JS1 95.068 <0.001 ✗ Not normal
## 10 Anderson-Darling JS2 77.457 <0.001 ✗ Not normal
## 11 Anderson-Darling JS3 64.392 <0.001 ✗ Not normal
## 12 Anderson-Darling JS4 63.501 <0.001 ✗ Not normal
## 13 Anderson-Darling JS5 62.547 <0.001 ✗ Not normal
## 14 Anderson-Darling WLB1 99.647 <0.001 ✗ Not normal
## 15 Anderson-Darling WLB2 53.882 <0.001 ✗ Not normal
## 16 Anderson-Darling WLB3 66.641 <0.001 ✗ Not normal
## 17 Anderson-Darling ITS1 117.994 <0.001 ✗ Not normal
## 18 Anderson-Darling ITS2 45.268 <0.001 ✗ Not normal
## 19 Anderson-Darling ITS3 56.843 <0.001 ✗ Not normal
## 20 Anderson-Darling ITS4 63.141 <0.001 ✗ Not normal
## 21 Anderson-Darling ITS5 75.848 <0.001 ✗ Not normal
par(mfrow = c(5, 5), mar = c(2, 2, 2, 1))
for (var in names(df_sem)) {
qqnorm(df_sem[[var]], main = paste("Q-Q:", var), cex.main = 0.8)
qqline(df_sem[[var]], col = "red")
}
par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)The Henze-Zirkler test was used to assess multivariate normality across all 21 observed indicators. The test statistic was 11.645 with a p-value below 0.001, indicating that the data do not follow a multivariate normal distribution. Individual univariate normality was also assessed using the Anderson-Darling test. All 21 indicators yielded p-values below 0.001, confirming that each variable deviates significantly from a normal distribution. These results are also visually supported by the Q-Q plots above, where most data points deviate from the diagonal reference line. Because the normality assumption is violated, the analysis proceeds using the Diagonally Weighted Least Squares (DWLS) estimator with robust standard errors, which is appropriate for ordinal and non-normally distributed data.
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(df_sem))
## Overall MSA = 0.95
## MSA for each item =
## POS1 POS2 POS3 POS4 POS5 POS6 POS7 POS8 JS1 JS2 JS3 JS4 JS5 WLB1 WLB2 WLB3
## 0.94 0.93 0.95 0.95 0.95 0.95 0.95 0.96 0.94 0.97 0.95 0.94 0.94 0.94 0.95 0.90
## ITS1 ITS2 ITS3 ITS4 ITS5
## 0.96 0.96 0.96 0.97 0.95
The Kaiser-Meyer-Olkin (KMO) test was conducted to evaluate whether the correlation structure among variables is adequate for factor analysis. The overall MSA value was 0.95, which far exceeds the minimum threshold of 0.50. MSA values for individual items ranged from 0.90 to 0.97, with all items surpassing the threshold. This indicates that the data are highly suitable for Structural Equation Modeling.
vif_model <- lm(ITS1 ~ POS1 + POS2 + POS3 + POS4 + POS5 + POS6 + POS7 + POS8 +
JS1 + JS2 + JS3 + JS4 + JS5 +
WLB1 + WLB2 + WLB3 +
ITS2 + ITS3 + ITS4 + ITS5,
data = df_sem)
vif_result <- vif(vif_model)
vif_df <- data.frame(
Variable = names(vif_result),
VIF = as.numeric(vif_result)
)
# Plot the results
ggplot(vif_df, aes(x = reorder(Variable, VIF), y = VIF, fill = VIF > 5)) +
geom_bar(stat = "identity", width = 0.7) +
geom_hline(yintercept = 5, linetype = "dashed", color = "red", linewidth = 0.8) +
geom_hline(yintercept = 3.3, linetype = "dotted", color = "orange", linewidth = 0.8) +
coord_flip() + # Memutar agar nama variabel mudah dibaca
scale_fill_manual(values = c("FALSE" = "#2ecc71", "TRUE" = "#e74c3c"),
guide = "none") +
theme_minimal() +
labs(title = "Multicollinearity Check (VIF Results)",
x = "Variables",
y = "Variance Inflation Factor (VIF)")The Variance Inflation Factor (VIF) was calculated for all indicators to check for multicollinearity. As shown in the chart, all VIF values fall well below 3.3, which is the threshold recommended as a safe boundary in formative indicator evaluation. No indicator approaches the more critical threshold of 5. This confirms that multicollinearity is not a concern in this dataset, and all indicators can be retained for further analysis.
model_pos_only <- paste("POS =~", paste0("POS", 1:8, collapse = " + "))
fit_pos <- cfa(model_pos_only, data = df_sem, ordered = paste0("POS", 1:8))## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= 1.708701e-16)
## is close to zero. This may be a symptom that the model is not identified.
## lavaan 0.6-21 ended normally after 22 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 30
##
## Number of observations 552
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 40.854 60.007
## Degrees of freedom 20 20
## P-value (Unknown) NA 0.000
## Scaling correction factor 0.699
## Shift parameter 1.591
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 9811.349 5669.005
## Degrees of freedom 28 28
## P-value NA 0.000
## Scaling correction factor 1.734
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.998 0.993
## Tucker-Lewis Index (TLI) 0.997 0.990
##
## Robust Comparative Fit Index (CFI) 0.940
## Robust Tucker-Lewis Index (TLI) 0.916
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.044 0.060
## 90 Percent confidence interval - lower 0.024 0.043
## 90 Percent confidence interval - upper 0.063 0.078
## P-value H_0: RMSEA <= 0.050 0.689 0.154
## P-value H_0: RMSEA >= 0.080 0.000 0.033
##
## Robust RMSEA 0.138
## 90 Percent confidence interval - lower 0.102
## 90 Percent confidence interval - upper 0.176
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 0.995
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.043 0.043
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS =~
## POS1 1.000 0.832 0.832
## POS2 0.982 0.032 30.312 0.000 0.817 0.817
## POS3 0.943 0.032 29.398 0.000 0.785 0.785
## POS4 0.986 0.035 28.152 0.000 0.821 0.821
## POS5 0.922 0.033 27.600 0.000 0.767 0.767
## POS6 1.024 0.034 30.068 0.000 0.852 0.852
## POS7 1.002 0.030 33.538 0.000 0.834 0.834
## POS8 0.997 0.034 29.295 0.000 0.830 0.830
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS1|t1 -0.939 0.063 -14.928 0.000 -0.939 -0.939
## POS1|t2 1.752 0.097 18.067 0.000 1.752 1.752
## POS2|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS2|t2 -0.831 0.061 -13.702 0.000 -0.831 -0.831
## POS2|t3 0.538 0.056 9.556 0.000 0.538 0.538
## POS3|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS3|t2 -0.953 0.063 -15.077 0.000 -0.953 -0.953
## POS3|t3 1.107 0.067 16.486 0.000 1.107 1.107
## POS4|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## POS4|t2 -0.932 0.063 -14.853 0.000 -0.932 -0.932
## POS4|t3 1.066 0.066 16.149 0.000 1.066 1.066
## POS5|t1 -0.967 0.064 -15.224 0.000 -0.967 -0.967
## POS5|t2 1.035 0.065 15.872 0.000 1.035 1.035
## POS6|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## POS6|t2 -0.864 0.061 -14.091 0.000 -0.864 -0.864
## POS6|t3 0.997 0.064 15.516 0.000 0.997 0.997
## POS7|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## POS7|t2 -0.831 0.061 -13.702 0.000 -0.831 -0.831
## POS7|t3 1.035 0.065 15.872 0.000 1.035 1.035
## POS8|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS8|t2 -0.884 0.062 -14.322 0.000 -0.884 -0.884
## POS8|t3 1.186 0.070 17.053 0.000 1.186 1.186
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .POS1 0.307 0.307 0.307
## .POS2 0.332 0.332 0.332
## .POS3 0.383 0.383 0.383
## .POS4 0.327 0.327 0.327
## .POS5 0.411 0.411 0.411
## .POS6 0.274 0.274 0.274
## .POS7 0.305 0.305 0.305
## .POS8 0.311 0.311 0.311
## POS 0.693 0.041 16.887 0.000 1.000 1.000
## POS
## 0.669
## For constructs with categorical indicators, Zumbo et al.`s (2007) "ordinal alpha" is calculated in addition to the standard alpha, which treats ordinal variables as numeric. See Chalmers (2018) for a critique of "alpha.ord" and the response by Zumbo & Kroc (2019). Likewise, average variance extracted is calculated from polychoric (polyserial) not Pearson correlations.
## POS
## alpha 0.8841865
## alpha.ord 0.9388436
## omega 0.8901366
## omega2 0.8901366
## omega3 0.9013816
## avevar 0.6685875
A single-construct CFA was run for the POS construct using the DWLS estimator with robust standard errors. The scaled model fit indices were CFI = 0.993 and TLI = 0.990, both meeting the threshold of >= 0.90. The scaled RMSEA was 0.060, which is below the threshold of 0.08, and the SRMR was 0.043, which is below the threshold of 0.05. The model fit is therefore considered acceptable.
All eight indicators showed standardized factor loadings above 0.70, ranging from 0.767 (POS5) to 0.852 (POS6), meeting the ideal loading threshold. The AVE was 0.669, exceeding the minimum threshold of 0.50, indicating adequate convergent validity. For reliability, the ordinal alpha was 0.939 and the composite reliability (omega) was 0.890, both well above the 0.70 threshold. The POS construct therefore demonstrates satisfactory validity and reliability.
model_js_only <- paste("JS =~", paste0("JS", 1:5, collapse = " + "))
fit_js <- cfa(model_js_only, data = df_sem, ordered = paste0("JS", 1:5))
summary(fit_js, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-21 ended normally after 16 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 17
##
## Number of observations 552
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 4.622 6.580
## Degrees of freedom 5 5
## P-value (Unknown) NA 0.254
## Scaling correction factor 0.709
## Shift parameter 0.058
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 2147.066 1722.463
## Degrees of freedom 10 10
## P-value NA 0.000
## Scaling correction factor 1.248
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000 0.999
## Tucker-Lewis Index (TLI) 1.000 0.998
##
## Robust Comparative Fit Index (CFI) 0.995
## Robust Tucker-Lewis Index (TLI) 0.991
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.024
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.057 0.067
## P-value H_0: RMSEA <= 0.050 0.910 0.802
## P-value H_0: RMSEA >= 0.080 0.003 0.012
##
## Robust RMSEA 0.051
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.155
## P-value H_0: Robust RMSEA <= 0.050 0.415
## P-value H_0: Robust RMSEA >= 0.080 0.405
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.025 0.025
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## JS =~
## JS1 1.000 0.818 0.818
## JS2 0.819 0.051 15.915 0.000 0.670 0.670
## JS3 0.969 0.049 19.857 0.000 0.793 0.793
## JS4 1.058 0.048 21.952 0.000 0.866 0.866
## JS5 0.962 0.049 19.645 0.000 0.787 0.787
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## JS1|t1 -0.224 0.054 -4.164 0.000 -0.224 -0.224
## JS1|t2 2.363 0.165 14.312 0.000 2.363 2.363
## JS2|t1 -0.401 0.055 -7.294 0.000 -0.401 -0.401
## JS2|t2 1.655 0.091 18.259 0.000 1.655 1.655
## JS3|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## JS3|t2 -0.333 0.054 -6.112 0.000 -0.333 -0.333
## JS3|t3 1.458 0.080 18.198 0.000 1.458 1.458
## JS4|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## JS4|t2 -0.333 0.054 -6.112 0.000 -0.333 -0.333
## JS4|t3 1.419 0.078 18.112 0.000 1.419 1.419
## JS5|t1 -0.367 0.055 -6.703 0.000 -0.367 -0.367
## JS5|t2 1.360 0.076 17.929 0.000 1.360 1.360
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .JS1 0.331 0.331 0.331
## .JS2 0.551 0.551 0.551
## .JS3 0.371 0.371 0.371
## .JS4 0.251 0.251 0.251
## .JS5 0.380 0.380 0.380
## JS 0.669 0.051 13.002 0.000 1.000 1.000
## JS
## 0.623
## For constructs with categorical indicators, Zumbo et al.`s (2007) "ordinal alpha" is calculated in addition to the standard alpha, which treats ordinal variables as numeric. See Chalmers (2018) for a critique of "alpha.ord" and the response by Zumbo & Kroc (2019). Likewise, average variance extracted is calculated from polychoric (polyserial) not Pearson correlations.
## JS
## alpha 0.8044338
## alpha.ord 0.8887029
## omega 0.8151876
## omega2 0.8151876
## omega3 0.8182810
## avevar 0.6233549
A single-construct CFA was estimated for the JS construct. The scaled CFI was 0.999 and TLI was 0.998, both above 0.90. The scaled RMSEA was 0.024 and the SRMR was 0.025, both meeting their respective thresholds. The model fit is excellent.
Standardized factor loadings ranged from 0.670 (JS2) to 0.866 (JS4). JS2 sits in the acceptable range of 0.50 to 0.70, while all other indicators meet the ideal threshold of >= 0.70. The AVE was 0.623, above the 0.50 threshold, confirming convergent validity. The ordinal alpha was 0.889 and omega was 0.815, both exceeding 0.70. The JS construct meets all reliability and validity criteria.
model_wlb_only <- paste("WLB =~", paste0("WLB", 1:3, collapse = " + "))
fit_wlb <- cfa(model_wlb_only, data = df_sem, ordered = paste0("WLB", 1:3))
summary(fit_wlb, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-21 ended normally after 17 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 12
##
## Number of observations 552
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 0.000 0.000
## Degrees of freedom 0 0
##
## Model Test Baseline Model:
##
## Test statistic 535.427 518.127
## Degrees of freedom 3 3
## P-value NA 0.000
## Scaling correction factor 1.034
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000 1.000
## Tucker-Lewis Index (TLI) 1.000 1.000
##
## Robust Comparative Fit Index (CFI) 1.000
## Robust Tucker-Lewis Index (TLI) 1.000
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.000
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.000 0.000
## P-value H_0: RMSEA <= 0.050 NA NA
## P-value H_0: RMSEA >= 0.080 NA NA
##
## Robust RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000 0.000
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## WLB =~
## WLB1 1.000 0.782 0.782
## WLB2 0.940 0.085 11.004 0.000 0.735 0.735
## WLB3 1.054 0.105 10.048 0.000 0.824 0.824
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## WLB1|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## WLB1|t2 -0.586 0.057 -10.303 0.000 -0.586 -0.586
## WLB1|t3 1.954 0.113 17.255 0.000 1.954 1.954
## WLB2|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## WLB2|t2 -0.276 0.054 -5.097 0.000 -0.276 -0.276
## WLB2|t3 1.233 0.071 17.343 0.000 1.233 1.233
## WLB3|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## WLB3|t2 -0.271 0.054 -5.012 0.000 -0.271 -0.271
## WLB3|t3 1.527 0.083 18.291 0.000 1.527 1.527
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .WLB1 0.388 0.388 0.388
## .WLB2 0.459 0.459 0.459
## .WLB3 0.321 0.321 0.321
## WLB 0.612 0.073 8.398 0.000 1.000 1.000
## WLB
## 0.61
## For constructs with categorical indicators, Zumbo et al.`s (2007) "ordinal alpha" is calculated in addition to the standard alpha, which treats ordinal variables as numeric. See Chalmers (2018) for a critique of "alpha.ord" and the response by Zumbo & Kroc (2019). Likewise, average variance extracted is calculated from polychoric (polyserial) not Pearson correlations.
## WLB
## alpha 0.6947285
## alpha.ord 0.8233283
## omega 0.7132284
## omega2 0.7132284
## omega3 0.7132284
## avevar 0.6103370
The WLB construct was estimated with three indicators. Because this is a just-identified model (degrees of freedom = 0), the model fit indices produce perfect values by definition, and no goodness-of-fit evaluation can be conducted on the model itself. This is a known limitation of three-indicator single-construct models.
Standardized factor loadings ranged from 0.735 (WLB2) to 0.824 (WLB3) in the first-order model output, with all loadings meeting the ideal threshold of >= 0.70. The AVE was 0.610, above 0.50. For reliability, the ordinal alpha was 0.823 and omega was 0.713, both above the 0.70 threshold. The WLB construct satisfies convergent validity and reliability requirements despite the model fit being non-evaluable at this stage.
model_its_only <- paste("ITS =~", paste0("ITS", 1:5, collapse = " + "))
fit_its <- cfa(model_its_only, data = df_sem, ordered = paste0("ITS", 1:5))## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= -3.558432e-17)
## is smaller than zero. This may be a symptom that the model is not
## identified.
## lavaan 0.6-21 ended normally after 16 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 20
##
## Number of observations 552
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 3.635 5.647
## Degrees of freedom 5 5
## P-value (Unknown) NA 0.342
## Scaling correction factor 0.654
## Shift parameter 0.090
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 4487.621 3398.378
## Degrees of freedom 10 10
## P-value NA 0.000
## Scaling correction factor 1.321
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000 1.000
## Tucker-Lewis Index (TLI) 1.001 1.000
##
## Robust Comparative Fit Index (CFI) 1.000
## Robust Tucker-Lewis Index (TLI) 1.004
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.015
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.050 0.063
## P-value H_0: RMSEA <= 0.050 0.949 0.857
## P-value H_0: RMSEA >= 0.080 0.001 0.007
##
## Robust RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.112
## P-value H_0: Robust RMSEA <= 0.050 0.691
## P-value H_0: Robust RMSEA >= 0.080 0.158
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.016 0.016
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ITS =~
## ITS1 1.000 0.900 0.900
## ITS2 0.847 0.033 25.549 0.000 0.762 0.762
## ITS3 0.875 0.032 27.197 0.000 0.787 0.787
## ITS4 0.911 0.033 27.492 0.000 0.820 0.820
## ITS5 0.970 0.032 30.276 0.000 0.872 0.872
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ITS1|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## ITS1|t2 -0.819 0.060 -13.545 0.000 -0.819 -0.819
## ITS1|t3 1.954 0.113 17.255 0.000 1.954 1.954
## ITS2|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## ITS2|t2 -0.635 0.057 -11.046 0.000 -0.635 -0.635
## ITS2|t3 0.763 0.059 12.831 0.000 0.763 0.763
## ITS3|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS3|t2 -0.733 0.059 -12.430 0.000 -0.733 -0.733
## ITS3|t3 0.997 0.064 15.516 0.000 0.997 0.997
## ITS4|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS4|t2 -0.703 0.058 -12.026 0.000 -0.703 -0.703
## ITS4|t3 1.150 0.068 16.808 0.000 1.150 1.150
## ITS5|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS5|t2 -0.775 0.060 -12.991 0.000 -0.775 -0.775
## ITS5|t3 1.305 0.074 17.707 0.000 1.305 1.305
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .ITS1 0.191 0.191 0.191
## .ITS2 0.419 0.419 0.419
## .ITS3 0.381 0.381 0.381
## .ITS4 0.328 0.328 0.328
## .ITS5 0.239 0.239 0.239
## ITS 0.809 0.044 18.484 0.000 1.000 1.000
## ITS
## 0.688
## For constructs with categorical indicators, Zumbo et al.`s (2007) "ordinal alpha" is calculated in addition to the standard alpha, which treats ordinal variables as numeric. See Chalmers (2018) for a critique of "alpha.ord" and the response by Zumbo & Kroc (2019). Likewise, average variance extracted is calculated from polychoric (polyserial) not Pearson correlations.
## ITS
## alpha 0.8359739
## alpha.ord 0.9146666
## omega 0.8406177
## omega2 0.8406177
## omega3 0.8432373
## avevar 0.6883316
The ITS construct CFA yielded excellent fit. The scaled CFI was 1.000 and TLI was 1.000, both exceeding 0.90. The scaled RMSEA was 0.015 and the SRMR was 0.016, both well within acceptable bounds. All five indicators had standardized loadings above 0.70, ranging from 0.762 (ITS2) to 0.900 (ITS1), meeting the ideal threshold. The AVE was 0.688, confirming convergent validity. Both ordinal alpha and omega exceeded 0.70, confirming adequate reliability. The ITS construct performs well across all measurement criteria.
outer_model <- paste(model_pos_only, model_js_only, model_wlb_only, model_its_only, sep = "\n")
fit_outer <- cfa(
outer_model,
data = df_sem,
ordered = ordinal_items
)## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= -2.782422e-16)
## is smaller than zero. This may be a symptom that the model is not
## identified.
## lavaan 0.6-21 ended normally after 44 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 85
##
## Number of observations 552
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 298.448 434.674
## Degrees of freedom 183 183
## P-value (Unknown) NA 0.000
## Scaling correction factor 0.768
## Shift parameter 46.011
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 37915.539 12360.651
## Degrees of freedom 210 210
## P-value NA 0.000
## Scaling correction factor 3.103
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.997 0.979
## Tucker-Lewis Index (TLI) 0.996 0.976
##
## Robust Comparative Fit Index (CFI) 0.892
## Robust Tucker-Lewis Index (TLI) 0.876
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.034 0.050
## 90 Percent confidence interval - lower 0.027 0.044
## 90 Percent confidence interval - upper 0.041 0.056
## P-value H_0: RMSEA <= 0.050 1.000 0.495
## P-value H_0: RMSEA >= 0.080 0.000 0.000
##
## Robust RMSEA 0.104
## 90 Percent confidence interval - lower 0.092
## 90 Percent confidence interval - upper 0.115
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.048 0.048
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS =~
## POS1 1.000 0.803 0.803
## POS2 0.992 0.036 27.182 0.000 0.797 0.797
## POS3 0.958 0.035 27.662 0.000 0.769 0.769
## POS4 1.038 0.039 26.918 0.000 0.834 0.834
## POS5 0.930 0.036 26.085 0.000 0.747 0.747
## POS6 1.061 0.040 26.609 0.000 0.853 0.853
## POS7 1.061 0.035 30.214 0.000 0.852 0.852
## POS8 1.080 0.039 28.035 0.000 0.868 0.868
## JS =~
## JS1 1.000 0.788 0.788
## JS2 0.930 0.046 20.263 0.000 0.733 0.733
## JS3 1.069 0.046 23.347 0.000 0.843 0.843
## JS4 1.060 0.046 22.914 0.000 0.835 0.835
## JS5 0.946 0.041 22.885 0.000 0.745 0.745
## WLB =~
## WLB1 1.000 0.848 0.848
## WLB2 0.903 0.042 21.412 0.000 0.766 0.766
## WLB3 0.867 0.041 21.262 0.000 0.735 0.735
## ITS =~
## ITS1 1.000 0.917 0.917
## ITS2 0.873 0.028 31.515 0.000 0.801 0.801
## ITS3 0.868 0.025 34.587 0.000 0.796 0.796
## ITS4 0.887 0.027 32.924 0.000 0.814 0.814
## ITS5 0.899 0.026 35.238 0.000 0.825 0.825
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS ~~
## JS 0.403 0.029 13.729 0.000 0.637 0.637
## WLB 0.383 0.032 11.893 0.000 0.562 0.562
## ITS 0.597 0.033 17.911 0.000 0.810 0.810
## JS ~~
## WLB 0.548 0.035 15.582 0.000 0.820 0.820
## ITS 0.513 0.033 15.431 0.000 0.709 0.709
## WLB ~~
## ITS 0.614 0.033 18.412 0.000 0.790 0.790
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS1|t1 -0.939 0.063 -14.928 0.000 -0.939 -0.939
## POS1|t2 1.752 0.097 18.067 0.000 1.752 1.752
## POS2|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS2|t2 -0.831 0.061 -13.702 0.000 -0.831 -0.831
## POS2|t3 0.538 0.056 9.556 0.000 0.538 0.538
## POS3|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS3|t2 -0.953 0.063 -15.077 0.000 -0.953 -0.953
## POS3|t3 1.107 0.067 16.486 0.000 1.107 1.107
## POS4|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## POS4|t2 -0.932 0.063 -14.853 0.000 -0.932 -0.932
## POS4|t3 1.066 0.066 16.149 0.000 1.066 1.066
## POS5|t1 -0.967 0.064 -15.224 0.000 -0.967 -0.967
## POS5|t2 1.035 0.065 15.872 0.000 1.035 1.035
## POS6|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## POS6|t2 -0.864 0.061 -14.091 0.000 -0.864 -0.864
## POS6|t3 0.997 0.064 15.516 0.000 0.997 0.997
## POS7|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## POS7|t2 -0.831 0.061 -13.702 0.000 -0.831 -0.831
## POS7|t3 1.035 0.065 15.872 0.000 1.035 1.035
## POS8|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS8|t2 -0.884 0.062 -14.322 0.000 -0.884 -0.884
## POS8|t3 1.186 0.070 17.053 0.000 1.186 1.186
## JS1|t1 -0.224 0.054 -4.164 0.000 -0.224 -0.224
## JS1|t2 2.363 0.165 14.312 0.000 2.363 2.363
## JS2|t1 -0.401 0.055 -7.294 0.000 -0.401 -0.401
## JS2|t2 1.655 0.091 18.259 0.000 1.655 1.655
## JS3|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## JS3|t2 -0.333 0.054 -6.112 0.000 -0.333 -0.333
## JS3|t3 1.458 0.080 18.198 0.000 1.458 1.458
## JS4|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## JS4|t2 -0.333 0.054 -6.112 0.000 -0.333 -0.333
## JS4|t3 1.419 0.078 18.112 0.000 1.419 1.419
## JS5|t1 -0.367 0.055 -6.703 0.000 -0.367 -0.367
## JS5|t2 1.360 0.076 17.929 0.000 1.360 1.360
## WLB1|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## WLB1|t2 -0.586 0.057 -10.303 0.000 -0.586 -0.586
## WLB1|t3 1.954 0.113 17.255 0.000 1.954 1.954
## WLB2|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## WLB2|t2 -0.276 0.054 -5.097 0.000 -0.276 -0.276
## WLB2|t3 1.233 0.071 17.343 0.000 1.233 1.233
## WLB3|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## WLB3|t2 -0.271 0.054 -5.012 0.000 -0.271 -0.271
## WLB3|t3 1.527 0.083 18.291 0.000 1.527 1.527
## ITS1|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## ITS1|t2 -0.819 0.060 -13.545 0.000 -0.819 -0.819
## ITS1|t3 1.954 0.113 17.255 0.000 1.954 1.954
## ITS2|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## ITS2|t2 -0.635 0.057 -11.046 0.000 -0.635 -0.635
## ITS2|t3 0.763 0.059 12.831 0.000 0.763 0.763
## ITS3|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS3|t2 -0.733 0.059 -12.430 0.000 -0.733 -0.733
## ITS3|t3 0.997 0.064 15.516 0.000 0.997 0.997
## ITS4|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS4|t2 -0.703 0.058 -12.026 0.000 -0.703 -0.703
## ITS4|t3 1.150 0.068 16.808 0.000 1.150 1.150
## ITS5|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS5|t2 -0.775 0.060 -12.991 0.000 -0.775 -0.775
## ITS5|t3 1.305 0.074 17.707 0.000 1.305 1.305
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .POS1 0.354 0.354 0.354
## .POS2 0.365 0.365 0.365
## .POS3 0.408 0.408 0.408
## .POS4 0.304 0.304 0.304
## .POS5 0.442 0.442 0.442
## .POS6 0.273 0.273 0.273
## .POS7 0.273 0.273 0.273
## .POS8 0.247 0.247 0.247
## .JS1 0.379 0.379 0.379
## .JS2 0.463 0.463 0.463
## .JS3 0.290 0.290 0.290
## .JS4 0.302 0.302 0.302
## .JS5 0.444 0.444 0.444
## .WLB1 0.280 0.280 0.280
## .WLB2 0.413 0.413 0.413
## .WLB3 0.459 0.459 0.459
## .ITS1 0.159 0.159 0.159
## .ITS2 0.359 0.359 0.359
## .ITS3 0.367 0.367 0.367
## .ITS4 0.338 0.338 0.338
## .ITS5 0.320 0.320 0.320
## POS 0.646 0.042 15.539 0.000 1.000 1.000
## JS 0.621 0.045 13.948 0.000 1.000 1.000
## WLB 0.720 0.052 13.858 0.000 1.000 1.000
## ITS 0.841 0.038 21.853 0.000 1.000 1.000
## POS JS WLB ITS
## 0.667 0.624 0.616 0.692
## For constructs with categorical indicators, Zumbo et al.`s (2007) "ordinal alpha" is calculated in addition to the standard alpha, which treats ordinal variables as numeric. See Chalmers (2018) for a critique of "alpha.ord" and the response by Zumbo & Kroc (2019). Likewise, average variance extracted is calculated from polychoric (polyserial) not Pearson correlations.
## POS JS WLB ITS
## alpha 0.8841865 0.8044338 0.6947285 0.8359739
## alpha.ord 0.9388436 0.8887029 0.8233283 0.9146666
## omega 0.8897832 0.8155385 0.7129708 0.8436781
## omega2 0.8897832 0.8155385 0.7129708 0.8436781
## omega3 0.8993271 0.8214535 0.7127891 0.8503059
## avevar 0.6667180 0.6243000 0.6157428 0.6917354
A combined first-order CFA was estimated with all four latent constructs modeled simultaneously. The scaled CFI was 0.979 and TLI was 0.976, both above the 0.90 threshold. The scaled RMSEA was 0.050, which sits exactly at the boundary of the threshold of < 0.08 and is considered acceptable. The SRMR was 0.048, just below the threshold of 0.05. Overall, the measurement model demonstrates acceptable fit.
All standardized factor loadi0.8890.889ngs met or exceeded 0.70 across the four constructs. For POS, loadings ranged from 0.747 to 0.868. For JS, loadings ranged from 0.733 to 0.843. For WLB, loadings ranged from 0.735 to 0.848. For ITS, loadings ranged from 0.796 to 0.917. All loadings satisfy the ideal threshold, confirming that each indicator adequately reflects its respective construct.
All constructs exceeded the recommended thresholds of 0.50 for AVE and 0.70 for reliability. Specifically, Perceived Organizational Support (POS) demonstrated strong validity and reliability (AVE = 0.667, Ordinal \(\alpha\) = 0.939, \(\omega\) = 0.890). Job Satisfaction (JS) also met all criteria (AVE = 0.624, Ordinal \(\alpha\) = 0.889, \(\omega\) = 0.816). Work-Life Balance (WLB) showed acceptable internal consistency and variance extraction (AVE = 0.616, Ordinal \(\alpha\) = 0.823, \(\omega\) = 0.713). Finally, Intention to Stay (ITS) exhibited robust psychometric properties (AVE = 0.692, Ordinal \(\alpha\) = 0.915, \(\omega\) = 0.844).
## POS JS WLB ITS
## POS 1.000
## JS 0.631 1.000
## WLB 0.535 0.824 1.000
## ITS 0.811 0.713 0.774 1.000
For the HTMT ratio, the highest value was observed between JS and WLB at 0.824, which remains below the stricter threshold of 0.85. All other HTMT values were also below 0.85. This confirms that each construct is empirically distinct from the others.
##
## --- Construct Correlations ---
## POS JS WLB ITS
## POS 1.000
## JS 0.637 1.000
## WLB 0.562 0.820 1.000
## ITS 0.810 0.709 0.790 1.000
ave_values <- semTools::AVE(fit_outer)
sqrt_ave <- sqrt(ave_values)
cat("\n--- Square Root of AVE ---\n")##
## --- Square Root of AVE ---
## POS JS WLB ITS
## 0.817 0.790 0.785 0.832
For the Fornell-Larcker criterion, the square roots of the AVE values were POS = 0.817, JS = 0.790, WLB = 0.785, and ITS = 0.832. In each case, the square root of a construct’s AVE exceeds its correlations with every other construct. For example, ITS correlates with POS at 0.810 and with WLB at 0.790, but the square root of the AVE for both ITS (0.832) and POS (0.817) and WLB (0.785) still satisfies the criterion when compared within their respective rows and columns.
Discriminant validity is therefore supported across all constructs, confirming that the four latent variables measure distinct dimensions.
inner_model <- "ITS ~ POS + JS + WLB"
full_model <- paste(outer_model, inner_model, sep = "\n")
fit_full <- sem(model = full_model,
data = df_sem,
ordered = ordinal_items)## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= -2.760201e-16)
## is smaller than zero. This may be a symptom that the model is not
## identified.
semPaths(fit_full,
whatLabels = "std",
layout = "tree",
rotation = 2,
intercepts = FALSE,
thresholds = FALSE,
sizeMan = 5,
sizeLat = 8,
edge.label.cex = 0.7,
color = list(lat = "lightblue", man = "lightgreen"),
edge.color = "black",
mar = c(2, 5, 2, 5))The full structural model was estimated using the DWLS estimator with all 21 ordinal indicators and the four latent constructs specified in the outer model. The inner model regresses Employee Intention to Stay (ITS) on Perceived Organizational Support (POS), Job Satisfaction (JS), and Work-Life Balance (WLB).
Acceptable thresholds are CFI and TLI >= 0.90, RMSEA < 0.08, and SRMR < 0.05.
## lavaan 0.6-21 ended normally after 50 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 85
##
## Number of observations 552
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 298.448 434.674
## Degrees of freedom 183 183
## P-value (Unknown) NA 0.000
## Scaling correction factor 0.768
## Shift parameter 46.011
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 37915.539 12360.651
## Degrees of freedom 210 210
## P-value NA 0.000
## Scaling correction factor 3.103
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.997 0.979
## Tucker-Lewis Index (TLI) 0.996 0.976
##
## Robust Comparative Fit Index (CFI) 0.892
## Robust Tucker-Lewis Index (TLI) 0.876
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.034 0.050
## 90 Percent confidence interval - lower 0.027 0.044
## 90 Percent confidence interval - upper 0.041 0.056
## P-value H_0: RMSEA <= 0.050 1.000 0.495
## P-value H_0: RMSEA >= 0.080 0.000 0.000
##
## Robust RMSEA 0.104
## 90 Percent confidence interval - lower 0.092
## 90 Percent confidence interval - upper 0.115
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.048 0.048
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS =~
## POS1 1.000 0.803 0.803
## POS2 0.992 0.036 27.182 0.000 0.797 0.797
## POS3 0.958 0.035 27.662 0.000 0.769 0.769
## POS4 1.038 0.039 26.918 0.000 0.834 0.834
## POS5 0.930 0.036 26.085 0.000 0.747 0.747
## POS6 1.061 0.040 26.609 0.000 0.853 0.853
## POS7 1.061 0.035 30.214 0.000 0.852 0.852
## POS8 1.080 0.039 28.035 0.000 0.868 0.868
## JS =~
## JS1 1.000 0.788 0.788
## JS2 0.930 0.046 20.263 0.000 0.733 0.733
## JS3 1.069 0.046 23.347 0.000 0.843 0.843
## JS4 1.060 0.046 22.914 0.000 0.835 0.835
## JS5 0.946 0.041 22.885 0.000 0.745 0.745
## WLB =~
## WLB1 1.000 0.848 0.848
## WLB2 0.903 0.042 21.412 0.000 0.766 0.766
## WLB3 0.867 0.041 21.262 0.000 0.735 0.735
## ITS =~
## ITS1 1.000 0.917 0.917
## ITS2 0.873 0.028 31.515 0.000 0.801 0.801
## ITS3 0.868 0.025 34.587 0.000 0.796 0.796
## ITS4 0.887 0.027 32.924 0.000 0.814 0.814
## ITS5 0.899 0.026 35.238 0.000 0.825 0.825
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ITS ~
## POS 0.644 0.066 9.719 0.000 0.564 0.564
## JS -0.133 0.123 -1.083 0.279 -0.114 -0.114
## WLB 0.612 0.111 5.500 0.000 0.566 0.566
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS ~~
## JS 0.403 0.029 13.729 0.000 0.637 0.637
## WLB 0.383 0.032 11.893 0.000 0.562 0.562
## JS ~~
## WLB 0.548 0.035 15.582 0.000 0.820 0.820
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## POS1|t1 -0.939 0.063 -14.928 0.000 -0.939 -0.939
## POS1|t2 1.752 0.097 18.067 0.000 1.752 1.752
## POS2|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS2|t2 -0.831 0.061 -13.702 0.000 -0.831 -0.831
## POS2|t3 0.538 0.056 9.556 0.000 0.538 0.538
## POS3|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS3|t2 -0.953 0.063 -15.077 0.000 -0.953 -0.953
## POS3|t3 1.107 0.067 16.486 0.000 1.107 1.107
## POS4|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## POS4|t2 -0.932 0.063 -14.853 0.000 -0.932 -0.932
## POS4|t3 1.066 0.066 16.149 0.000 1.066 1.066
## POS5|t1 -0.967 0.064 -15.224 0.000 -0.967 -0.967
## POS5|t2 1.035 0.065 15.872 0.000 1.035 1.035
## POS6|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## POS6|t2 -0.864 0.061 -14.091 0.000 -0.864 -0.864
## POS6|t3 0.997 0.064 15.516 0.000 0.997 0.997
## POS7|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## POS7|t2 -0.831 0.061 -13.702 0.000 -0.831 -0.831
## POS7|t3 1.035 0.065 15.872 0.000 1.035 1.035
## POS8|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## POS8|t2 -0.884 0.062 -14.322 0.000 -0.884 -0.884
## POS8|t3 1.186 0.070 17.053 0.000 1.186 1.186
## JS1|t1 -0.224 0.054 -4.164 0.000 -0.224 -0.224
## JS1|t2 2.363 0.165 14.312 0.000 2.363 2.363
## JS2|t1 -0.401 0.055 -7.294 0.000 -0.401 -0.401
## JS2|t2 1.655 0.091 18.259 0.000 1.655 1.655
## JS3|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## JS3|t2 -0.333 0.054 -6.112 0.000 -0.333 -0.333
## JS3|t3 1.458 0.080 18.198 0.000 1.458 1.458
## JS4|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## JS4|t2 -0.333 0.054 -6.112 0.000 -0.333 -0.333
## JS4|t3 1.419 0.078 18.112 0.000 1.419 1.419
## JS5|t1 -0.367 0.055 -6.703 0.000 -0.367 -0.367
## JS5|t2 1.360 0.076 17.929 0.000 1.360 1.360
## WLB1|t1 -2.685 0.236 -11.374 0.000 -2.685 -2.685
## WLB1|t2 -0.586 0.057 -10.303 0.000 -0.586 -0.586
## WLB1|t3 1.954 0.113 17.255 0.000 1.954 1.954
## WLB2|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## WLB2|t2 -0.276 0.054 -5.097 0.000 -0.276 -0.276
## WLB2|t3 1.233 0.071 17.343 0.000 1.233 1.233
## WLB3|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## WLB3|t2 -0.271 0.054 -5.012 0.000 -0.271 -0.271
## WLB3|t3 1.527 0.083 18.291 0.000 1.527 1.527
## ITS1|t1 -2.909 0.313 -9.306 0.000 -2.909 -2.909
## ITS1|t2 -0.819 0.060 -13.545 0.000 -0.819 -0.819
## ITS1|t3 1.954 0.113 17.255 0.000 1.954 1.954
## ITS2|t1 -2.547 0.201 -12.664 0.000 -2.547 -2.547
## ITS2|t2 -0.635 0.057 -11.046 0.000 -0.635 -0.635
## ITS2|t3 0.763 0.059 12.831 0.000 0.763 0.763
## ITS3|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS3|t2 -0.733 0.059 -12.430 0.000 -0.733 -0.733
## ITS3|t3 0.997 0.064 15.516 0.000 0.997 0.997
## ITS4|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS4|t2 -0.703 0.058 -12.026 0.000 -0.703 -0.703
## ITS4|t3 1.150 0.068 16.808 0.000 1.150 1.150
## ITS5|t1 -2.363 0.165 -14.312 0.000 -2.363 -2.363
## ITS5|t2 -0.775 0.060 -12.991 0.000 -0.775 -0.775
## ITS5|t3 1.305 0.074 17.707 0.000 1.305 1.305
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .POS1 0.354 0.354 0.354
## .POS2 0.365 0.365 0.365
## .POS3 0.408 0.408 0.408
## .POS4 0.304 0.304 0.304
## .POS5 0.442 0.442 0.442
## .POS6 0.273 0.273 0.273
## .POS7 0.273 0.273 0.273
## .POS8 0.247 0.247 0.247
## .JS1 0.379 0.379 0.379
## .JS2 0.463 0.463 0.463
## .JS3 0.290 0.290 0.290
## .JS4 0.302 0.302 0.302
## .JS5 0.444 0.444 0.444
## .WLB1 0.280 0.280 0.280
## .WLB2 0.413 0.413 0.413
## .WLB3 0.459 0.459 0.459
## .ITS1 0.159 0.159 0.159
## .ITS2 0.359 0.359 0.359
## .ITS3 0.367 0.367 0.367
## .ITS4 0.338 0.338 0.338
## .ITS5 0.320 0.320 0.320
## POS 0.646 0.042 15.539 0.000 1.000 1.000
## JS 0.621 0.045 13.948 0.000 1.000 1.000
## WLB 0.720 0.052 13.858 0.000 1.000 1.000
## .ITS 0.149 0.034 4.441 0.000 0.177 0.177
The scaled CFI was 0.979 and TLI was 0.976, both exceeding the 0.90 threshold. The scaled RMSEA was 0.050, sitting exactly at the boundary of the acceptable range, with a 90% confidence interval of [0.044, 0.056]. The SRMR was 0.048, falling just below the 0.05 threshold. Overall, the scaled fit indicators suggest that the full structural model demonstrates acceptable fit to the data.
A path is considered statistically significant at p < 0.05. Standardized coefficients (Std.all) are interpreted as the expected standard deviation change in ITS for a one standard deviation increase in the predictor, holding all other predictors constant.
The table below presents the standardized structural path coefficients estimating the direct effect of each predictor on ITS.
parameterEstimates(fit_full, standardized = TRUE) |>
filter(op == "~") |>
select(lhs, op, rhs, std.all, se, z, pvalue) |>
rename(
Outcome = lhs,
` ` = op,
Predictor = rhs,
`Std. Beta` = std.all,
SE = se,
`z-value` = z,
`p-value` = pvalue
)POS showed a positive and statistically significant effect on ITS (β = 0.564, SE = 0.066, z = 9.719, p < 0.001), indicating that higher perceived organizational support is associated with stronger employee intention to stay. WLB similarly showed a positive and significant effect (β = 0.566, SE = 0.111, z = 5.500, p < 0.001), suggesting that better work-life balance significantly strengthens retention intention. JS, however, showed a negative and non-significant effect (β = -0.114, SE = 0.123, z = -1.083, p = 0.279), indicating that job satisfaction does not single-handedly predict intention to stay when POS and WLB are controlled for.
The R² value indicates the proportion of variance in Employee Intention to Stay (ITS) that is jointly explained by POS, JS, and WLB in the structural model. A higher R² reflects greater collective explanatory power of the predictors.
## [1] 0.8225456
POS, JS, and WLB combined explain 82.3% of the variance in Employee Intention to Stay, indicating that the structural model has very strong explanatory power. This suggests that perceived organizational support, job satisfaction, and work-life balance collectively responsible for the large majority of what drives employees’ intention to remain with the organization.
Cohen’s f² measures the practical significance of each predictor’s unique contribution to explaining ITS variance. Thresholds follow Cohen (1988): f² >= 0.02 small, >= 0.15 medium, >= 0.35 large.
fit_no_pos <- sem(
model = paste(outer_model, "ITS ~ JS + WLB", sep = "\n"),
data = df_sem,
ordered = ordinal_items
)## Warning: lavaan->lav_lavaan_step11_estoptim():
## the optimizer (NLMINB) claimed the model converged, but not all elements
## of the gradient are (near) zero; the optimizer may not have found a local
## solution use check.gradient = FALSE to skip this check.
fit_no_js <- sem(
model = paste(outer_model, "ITS ~ POS + WLB", sep = "\n"),
data = df_sem,
ordered = ordinal_items
)## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= -2.682349e-16)
## is smaller than zero. This may be a symptom that the model is not
## identified.
fit_no_wlb <- sem(
model = paste(outer_model, "ITS ~ POS + JS", sep = "\n"),
data = df_sem,
ordered = ordinal_items
)## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= -2.837282e-16)
## is smaller than zero. This may be a symptom that the model is not
## identified.
r2_full <- lavInspect(fit_full, what = "r2")[["ITS"]]
r2_no_pos <- lavInspect(fit_no_pos, what = "r2")[["ITS"]]
r2_no_js <- lavInspect(fit_no_js, what = "r2")[["ITS"]]
r2_no_wlb <- lavInspect(fit_no_wlb, what = "r2")[["ITS"]]
f2 <- function(r2_full, r2_reduced) (r2_full - r2_reduced) / (1 - r2_full)
data.frame(
Predictor = c("POS", "JS", "WLB"),
R2_full = round(r2_full, 3),
R2_reduced = round(c(r2_no_pos, r2_no_js, r2_no_wlb), 3),
f2 = round(c(
f2(r2_full, r2_no_pos),
f2(r2_full, r2_no_js),
f2(r2_full, r2_no_wlb)
), 3),
Interpretation = c(
ifelse(f2(r2_full, r2_no_pos) >= 0.35, "Large",
ifelse(f2(r2_full, r2_no_pos) >= 0.15, "Medium", "Small")),
ifelse(f2(r2_full, r2_no_js) >= 0.35, "Large",
ifelse(f2(r2_full, r2_no_js) >= 0.15, "Medium", "Small")),
ifelse(f2(r2_full, r2_no_wlb) >= 0.35, "Large",
ifelse(f2(r2_full, r2_no_wlb) >= 0.15, "Medium", "Small"))
)
)The f² for POS demonstrate a rather impossible result (f² = 111.856) because the reduced model excluding POS failed to converge. Without POS, the structural paths for JS and WLB inflated to extreme values (141.356 and -144.078 respectively), and the residual variance of ITS blew up to 16.846, producing a negative R² of -19.027. This instability is possibly due to the high intercorrelation between JS and WLB (r = 0.820), which makes the model unidentified when POS is removed. Consequently, the f² value for POS cannot be meaningfully interpreted and is excluded from the discussion. For JS, the f² was 0.124, indicating a small effect. For WLB, the f² was 0.267, indicating a medium effect. These results suggest that while both JS and WLB contribute uniquely to explaining ITS variance beyond the other predictors, WLB carries greater incremental explanatory power.
The mediation analysis tests whether JS and WLB act as mediators in the relationship between POS and ITS. POS is treated as the exogenous predictor, JS and WLB as potential mediators, and ITS as the outcome. Indirect effects are estimated using bootstrap resampling (1000 iterations) with bias-corrected accelerated (BCa) 95% confidence intervals. A mediation effect is considered significant if the confidence interval excludes zero.
mediation_model <- "
JS ~ a1 * POS
WLB ~ a2 * POS
ITS ~ b1 * JS + b2 * WLB + c * POS
indirect_js := a1 * b1
indirect_wlb := a2 * b2
total_indirect := a1 * b1 + a2 * b2
total := c + a1 * b1 + a2 * b2
"
full_mediation_model <- paste(outer_model, mediation_model, sep = "\n")
fit_mediation <- sem(
model = full_mediation_model,
data = df_sem,
se = "bootstrap",
bootstrap = 1000
)
parameterEstimates(fit_mediation, boot.ci.type = "bca.simple", level = 0.95) |>
filter(op == ":=") |>
select(label, est, se, z, pvalue, ci.lower, ci.upper) |>
rename(
Effect = label,
Estimate = est,
SE = se,
`z-value` = z,
`p-value` = pvalue,
`95% CI Lower` = ci.lower,
`95% CI Upper` = ci.upper
)The mediation analysis tested whether JS and WLB mediate the relationship between POS and ITS using bootstrap resampling with BCa 95% confidence intervals. A mediation effect is considered significant when the confidence interval excludes zero.
The indirect effect of POS on ITS through JS was non-significant (β = 0.048, SE = 0.052, z = 0.919, p = 0.358, 95% CI [-0.066, 0.137]), as the confidence interval includes zero. This indicates that JS does not function as a mediator in the POS-ITS relationship. The indirect effect through WLB, however, was positive and significant (β = 0.307, SE = 0.060, z = 5.102, p < 0.001, 95% CI [0.203, 0.441]), confirming that WLB partially mediates the effect of POS on ITS. The total indirect effect was also significant (β = 0.355, SE = 0.064, z = 5.546, p < 0.001, 95% CI [0.231, 0.487]), driven primarily by the WLB pathway. The total effect of POS on ITS remained significant (β = 0.976, SE = 0.095, z = 10.244, p < 0.001, 95% CI [0.809, 1.197]), confirming that POS retains a strong direct effect on ITS independent of the mediation pathways.
This study used CB-SEM to examine how Perceived Organizational Support (POS), Job Satisfaction (JS), and Work-Life Balance (WLB) affect employees’ Intention to Stay (ITS). The analysis covered three stages: checking assumptions, validating the measurement model, and testing structural paths with mediation.
Before modeling, multivariate outliers were removed using Mahalanobis distance. The normality assumption was violated across all 21 indicators, so the DWLS estimator with robust standard errors was used throughout. The KMO value of 0.95 confirmed the data were well-suited for factor analysis, and all VIF values stayed below 3.3, ruling out multicollinearity.
The measurement model performed well. All four constructs met the AVE threshold of 0.50, reliability thresholds of 0.70, and factor loadings at or above 0.70. Discriminant validity held under both the HTMT criterion (highest ratio = 0.824, below 0.85) and the Fornell-Larcker criterion. The combined CFA fit was acceptable (CFI = 0.979, TLI = 0.976, RMSEA = 0.050, SRMR = 0.048).
In the structural model, POS and WLB each had a significant positive effect on ITS (POS: β = 0.564, p < 0.001; WLB: β = 0.566, p < 0.001). JS did not show a significant direct effect once POS and WLB were controlled for (β = -0.114, p = 0.279). The three predictors together explained 82.3% of the variance in ITS. Effect sizes showed WLB had a medium incremental contribution (f² = 0.267) and JS had a small one (f² = 0.124). POS could not be evaluated for effect size due to model instability upon its removal, which itself reflects how central POS is to the model structure.
The mediation analysis found that WLB significantly mediates the effect of POS on ITS (indirect effect β = 0.307, 95% CI [0.203, 0.441]). This means part of how organizational support increases retention is by improving employees’ work-life balance. JS did not mediate this relationship (indirect effect β = 0.048, 95% CI [-0.066, 0.137]). POS also retained a strong direct effect on ITS (total effect β = 0.976, p < 0.001), independent of the mediation paths.
Overall, the findings suggest that POS and WLB are the most important drivers of employee retention. Organizations looking to retain staff should focus on making employees feel supported and on creating conditions that allow for sustainable work-life balance. JS, while valuable, does not independently predict retention once support and balance are in place. This study shows that CB-SEM, when built on solid measurement validation, can produce both statistically sound and practically useful insights into what keeps employees in an organization.