library(readxl)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.5.3
library(mice)
## Warning: package 'mice' was built under R version 4.5.3
##
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
##
## filter
## The following objects are masked from 'package:base':
##
## cbind, rbind
library(MVN)
## Warning: package 'MVN' was built under R version 4.5.3
library(psych)
##
## Attaching package: 'psych'
## The following object is masked from 'package:MVN':
##
## mardia
library(lavaan)
## 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.
##
## Attaching package: 'lavaan'
## The following object is masked from 'package:psych':
##
## cor2cov
library(semPlot)
## Warning: package 'semPlot' was built under R version 4.5.3
library(semTools)
## 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 objects are masked from 'package:psych':
##
## reliability, skew
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
library(corrplot)
## corrplot 0.95 loaded
library(flextable)
## Warning: package 'flextable' was built under R version 4.5.3
library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.5.3
cat("\n========== TAHAP 1: PERSIAPAN & SELEKSI VARIABEL ==========\n")
##
## ========== TAHAP 1: PERSIAPAN & SELEKSI VARIABEL ==========
# --- 1.1 Load data -----------------------------------------------------------
df_raw <- read_excel("data6.xlsx", sheet = "Data")
cat("Dimensi data raw:", nrow(df_raw), "baris x", ncol(df_raw), "kolom\n")
## Dimensi data raw: 200 baris x 255 kolom
indicators <- list(
JOB_DEMANDS = paste0("JDR", 11:20), # 10 item demands
JOB_RESOURCES = paste0("JDR", 21:30), # 10 item resources
WELLBEING = paste0("FW", 1:10), # 10 item wellbeing
WORK_ENGAGE = paste0("WMI", 1:10), # 10 item engagement
ITL = c("ITL1", "ITL2", "ITL3"), # 3 item intention to leave
JOB_PERF = c("JP1", "JP2", "JP3", "JP4") # 4 item job performance
)
all_items <- unlist(indicators, use.names = FALSE)
cat("Total item dipilih:", length(all_items), "\n")
## Total item dipilih: 47
cat("Konstruk:", paste(names(indicators), collapse = ", "), "\n")
## Konstruk: JOB_DEMANDS, JOB_RESOURCES, WELLBEING, WORK_ENGAGE, ITL, JOB_PERF
# --- 1.3 Subset data ---------------------------------------------------------
df <- df_raw %>%
select(No, Gender, Age, all_of(all_items)) %>%
mutate(across(all_of(all_items), as.numeric))
cat("Dimensi data subset:", nrow(df), "baris x", ncol(df), "kolom\n")
## Dimensi data subset: 200 baris x 50 kolom
# --- 1.4 Statistik deskriptif awal -------------------------------------------
cat("\n--- Statistik Deskriptif Awal ---\n")
##
## --- Statistik Deskriptif Awal ---
desc_stats <- df %>%
select(all_of(all_items)) %>%
describe() %>% # dari package psych
round(3)
print(desc_stats[, c("n", "mean", "sd", "min", "max", "skew", "kurtosis")])
## n mean sd min max skew kurtosis
## JDR11 200 3.86 0.97 1 5 -0.34 -0.57
## JDR12 197 3.46 1.14 1 5 -0.34 -0.61
## JDR13 200 3.69 1.08 1 5 -0.40 -0.50
## JDR14 199 3.72 1.03 1 5 -0.54 -0.15
## JDR15 198 4.31 0.82 2 5 -1.00 0.21
## JDR16 199 4.10 1.00 1 5 -0.89 0.18
## JDR17 199 2.77 1.29 1 5 0.17 -1.00
## JDR18 200 3.07 1.28 1 5 -0.12 -1.01
## JDR19 197 2.81 1.17 1 5 0.08 -0.70
## JDR20 200 3.69 1.10 1 5 -0.60 -0.24
## JDR21 200 4.04 1.01 1 5 -0.77 -0.19
## JDR22 200 4.04 0.97 1 5 -0.91 0.44
## JDR23 199 4.32 0.85 1 5 -1.20 1.23
## JDR24 200 4.07 1.00 1 5 -0.93 0.29
## JDR25 200 3.95 1.01 1 5 -0.54 -0.65
## JDR26 199 4.00 1.18 1 5 -0.97 -0.02
## JDR27 200 4.32 0.93 1 5 -1.23 0.80
## JDR28 199 3.92 1.14 1 5 -0.77 -0.36
## JDR29 199 3.86 1.03 1 5 -0.56 -0.52
## JDR30 199 4.35 0.83 1 5 -1.19 0.98
## FW1 200 4.74 1.02 1 6 -1.47 2.92
## FW2 199 4.54 1.11 1 6 -1.13 1.32
## FW3 199 5.04 1.03 1 6 -1.39 2.48
## FW4 200 2.75 1.31 1 6 0.73 -0.18
## FW5 200 1.99 1.23 1 6 1.46 1.75
## FW6 200 2.02 1.25 1 6 1.25 0.95
## FW7 199 4.54 1.29 1 6 -1.18 0.96
## FW8 200 4.72 1.30 1 6 -1.22 0.98
## FW9 199 4.97 0.97 1 6 -1.05 1.26
## FW10 195 4.58 1.12 1 6 -0.94 0.49
## WMI1 175 4.37 0.80 1 5 -1.49 2.87
## WMI2 198 4.28 0.98 1 6 -1.66 2.80
## WMI3 198 4.11 1.45 1 5 -1.35 0.17
## WMI4 198 4.37 0.91 1 5 -1.93 4.13
## WMI5 198 4.48 0.80 1 6 -2.13 6.00
## WMI6 196 4.59 0.82 1 6 -2.27 6.57
## WMI7 198 4.37 0.86 1 6 -1.79 4.38
## WMI8 197 4.26 0.95 1 6 -1.39 2.14
## WMI9 198 4.35 0.86 1 6 -1.45 2.94
## WMI10 198 4.59 0.79 1 6 -2.18 6.48
## ITL1 199 2.43 1.30 1 6 0.65 -0.44
## ITL2 199 2.14 1.08 1 5 0.81 0.15
## ITL3 199 2.28 1.26 1 5 0.76 -0.44
## JP1 194 7.47 1.66 1 10 -0.89 1.54
## JP2 194 7.65 1.68 1 10 -0.72 0.38
## JP3 195 7.80 1.55 1 10 -0.78 0.98
## JP4 195 7.91 1.65 1 10 -1.17 1.81
# Simpan ke Excel
wb <- createWorkbook()
addWorksheet(wb, "Deskriptif")
writeData(wb, "Deskriptif", desc_stats, rowNames = TRUE)
# =============================================================================
# TAHAP 2: PEMBERSIHAN DATA & MISSING VALUES
# =============================================================================
cat("\n========== TAHAP 2: PEMBERSIHAN DATA & MISSING VALUES ==========\n")
##
## ========== TAHAP 2: PEMBERSIHAN DATA & MISSING VALUES ==========
# --- 2.1 Identifikasi missing values -----------------------------------------
df_items <- df %>% select(all_of(all_items))
missing_summary <- data.frame(
Variabel = names(df_items),
N_Missing = colSums(is.na(df_items)),
Pct_Missing = round(colMeans(is.na(df_items)) * 100, 2)
) %>% filter(N_Missing > 0) %>% arrange(desc(Pct_Missing))
cat("Variabel dengan missing values:\n")
## Variabel dengan missing values:
print(missing_summary)
## Variabel N_Missing Pct_Missing
## WMI1 WMI1 25 12.5
## JP1 JP1 6 3.0
## JP2 JP2 6 3.0
## FW10 FW10 5 2.5
## JP3 JP3 5 2.5
## JP4 JP4 5 2.5
## WMI6 WMI6 4 2.0
## JDR12 JDR12 3 1.5
## JDR19 JDR19 3 1.5
## WMI8 WMI8 3 1.5
## JDR15 JDR15 2 1.0
## WMI2 WMI2 2 1.0
## WMI3 WMI3 2 1.0
## WMI4 WMI4 2 1.0
## WMI5 WMI5 2 1.0
## WMI7 WMI7 2 1.0
## WMI9 WMI9 2 1.0
## WMI10 WMI10 2 1.0
## JDR14 JDR14 1 0.5
## JDR16 JDR16 1 0.5
## JDR17 JDR17 1 0.5
## JDR23 JDR23 1 0.5
## JDR26 JDR26 1 0.5
## JDR28 JDR28 1 0.5
## JDR29 JDR29 1 0.5
## JDR30 JDR30 1 0.5
## FW2 FW2 1 0.5
## FW3 FW3 1 0.5
## FW7 FW7 1 0.5
## FW9 FW9 1 0.5
## ITL1 ITL1 1 0.5
## ITL2 ITL2 1 0.5
## ITL3 ITL3 1 0.5
cat("Total missing:", sum(is.na(df_items)), "dari",
nrow(df_items) * ncol(df_items), "sel (",
round(mean(is.na(df_items)) * 100, 2), "%)\n")
## Total missing: 96 dari 9400 sel ( 1.02 %)
# Visualisasi pola missing
missing_pct_per_row <- rowMeans(is.na(df_items)) * 100
cat("Baris dengan ANY missing:", sum(rowMeans(is.na(df_items)) > 0), "\n")
## Baris dengan ANY missing: 25
cat("Baris dengan > 10% missing:", sum(missing_pct_per_row > 10), "\n")
## Baris dengan > 10% missing: 5
# --- 2.2 Multiple Imputation dengan MICE -------------------------------------
# Metode: pmm (predictive mean matching) — cocok untuk data Likert ordinal
# m = 5 → buat 5 dataset imputasi, lalu pool hasilnya
cat("\nMelakukan Multiple Imputation (MICE)...\n")
##
## Melakukan Multiple Imputation (MICE)...
set.seed(42)
imp_result <- mice(
data = df_items,
m = 5, # 5 dataset imputasi
method = "pmm", # predictive mean matching
maxit = 10, # iterasi maksimal
printFlag = FALSE
)
# Ambil dataset komplit pertama (untuk analisis awal)
# Untuk analisis final sebaiknya gunakan pool() dari semua imputasi
df_clean <- complete(imp_result, action = 1)
# Bulatkan ke integer (karena Likert adalah ordinal)
df_clean <- df_clean %>% mutate(across(everything(), ~ round(.x)))
cat("\nNormalisasi JP1-JP4 dari skala 1-10 ke 1-5...\n")
##
## Normalisasi JP1-JP4 dari skala 1-10 ke 1-5...
jp_items <- c("JP1", "JP2", "JP3", "JP4")
df_clean[, jp_items] <- lapply(df_clean[, jp_items], function(x) {
x_norm <- ceiling((x / 10) * 5)
x_norm <- pmax(1, pmin(5, x_norm)) # clamp ke [1,5]
x_norm
})
cat("Range JP setelah normalisasi:\n")
## Range JP setelah normalisasi:
print(sapply(df_clean[, jp_items], range))
## JP1 JP2 JP3 JP4
## [1,] 1 1 1 1
## [2,] 5 5 5 5
# Verifikasi: tidak ada missing lagi
cat("Missing setelah imputasi:", sum(is.na(df_clean)), "\n")
## Missing setelah imputasi: 0
# --- 2.3 Deteksi outlier multivariat (Mahalanobis distance) ------------------
cat("\nDeteksi outlier dengan Mahalanobis distance...\n")
##
## Deteksi outlier dengan Mahalanobis distance...
mah_dist <- mahalanobis(
x = df_clean,
center = colMeans(df_clean),
cov = cov(df_clean)
)
# Threshold: chi-square p < 0.001 dengan df = jumlah variabel
df_mah <- ncol(df_clean)
threshold <- qchisq(0.999, df = df_mah)
outlier_idx <- which(mah_dist > threshold)
cat("Jumlah outlier multivariat:", length(outlier_idx), "\n")
## Jumlah outlier multivariat: 6
if (length(outlier_idx) > 0) {
cat("Indeks baris outlier:", outlier_idx, "\n")
}
## Indeks baris outlier: 6 66 81 110 166 190
# Simpan data bersih
addWorksheet(wb, "Data_Clean")
writeData(wb, "Data_Clean", df_clean)
# =============================================================================
# TAHAP 3: UJI ASUMSI
# =============================================================================
cat("\n========== TAHAP 3: UJI ASUMSI ==========\n")
##
## ========== TAHAP 3: UJI ASUMSI ==========
# --- 3.1 Uji Normalitas (tanpa bergantung versi MVN) ---
cat("\n--- Uji Normalitas Multivariat ---\n")
##
## --- Uji Normalitas Multivariat ---
wellbeing_data <- df_clean[, indicators$WELLBEING]
# Shapiro-Wilk per indikator (univariat)
cat("\nHasil Shapiro-Wilk per indikator WELLBEING:\n")
##
## Hasil Shapiro-Wilk per indikator WELLBEING:
sw_results <- sapply(wellbeing_data, function(x) {
test <- shapiro.test(x)
c(W = round(test$statistic, 4), p.value = round(test$p.value, 4))
})
print(t(sw_results))
## W.W p.value
## FW1 0.7907 0
## FW2 0.8403 0
## FW3 0.7984 0
## FW4 0.8872 0
## FW5 0.7692 0
## FW6 0.7881 0
## FW7 0.8301 0
## FW8 0.8191 0
## FW9 0.8318 0
## FW10 0.8607 0
# Mardia's test manual via psych package (lebih stabil)
if (!requireNamespace("psych", quietly = TRUE)) install.packages("psych")
library(psych)
mardia_result <- mardia(wellbeing_data, plot = FALSE)
cat("\nMardia's Skewness:\n")
##
## Mardia's Skewness:
cat(sprintf(" b1p = %.4f, chi2 = %.4f, p = %.4f\n",
mardia_result$b1p,
mardia_result$chi.skew,
mardia_result$p.skew))
cat("\nMardia's Kurtosis:\n")
##
## Mardia's Kurtosis:
cat(sprintf(" b2p = %.4f, z = %.4f, p = %.4f\n",
mardia_result$b2p,
mardia_result$kurtosis,
mardia_result$p.kurt))
## b2p = 162.4572, z = 19.3790, p = 0.0000
# Interpretasi otomatis
if (mardia_result$p.skew < 0.05 | mardia_result$p.kurt < 0.05) {
cat("\n⚠ Normalitas multivariat TIDAK terpenuhi → Estimator MLR sudah tepat digunakan.\n")
} else {
cat("\n✓ Normalitas multivariat terpenuhi (p > 0.05).\n")
}
##
## ⚠ Normalitas multivariat TIDAK terpenuhi → Estimator MLR sudah tepat digunakan.
cat("\nCatatan: Script ini menggunakan estimator MLR (robust ML)\n")
##
## Catatan: Script ini menggunakan estimator MLR (robust ML)
cat("yang tidak mengasumsikan normalitas multivariat sempurna.\n")
## yang tidak mengasumsikan normalitas multivariat sempurna.
# --- 3.2 Heatmap korelasi antar konstruk ------------------------------------
cat("\n--- Korelasi Antar Konstruk ---\n")
##
## --- Korelasi Antar Konstruk ---
# Hitung rata-rata per konstruk sebagai skor komposit sementara
composite_scores <- data.frame(
JOB_DEMANDS = rowMeans(df_clean[, indicators$JOB_DEMANDS]),
JOB_RESOURCES = rowMeans(df_clean[, indicators$JOB_RESOURCES]),
WELLBEING = rowMeans(df_clean[, indicators$WELLBEING]),
WORK_ENGAGE = rowMeans(df_clean[, indicators$WORK_ENGAGE]),
ITL = rowMeans(df_clean[, indicators$ITL]),
JOB_PERF = rowMeans(df_clean[, indicators$JOB_PERF])
)
cor_matrix <- cor(composite_scores, use = "pairwise.complete.obs")
cat("Matriks korelasi antar konstruk:\n")
## Matriks korelasi antar konstruk:
print(round(cor_matrix, 3))
## JOB_DEMANDS JOB_RESOURCES WELLBEING WORK_ENGAGE ITL JOB_PERF
## JOB_DEMANDS 1.000 0.027 -0.026 -0.028 0.060 0.073
## JOB_RESOURCES 0.027 1.000 0.190 0.143 -0.104 0.076
## WELLBEING -0.026 0.190 1.000 0.025 -0.116 0.011
## WORK_ENGAGE -0.028 0.143 0.025 1.000 -0.129 0.026
## ITL 0.060 -0.104 -0.116 -0.129 1.000 0.080
## JOB_PERF 0.073 0.076 0.011 0.026 0.080 1.000
# Plot heatmap korelasi
png("plot_korelasi_konstruk.png", width = 800, height = 700, res = 120)
corrplot(
cor_matrix,
method = "color",
type = "upper",
addCoef.col = "black",
number.cex = 0.8,
tl.cex = 0.9,
tl.col = "black",
col = colorRampPalette(c("#3B8BD4", "white", "#D85A30"))(200),
title = "Korelasi Antar Konstruk",
mar = c(0, 0, 2, 0)
)
dev.off()
## png
## 2
cat("Plot korelasi disimpan: plot_korelasi_konstruk.png\n")
## Plot korelasi disimpan: plot_korelasi_konstruk.png
# --- 3.3 Uji multikolinearitas (VIF) ----------------------------------------
cat("\n--- Cek Multikolinearitas ---\n")
##
## --- Cek Multikolinearitas ---
# Korelasi > 0.85 antar konstruk mengindikasikan masalah diskriminan
high_cor <- which(abs(cor_matrix) > 0.85 & cor_matrix != 1, arr.ind = TRUE)
if (nrow(high_cor) > 0) {
cat("PERHATIAN: Pasangan konstruk dengan korelasi > 0.85:\n")
print(high_cor)
} else {
cat("Tidak ada masalah multikolinearitas (semua korelasi < 0.85)\n")
}
## Tidak ada masalah multikolinearitas (semua korelasi < 0.85)
# Simpan matriks korelasi
addWorksheet(wb, "Korelasi_Konstruk")
writeData(wb, "Korelasi_Konstruk", as.data.frame(round(cor_matrix, 3)),
rowNames = TRUE)
# =============================================================================
# TAHAP 4: MEASUREMENT MODEL (CFA)
# =============================================================================
cat("\n========== TAHAP 4: MEASUREMENT MODEL (CFA) ==========\n")
##
## ========== TAHAP 4: MEASUREMENT MODEL (CFA) ==========
# --- 4.1 Spesifikasi model CFA -----------------------------------------------
model_cfa <- '
# Konstruk 1: Job Demands (item JDR11-JDR20)
JOB_DEMANDS =~ JDR11 + JDR12 + JDR13 + JDR14 + JDR15 +
JDR16 + JDR17 + JDR18 + JDR19 + JDR20
# Konstruk 2: Job Resources (item JDR21-JDR30)
JOB_RESOURCES =~ JDR21 + JDR22 + JDR23 + JDR24 + JDR25 +
JDR26 + JDR27 + JDR28 + JDR29 + JDR30
# Konstruk 3: Wellbeing (item FW1-FW10)
WELLBEING =~ FW1 + FW2 + FW3 + FW4 + FW5 +
FW6 + FW7 + FW8 + FW9 + FW10
# Konstruk 4: Work Engagement (item WMI1-WMI10)
WORK_ENGAGE =~ WMI1 + WMI2 + WMI3 + WMI4 + WMI5 +
WMI6 + WMI7 + WMI8 + WMI9 + WMI10
# Konstruk 5: Intention to Leave (item ITL1-ITL3)
ITL =~ ITL1 + ITL2 + ITL3
# Konstruk 6: Job Performance (item JP1-JP4)
JOB_PERF =~ JP1 + JP2 + JP3 + JP4
'
# --- 4.2 Estimasi CFA --------------------------------------------------------
cat("Mengestimasi model CFA...\n")
## Mengestimasi model CFA...
fit_cfa <- cfa(
model = model_cfa,
data = df_clean,
estimator = "MLR",
std.lv = TRUE, # <-- ubah ini
missing = "fiml",
control = list(iter.max = 10000)
)
cfa_converged <- lavInspect(fit_cfa, "converged")
cat("CFA converged:", cfa_converged, "\n")
## CFA converged: TRUE
if (!cfa_converged) {
stop("Model CFA tidak konvergen. Periksa spesifikasi model atau data.")
}
# --- 4.3 Evaluasi fit CFA ----------------------------------------------------
cat("\n--- Indeks Fit CFA ---\n")
##
## --- Indeks Fit CFA ---
fit_indices_cfa <- fitMeasures(fit_cfa, c(
"chisq.scaled", "df.scaled", "pvalue.scaled",
"cfi.scaled", "tli.scaled",
"rmsea.scaled", "rmsea.ci.lower.scaled", "rmsea.ci.upper.scaled",
"srmr"
))
print(round(fit_indices_cfa, 4))
## chisq.scaled df.scaled pvalue.scaled
## 1835.531 1019.000 0.000
## cfi.scaled tli.scaled rmsea.scaled
## 0.661 0.640 0.063
## rmsea.ci.lower.scaled rmsea.ci.upper.scaled srmr
## 0.059 0.068 0.083
# Evaluasi berdasarkan kriteria
cat("\n--- Evaluasi Fit ---\n")
##
## --- Evaluasi Fit ---
cfi_val <- fit_indices_cfa["cfi.scaled"]
rmsea_val <- fit_indices_cfa["rmsea.scaled"]
srmr_val <- fit_indices_cfa["srmr"]
cat("CFI :", round(cfi_val, 3),
ifelse(cfi_val >= 0.95, "[EXCELLENT]",
ifelse(cfi_val >= 0.90, "[ACCEPTABLE]", "[POOR - perlu revisi]")), "\n")
## CFI : 0.661 [POOR - perlu revisi]
cat("RMSEA:", round(rmsea_val, 3),
ifelse(rmsea_val <= 0.06, "[EXCELLENT]",
ifelse(rmsea_val <= 0.08, "[ACCEPTABLE]", "[POOR - perlu revisi]")), "\n")
## RMSEA: 0.063 [ACCEPTABLE]
cat("SRMR :", round(srmr_val, 3),
ifelse(srmr_val <= 0.08, "[ACCEPTABLE]", "[POOR - perlu revisi]"), "\n")
## SRMR : 0.083 [POOR - perlu revisi]
# --- 4.4 Factor loadings (standardized) -------------------------------------
cat("\n--- Factor Loadings (Standardized) ---\n")
##
## --- Factor Loadings (Standardized) ---
params_cfa <- parameterEstimates(fit_cfa, standardized = TRUE) %>%
filter(op == "=~") %>%
select(Konstruk = lhs, Indikator = rhs,
Loading_Std = std.all, SE = se,
z = z, pvalue = pvalue) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)),
Signifikan = ifelse(pvalue < 0.05, "Ya", "Tidak"),
Evaluasi = ifelse(abs(Loading_Std) >= 0.50, "OK (>= 0.5)",
"Rendah (< 0.5)"))
print(params_cfa)
## Konstruk Indikator Loading_Std SE z pvalue Signifikan
## 1 JOB_DEMANDS JDR11 0.496 0.099 4.848 0.000 Ya
## 2 JOB_DEMANDS JDR12 0.057 0.103 0.634 0.526 Tidak
## 3 JOB_DEMANDS JDR13 0.468 0.155 3.245 0.001 Ya
## 4 JOB_DEMANDS JDR14 0.332 0.111 3.109 0.002 Ya
## 5 JOB_DEMANDS JDR15 0.083 0.088 0.782 0.434 Tidak
## 6 JOB_DEMANDS JDR16 0.433 0.120 3.605 0.000 Ya
## 7 JOB_DEMANDS JDR17 0.442 0.142 4.009 0.000 Ya
## 8 JOB_DEMANDS JDR18 0.654 0.166 5.027 0.000 Ya
## 9 JOB_DEMANDS JDR19 0.363 0.111 3.840 0.000 Ya
## 10 JOB_DEMANDS JDR20 -0.168 0.121 -1.535 0.125 Tidak
## 11 JOB_RESOURCES JDR21 0.621 0.086 7.274 0.000 Ya
## 12 JOB_RESOURCES JDR22 0.782 0.076 10.034 0.000 Ya
## 13 JOB_RESOURCES JDR23 0.365 0.074 4.197 0.000 Ya
## 14 JOB_RESOURCES JDR24 0.582 0.075 7.797 0.000 Ya
## 15 JOB_RESOURCES JDR25 0.639 0.058 11.076 0.000 Ya
## 16 JOB_RESOURCES JDR26 0.016 0.107 0.173 0.862 Tidak
## 17 JOB_RESOURCES JDR27 0.638 0.070 8.395 0.000 Ya
## 18 JOB_RESOURCES JDR28 0.381 0.112 3.908 0.000 Ya
## 19 JOB_RESOURCES JDR29 0.215 0.086 2.569 0.010 Ya
## 20 JOB_RESOURCES JDR30 0.445 0.076 4.858 0.000 Ya
## 21 WELLBEING FW1 0.650 0.105 6.308 0.000 Ya
## 22 WELLBEING FW2 0.355 0.104 3.777 0.000 Ya
## 23 WELLBEING FW3 0.302 0.101 3.066 0.002 Ya
## 24 WELLBEING FW4 -0.547 0.118 -6.049 0.000 Ya
## 25 WELLBEING FW5 -0.777 0.115 -8.240 0.000 Ya
## 26 WELLBEING FW6 -0.583 0.108 -6.756 0.000 Ya
## 27 WELLBEING FW7 0.259 0.134 2.537 0.011 Ya
## 28 WELLBEING FW8 0.432 0.131 4.299 0.000 Ya
## 29 WELLBEING FW9 0.317 0.092 3.416 0.001 Ya
## 30 WELLBEING FW10 0.068 0.090 0.873 0.383 Tidak
## 31 WORK_ENGAGE WMI1 0.026 0.057 0.361 0.718 Tidak
## 32 WORK_ENGAGE WMI2 -0.386 0.098 -3.834 0.000 Ya
## 33 WORK_ENGAGE WMI3 -0.082 0.098 -1.212 0.226 Tidak
## 34 WORK_ENGAGE WMI4 -0.699 0.090 -7.017 0.000 Ya
## 35 WORK_ENGAGE WMI5 -0.784 0.094 -6.705 0.000 Ya
## 36 WORK_ENGAGE WMI6 -0.241 0.103 -1.908 0.056 Tidak
## 37 WORK_ENGAGE WMI7 -0.848 0.091 -7.953 0.000 Ya
## 38 WORK_ENGAGE WMI8 -0.224 0.105 -2.008 0.045 Ya
## 39 WORK_ENGAGE WMI9 -0.763 0.082 -7.905 0.000 Ya
## 40 WORK_ENGAGE WMI10 -0.656 0.091 -5.666 0.000 Ya
## 41 ITL ITL1 0.714 0.087 10.629 0.000 Ya
## 42 ITL ITL2 0.949 0.070 14.664 0.000 Ya
## 43 ITL ITL3 0.796 0.072 13.898 0.000 Ya
## 44 JOB_PERF JP1 0.270 0.067 3.118 0.002 Ya
## 45 JOB_PERF JP2 0.331 0.080 3.304 0.001 Ya
## 46 JOB_PERF JP3 0.794 0.065 9.411 0.000 Ya
## 47 JOB_PERF JP4 0.899 0.072 10.461 0.000 Ya
## Evaluasi
## 1 Rendah (< 0.5)
## 2 Rendah (< 0.5)
## 3 Rendah (< 0.5)
## 4 Rendah (< 0.5)
## 5 Rendah (< 0.5)
## 6 Rendah (< 0.5)
## 7 Rendah (< 0.5)
## 8 OK (>= 0.5)
## 9 Rendah (< 0.5)
## 10 Rendah (< 0.5)
## 11 OK (>= 0.5)
## 12 OK (>= 0.5)
## 13 Rendah (< 0.5)
## 14 OK (>= 0.5)
## 15 OK (>= 0.5)
## 16 Rendah (< 0.5)
## 17 OK (>= 0.5)
## 18 Rendah (< 0.5)
## 19 Rendah (< 0.5)
## 20 Rendah (< 0.5)
## 21 OK (>= 0.5)
## 22 Rendah (< 0.5)
## 23 Rendah (< 0.5)
## 24 OK (>= 0.5)
## 25 OK (>= 0.5)
## 26 OK (>= 0.5)
## 27 Rendah (< 0.5)
## 28 Rendah (< 0.5)
## 29 Rendah (< 0.5)
## 30 Rendah (< 0.5)
## 31 Rendah (< 0.5)
## 32 Rendah (< 0.5)
## 33 Rendah (< 0.5)
## 34 OK (>= 0.5)
## 35 OK (>= 0.5)
## 36 Rendah (< 0.5)
## 37 OK (>= 0.5)
## 38 Rendah (< 0.5)
## 39 OK (>= 0.5)
## 40 OK (>= 0.5)
## 41 OK (>= 0.5)
## 42 OK (>= 0.5)
## 43 OK (>= 0.5)
## 44 Rendah (< 0.5)
## 45 Rendah (< 0.5)
## 46 OK (>= 0.5)
## 47 OK (>= 0.5)
# Tandai item dengan loading rendah
low_loading <- params_cfa %>% filter(abs(Loading_Std) < 0.50)
if (nrow(low_loading) > 0) {
cat("\nPERHATIAN - Item dengan loading < 0.50 (pertimbangkan dihapus):\n")
print(low_loading %>% select(Konstruk, Indikator, Loading_Std))
}
##
## PERHATIAN - Item dengan loading < 0.50 (pertimbangkan dihapus):
## Konstruk Indikator Loading_Std
## 1 JOB_DEMANDS JDR11 0.496
## 2 JOB_DEMANDS JDR12 0.057
## 3 JOB_DEMANDS JDR13 0.468
## 4 JOB_DEMANDS JDR14 0.332
## 5 JOB_DEMANDS JDR15 0.083
## 6 JOB_DEMANDS JDR16 0.433
## 7 JOB_DEMANDS JDR17 0.442
## 8 JOB_DEMANDS JDR19 0.363
## 9 JOB_DEMANDS JDR20 -0.168
## 10 JOB_RESOURCES JDR23 0.365
## 11 JOB_RESOURCES JDR26 0.016
## 12 JOB_RESOURCES JDR28 0.381
## 13 JOB_RESOURCES JDR29 0.215
## 14 JOB_RESOURCES JDR30 0.445
## 15 WELLBEING FW2 0.355
## 16 WELLBEING FW3 0.302
## 17 WELLBEING FW7 0.259
## 18 WELLBEING FW8 0.432
## 19 WELLBEING FW9 0.317
## 20 WELLBEING FW10 0.068
## 21 WORK_ENGAGE WMI1 0.026
## 22 WORK_ENGAGE WMI2 -0.386
## 23 WORK_ENGAGE WMI3 -0.082
## 24 WORK_ENGAGE WMI6 -0.241
## 25 WORK_ENGAGE WMI8 -0.224
## 26 JOB_PERF JP1 0.270
## 27 JOB_PERF JP2 0.331
# --- 4.5 Reliabilitas & validitas konstruk ----------------------------------
cat("\n--- Composite Reliability (CR) & Average Variance Extracted (AVE) ---\n")
##
## --- Composite Reliability (CR) & Average Variance Extracted (AVE) ---
# Hitung CR dan AVE secara manual dari factor loadings
hitung_cr_ave <- function(loadings) {
lam <- loadings # factor loadings
err <- 1 - lam^2 # error variance
cr <- sum(lam)^2 / (sum(lam)^2 + sum(err))
ave <- sum(lam^2) / (sum(lam^2) + sum(err))
return(c(CR = round(cr, 3), AVE = round(ave, 3)))
}
std_loads <- parameterEstimates(fit_cfa, standardized = TRUE) %>%
filter(op == "=~") %>%
select(lhs, rhs, std.all)
reliability_table <- std_loads %>%
group_by(Konstruk = lhs) %>%
summarise(
N_Item = n(),
CR = hitung_cr_ave(std.all)["CR"],
AVE = hitung_cr_ave(std.all)["AVE"],
.groups = "drop"
) %>%
mutate(
CR_OK = ifelse(CR >= 0.70, "OK (>= 0.7)", "Rendah"),
AVE_OK = ifelse(AVE >= 0.50, "OK (>= 0.5)", "Rendah")
)
print(reliability_table)
## # A tibble: 6 × 6
## Konstruk N_Item CR AVE CR_OK AVE_OK
## <chr> <int> <dbl> <dbl> <chr> <chr>
## 1 ITL 3 0.864 0.681 OK (>= 0.7) OK (>= 0.5)
## 2 JOB_DEMANDS 10 0.542 0.156 Rendah Rendah
## 3 JOB_PERF 4 0.689 0.405 Rendah Rendah
## 4 JOB_RESOURCES 10 0.75 0.267 OK (>= 0.7) Rendah
## 5 WELLBEING 10 0.028 0.224 Rendah Rendah
## 6 WORK_ENGAGE 10 0.759 0.31 OK (>= 0.7) Rendah
# --- 4.6 Discriminant validity (HTMT / sqrt AVE) ----------------------------
cat("\n--- Discriminant Validity (sqrt AVE vs korelasi konstruk) ---\n")
##
## --- Discriminant Validity (sqrt AVE vs korelasi konstruk) ---
# Dapatkan korelasi konstruk dari CFA
cor_latent <- lavInspect(fit_cfa, what = "cor.lv")
sqrt_ave <- sqrt(reliability_table$AVE)
names(sqrt_ave) <- reliability_table$Konstruk
cat("sqrt(AVE) per konstruk:\n")
## sqrt(AVE) per konstruk:
print(round(sqrt_ave, 3))
## ITL JOB_DEMANDS JOB_PERF JOB_RESOURCES WELLBEING
## 0.825 0.395 0.636 0.517 0.473
## WORK_ENGAGE
## 0.557
cat("\nKorelasi antar konstruk laten:\n")
##
## Korelasi antar konstruk laten:
print(round(cor_latent, 3))
## JOB_DE JOB_RE WELLBE WORK_E ITL JOB_PE
## JOB_DEMANDS 1.000
## JOB_RESOURCES -0.133 1.000
## WELLBEING -0.354 0.572 1.000
## WORK_ENGAGE -0.003 -0.122 -0.103 1.000
## ITL 0.115 -0.051 -0.171 0.169 1.000
## JOB_PERF 0.102 0.106 0.099 -0.057 0.129 1.000
cat("\nDiscriminant validity terpenuhi jika sqrt(AVE) > korelasi dengan konstruk lain\n")
##
## Discriminant validity terpenuhi jika sqrt(AVE) > korelasi dengan konstruk lain
# --- 4.7 Modification indices (jika fit kurang baik) -------------------------
cat("\n--- Modification Indices (top 10) ---\n")
##
## --- Modification Indices (top 10) ---
mi <- modificationIndices(fit_cfa, sort. = TRUE, maximum.number = 10)
print(mi[, c("lhs", "op", "rhs", "mi", "epc")])
## lhs op rhs mi epc
## 1479 JP1 ~~ JP2 52.336 0.292
## 215 JOB_RESOURCES =~ JDR20 49.612 0.605
## 1276 FW7 ~~ FW9 40.084 0.545
## 782 JDR20 ~~ JDR21 39.015 0.401
## 1420 WMI6 ~~ WMI8 34.267 0.302
## 819 JDR21 ~~ JDR22 29.805 0.265
## 1449 WMI9 ~~ WMI10 29.784 0.154
## 1484 JP3 ~~ JP4 23.374 1.170
## 252 WELLBEING =~ JDR20 21.168 0.438
## 1050 JDR28 ~~ JDR29 20.874 0.355
# Simpan hasil CFA ke Excel
addWorksheet(wb, "CFA_Factor_Loadings")
writeData(wb, "CFA_Factor_Loadings", params_cfa)
addWorksheet(wb, "CFA_CR_AVE")
writeData(wb, "CFA_CR_AVE", reliability_table)
addWorksheet(wb, "CFA_Fit_Indices")
fit_df <- as.data.frame(t(round(fit_indices_cfa, 4)))
writeData(wb, "CFA_Fit_Indices", fit_df)
# =============================================================================
# TAHAP 5: STRUCTURAL MODEL (SEM PENUH)
# =============================================================================
cat("\n========== TAHAP 5: STRUCTURAL MODEL (SEM PENUH) ==========\n")
##
## ========== TAHAP 5: STRUCTURAL MODEL (SEM PENUH) ==========
# --- 5.1 Spesifikasi structural model ----------------------------------------
#
# Hipotesis:
# H1: Job Demands (-) → Wellbeing
# H2: Job Resources (+) → Wellbeing
# H3: Work Engagement (+) → Wellbeing
# H4: Wellbeing (-) → ITL
# H5: Job Demands (+) → ITL (efek langsung)
# H6: Wellbeing (+) → Job Perf
# H7: Work Engagement (+) → Job Perf
# H8: Mediasi Wellbeing: JOB_DEMANDS → WELLBEING → ITL
#
model_sem <- '
# ---- MEASUREMENT MODEL ----
JOB_DEMANDS =~ JDR11 + JDR12 + JDR13 + JDR14 + JDR15 +
JDR16 + JDR17 + JDR18 + JDR19 + JDR20
JOB_RESOURCES =~ JDR21 + JDR22 + JDR23 + JDR24 + JDR25 +
JDR26 + JDR27 + JDR28 + JDR29 + JDR30
WELLBEING =~ FW1 + FW2 + FW3 + FW4 + FW5 +
FW6 + FW7 + FW8 + FW9 + FW10
WORK_ENGAGE =~ WMI1 + WMI2 + WMI3 + WMI4 + WMI5 +
WMI6 + WMI7 + WMI8 + WMI9 + WMI10
ITL =~ ITL1 + ITL2 + ITL3
JOB_PERF =~ JP1 + JP2 + JP3 + JP4
# ---- STRUCTURAL PATHS (Hipotesis) ----
# H1 & H2 & H3: Prediktor Wellbeing
WELLBEING ~ a1*JOB_DEMANDS + a2*JOB_RESOURCES + a3*WORK_ENGAGE
# H4 & H5: Prediktor ITL
ITL ~ b1*WELLBEING + b2*JOB_DEMANDS
# H6 & H7: Prediktor Job Performance
JOB_PERF ~ c1*WELLBEING + c2*WORK_ENGAGE
# ---- EFEK TIDAK LANGSUNG (Mediasi) ----
# Indirect effect: JOB_DEMANDS → WELLBEING → ITL
ind_dem_itl := a1 * b1
# Indirect effect: JOB_RESOURCES → WELLBEING → ITL
ind_res_itl := a2 * b1
# Total effect: JOB_DEMANDS → ITL
tot_dem_itl := a1 * b1 + b2
'
# --- 5.2 Estimasi SEM --------------------------------------------------------
cat("Mengestimasi model SEM...\n")
## Mengestimasi model SEM...
fit_sem <- sem(
model = model_sem,
data = df_clean,
estimator = "MLR", # Robust untuk non-normalitas
std.lv = FALSE,
missing = "fiml"
# Tidak pakai se="bootstrap" di sini
)
cat("Estimasi MLR selesai.\n")
## Estimasi MLR selesai.
sem_converged <- lavInspect(fit_sem, "converged")
cat("SEM converged:", sem_converged, "\n")
## SEM converged: TRUE
if (!sem_converged) {
# --- DIAGNOSIS otomatis ---
cat("\n========== DIAGNOSIS KONVERGENSI ==========\n")
cat("N observasi:", nrow(df_clean), "\n")
n_par <- tryCatch(
lavInspect(sem(model_sem, data=df_clean, estimator="ML", do.fit=FALSE), "npar"),
error = function(e) NA
)
cat("N parameter model:", n_par, "\n")
if (!is.na(n_par)) cat("Rasio N/par:", round(nrow(df_clean)/n_par, 1),
"(disarankan >= 5)\n")
var_check <- sapply(df_clean[, all_items], var, na.rm=TRUE)
low_var <- var_check[var_check < 0.1]
if (length(low_var) > 0) {
cat("\nItem dengan varians sangat rendah (< 0.1) — kemungkinan penyebab:\n")
print(round(low_var, 4))
} else {
cat("\nSemua item memiliki varians memadai (>= 0.1)\n")
}
cor_items <- cor(df_clean[, all_items], use="pairwise.complete.obs")
diag(cor_items) <- NA
high_cor_idx <- which(abs(cor_items) > 0.95, arr.ind=TRUE)
if (nrow(high_cor_idx) > 0) {
cat("\nPasangan item dengan korelasi > 0.95 (multikolinearitas):\n")
print(high_cor_idx)
} else {
cat("Tidak ada multikolinearitas ekstrem (korelasi < 0.95)\n")
}
cat("===========================================\n\n")
# Coba ML biasa
warning("Model SEM tidak konvergen dengan MLR. Mencoba estimator ML...")
fit_sem <- sem(
model = model_sem,
data = df_clean,
estimator = "ML",
std.lv = FALSE,
missing = "fiml",
control = list(iter.max = 10000)
)
sem_converged <- lavInspect(fit_sem, "converged")
cat("SEM converged (ML):", sem_converged, "\n")
if (!sem_converged) {
warning("Model SEM tetap tidak konvergen. Output mungkin tidak reliable, tapi analisis dilanjutkan.")
}
}
# --- 5.3 Evaluasi fit SEM ----------------------------------------------------
cat("\n--- Indeks Fit SEM ---\n")
##
## --- Indeks Fit SEM ---
fit_indices_sem <- tryCatch(
fitMeasures(fit_sem, c(
"chisq.scaled", "df.scaled", "pvalue.scaled",
"cfi.scaled", "tli.scaled",
"rmsea.scaled", "rmsea.ci.lower.scaled", "rmsea.ci.upper.scaled",
"srmr", "aic", "bic"
)),
error = function(e) {
cat("⚠ fitMeasures tidak tersedia (model belum konvergen):", conditionMessage(e), "\n")
NULL
}
)
if (!is.null(fit_indices_sem)) print(round(fit_indices_sem, 4))
## chisq.scaled df.scaled pvalue.scaled
## 1841.853 1023.000 0.000
## cfi.scaled tli.scaled rmsea.scaled
## 0.660 0.641 0.063
## rmsea.ci.lower.scaled rmsea.ci.upper.scaled srmr
## 0.059 0.068 0.084
## aic bic
## 25512.275 26013.619
# Evaluasi
cat("\n--- Evaluasi Fit Model SEM ---\n")
##
## --- Evaluasi Fit Model SEM ---
if (!is.null(fit_indices_sem)) {
for (idx in c("cfi.scaled", "tli.scaled", "rmsea.scaled", "srmr")) {
val <- fit_indices_sem[idx]
label <- switch(idx,
"cfi.scaled" = paste("CFI:", round(val, 3),
ifelse(val >= 0.95, "[EXCELLENT]",
ifelse(val >= 0.90, "[ACCEPTABLE]", "[POOR]"))),
"tli.scaled" = paste("TLI:", round(val, 3),
ifelse(val >= 0.95, "[EXCELLENT]",
ifelse(val >= 0.90, "[ACCEPTABLE]", "[POOR]"))),
"rmsea.scaled"= paste("RMSEA:", round(val, 3),
ifelse(val <= 0.06, "[EXCELLENT]",
ifelse(val <= 0.08, "[ACCEPTABLE]", "[POOR]"))),
"srmr" = paste("SRMR:", round(val, 3),
ifelse(val <= 0.08, "[ACCEPTABLE]", "[POOR]"))
)
cat(label, "\n")
}
} else {
cat("⚠ Evaluasi fit dilewati karena model tidak konvergen.\n")
}
## CFI: 0.66 [POOR]
## TLI: 0.641 [POOR]
## RMSEA: 0.063 [ACCEPTABLE]
## SRMR: 0.084 [POOR]
# --- 5.4 Path coefficients (standardized) ------------------------------------
cat("\n--- Path Coefficients Struktural (Standardized) ---\n")
##
## --- Path Coefficients Struktural (Standardized) ---
cat("\n--- Path Coefficients Struktural (Standardized, MLR) ---\n")
##
## --- Path Coefficients Struktural (Standardized, MLR) ---
params_sem <- parameterEstimates(
fit_sem,
standardized = TRUE
) %>%
filter(op == "~" | op == ":=") %>%
select(
Dari = lhs,
Arah = op,
Ke = rhs,
Beta_Std = std.all,
SE = se,
z = z,
pvalue = pvalue,
CI_Low = ci.lower,
CI_High = ci.upper
) %>%
mutate(
across(where(is.numeric), ~ round(.x, 3)),
Signifikan = case_when(
pvalue < 0.001 ~ "*** p<.001",
pvalue < 0.01 ~ "** p<.01",
pvalue < 0.05 ~ "* p<.05",
TRUE ~ "ns"
)
)
print(params_sem)
## Dari Arah Ke Beta_Std SE z pvalue CI_Low CI_High
## 1 WELLBEING ~ JOB_DEMANDS -0.278 0.212 -1.764 0.078 -0.791 0.042
## 2 WELLBEING ~ JOB_RESOURCES 0.530 0.147 3.844 0.000 0.276 0.851
## 3 WELLBEING ~ WORK_ENGAGE -0.042 5.467 -0.261 0.794 -12.141 9.291
## 4 ITL ~ WELLBEING -0.149 0.176 -1.186 0.236 -0.552 0.136
## 5 ITL ~ JOB_DEMANDS 0.051 0.237 0.407 0.684 -0.368 0.561
## 6 JOB_PERF ~ WELLBEING 0.091 0.032 0.906 0.365 -0.033 0.090
## 7 JOB_PERF ~ WORK_ENGAGE -0.073 2.440 -0.315 0.753 -5.549 4.014
## 8 ind_dem_itl := a1*b1 0.042 0.085 0.921 0.357 -0.088 0.244
## 9 ind_res_itl := a2*b1 -0.079 0.093 -1.265 0.206 -0.299 0.064
## 10 tot_dem_itl := a1*b1+b2 0.093 0.190 0.920 0.358 -0.197 0.546
## Signifikan
## 1 ns
## 2 *** p<.001
## 3 ns
## 4 ns
## 5 ns
## 6 ns
## 7 ns
## 8 ns
## 9 ns
## 10 ns
# --- 5.5 R-squared per konstruk endogen -------------------------------------
cat("\n--- R-squared (Varians yang Dijelaskan) ---\n")
##
## --- R-squared (Varians yang Dijelaskan) ---
r2 <- tryCatch(lavInspect(fit_sem, what = "r2"), error = function(e) NULL)
if (!is.null(r2)) {
r2_df <- data.frame(
Konstruk = names(r2),
R_squared = round(r2, 3),
Pct = paste0(round(r2 * 100, 1), "%")
)
print(r2_df)
} else {
cat("⚠ R-squared tidak tersedia karena model tidak konvergen.\n")
r2_df <- data.frame(Konstruk = character(), R_squared = numeric(), Pct = character())
}
## Konstruk R_squared Pct
## JDR11 JDR11 0.257 25.7%
## JDR12 JDR12 0.004 0.4%
## JDR13 JDR13 0.229 22.9%
## JDR14 JDR14 0.116 11.6%
## JDR15 JDR15 0.007 0.7%
## JDR16 JDR16 0.189 18.9%
## JDR17 JDR17 0.195 19.5%
## JDR18 JDR18 0.410 41%
## JDR19 JDR19 0.125 12.5%
## JDR20 JDR20 0.027 2.7%
## JDR21 JDR21 0.383 38.3%
## JDR22 JDR22 0.604 60.4%
## JDR23 JDR23 0.137 13.7%
## JDR24 JDR24 0.342 34.2%
## JDR25 JDR25 0.410 41%
## JDR26 JDR26 0.000 0%
## JDR27 JDR27 0.407 40.7%
## JDR28 JDR28 0.145 14.5%
## JDR29 JDR29 0.047 4.7%
## JDR30 JDR30 0.201 20.1%
## FW1 FW1 0.424 42.4%
## FW2 FW2 0.125 12.5%
## FW3 FW3 0.091 9.1%
## FW4 FW4 0.299 29.9%
## FW5 FW5 0.605 60.5%
## FW6 FW6 0.341 34.1%
## FW7 FW7 0.067 6.7%
## FW8 FW8 0.186 18.6%
## FW9 FW9 0.101 10.1%
## FW10 FW10 0.005 0.5%
## WMI1 WMI1 0.001 0.1%
## WMI2 WMI2 0.148 14.8%
## WMI3 WMI3 0.007 0.7%
## WMI4 WMI4 0.491 49.1%
## WMI5 WMI5 0.614 61.4%
## WMI6 WMI6 0.059 5.9%
## WMI7 WMI7 0.720 72%
## WMI8 WMI8 0.050 5%
## WMI9 WMI9 0.583 58.3%
## WMI10 WMI10 0.430 43%
## ITL1 ITL1 0.508 50.8%
## ITL2 ITL2 0.907 90.7%
## ITL3 ITL3 0.629 62.9%
## JP1 JP1 0.073 7.3%
## JP2 JP2 0.114 11.4%
## JP3 JP3 0.656 65.6%
## JP4 JP4 0.778 77.8%
## WELLBEING WELLBEING 0.404 40.4%
## ITL ITL 0.030 3%
## JOB_PERF JOB_PERF 0.015 1.5%
# --- 5.6 Efek mediasi (indirect effects) ------------------------------------
cat("\n--- Uji Mediasi (Bootstrap CI) ---\n")
##
## --- Uji Mediasi (Bootstrap CI) ---
# Estimasi ulang dengan bootstrap untuk CI mediasi
cat("Mengestimasi bootstrap SEM (1000 iterasi, harap tunggu)...\n")
## Mengestimasi bootstrap SEM (1000 iterasi, harap tunggu)...
set.seed(123)
fit_sem_boot <- tryCatch(
sem(
model = model_sem,
data = df_clean,
estimator = "MLR",
std.lv = FALSE,
missing = "fiml",
se = "bootstrap",
bootstrap = 1000
),
error = function(e) { cat("Bootstrap gagal:", conditionMessage(e), "\n"); NULL }
)
## Bootstrap gagal: lavaan->lav_options_est_ml():
## use ML estimator for bootstrap
if (is.null(fit_sem_boot) || !lavInspect(fit_sem_boot, "converged")) {
cat("⚠ Bootstrap SEM tidak konvergen, menggunakan hasil MLR biasa untuk mediasi.\n")
fit_sem_boot <- fit_sem
}
## ⚠ Bootstrap SEM tidak konvergen, menggunakan hasil MLR biasa untuk mediasi.
indirect_boot <- parameterEstimates(
fit_sem_boot,
standardized = TRUE,
boot.ci.type = "bca.simple" # Bias-corrected accelerated
) %>%
filter(op == ":=") %>%
select(
Efek = lhs,
Beta_Std = std.all,
SE = se,
CI_Low = ci.lower,
CI_High = ci.upper,
pvalue = pvalue
) %>%
mutate(
across(where(is.numeric), ~ round(.x, 3)),
Signifikan = ifelse(pvalue < 0.05, "Ya*", "Tidak"),
CI_Tidak_Nol = ifelse(CI_Low > 0 | CI_High < 0, "✓ Mediasi", "✗ Tidak"),
Interpretasi = case_when(
CI_Low > 0 | CI_High < 0 ~ "Mediasi signifikan",
TRUE ~ "Mediasi tidak signifikan"
)
)
print(indirect_boot)
## Efek Beta_Std SE CI_Low CI_High pvalue Signifikan CI_Tidak_Nol
## 1 ind_dem_itl 0.042 0.085 -0.088 0.244 0.357 Tidak ✗ Tidak
## 2 ind_res_itl -0.079 0.093 -0.299 0.064 0.206 Tidak ✗ Tidak
## 3 tot_dem_itl 0.093 0.190 -0.197 0.546 0.358 Tidak ✗ Tidak
## Interpretasi
## 1 Mediasi tidak signifikan
## 2 Mediasi tidak signifikan
## 3 Mediasi tidak signifikan
cat("\nKeterangan: Mediasi signifikan jika 95% CI tidak mencakup angka 0\n")
##
## Keterangan: Mediasi signifikan jika 95% CI tidak mencakup angka 0
# --- 5.7 Perbandingan model (opsional: model alternatif) --------------------
# Model alternatif: tanpa jalur langsung JOB_DEMANDS → ITL
model_alt <- '
JOB_DEMANDS =~ JDR11+JDR12+JDR13+JDR14+JDR15+JDR16+JDR17+JDR18+JDR19+JDR20
JOB_RESOURCES =~ JDR21+JDR22+JDR23+JDR24+JDR25+JDR26+JDR27+JDR28+JDR29+JDR30
WELLBEING =~ FW1+FW2+FW3+FW4+FW5+FW6+FW7+FW8+FW9+FW10
WORK_ENGAGE =~ WMI1+WMI2+WMI3+WMI4+WMI5+WMI6+WMI7+WMI8+WMI9+WMI10
ITL =~ ITL1+ITL2+ITL3
JOB_PERF =~ JP1+JP2+JP3+JP4
WELLBEING ~ JOB_DEMANDS + JOB_RESOURCES + WORK_ENGAGE
ITL ~ WELLBEING # tanpa jalur langsung dari JOB_DEMANDS
JOB_PERF ~ WELLBEING + WORK_ENGAGE
'
fit_alt <- sem(
model_alt,
data = df_clean,
estimator = "MLR",
missing = "fiml"
)
## Warning: lavaan->lav_lavaan_step11_estoptim():
## the optimizer warns that a solution has NOT been found!
cat("\n--- Perbandingan Model Utama vs Model Alternatif ---\n")
##
## --- Perbandingan Model Utama vs Model Alternatif ---
tryCatch({
comp <- compareFit(fit_sem, fit_alt)
summary(comp)
}, error = function(e) {
cat("⚠ Perbandingan model gagal:", conditionMessage(e), "\n")
cat("Tip: Pastikan kedua model konvergen sebelum membandingkan.\n")
})
## The following models did not converge, so they are ignored:
## fit_alt
## ################### Nested Model Comparison #########################
## Chi-Squared Test Statistic (scaled)
## Yuan-Bentler correction (Mplus variant)
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## Saturated 0 0.0
## Model 1023 25512 26014 1841.9 1841.8 1023 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ####################### Model Fit Indices ###########################
## ⚠ Perbandingan model gagal: invalid 'row.names' length
## Tip: Pastikan kedua model konvergen sebelum membandingkan.
# Simpan hasil SEM ke Excel
addWorksheet(wb, "SEM_Path_Coefficients")
writeData(wb, "SEM_Path_Coefficients", params_sem)
addWorksheet(wb, "SEM_Indirect_Bootstrap")
writeData(wb, "SEM_Indirect_Bootstrap", indirect_boot)
# =============================================================================
# TAHAP 6: INTERPRETASI & PELAPORAN + VISUALISASI LENGKAP
# =============================================================================
cat("\n========== TAHAP 6: INTERPRETASI & PELAPORAN ==========\n")
##
## ========== TAHAP 6: INTERPRETASI & PELAPORAN ==========
# Helper aman untuk fitMeasures
safe_fit <- function(fit_obj, measure) {
tryCatch(round(fitMeasures(fit_obj, measure), 3),
error = function(e) NA_real_)
}
# --- 6.1 Path diagram SEM (inline + file) ------------------------------------
cat("Membuat path diagram...\n")
## Membuat path diagram...
semPaths(
object = fit_sem,
what = "std",
whatLabels = "std",
layout = "tree2",
rotation = 2,
edge.label.cex = 0.65,
node.label.cex = 0.75,
edge.color = "black",
style = "ram",
residuals = TRUE,
intercepts = FALSE,
nCharNodes = 10,
title = TRUE,
title.cex = 1.1,
mar = c(2, 2, 2, 2),
color = list(lat = "#CECBF6", man = "#9FE1CB")
)
## Warning in qgraph::qgraph(Edgelist, labels = nLab, bidirectional = Bidir, : The
## following arguments are not documented and likely not arguments of qgraph and
## thus ignored: node.label.cex
png("plot_path_diagram_SEM.png", width = 1400, height = 1000, res = 140)
semPaths(fit_sem, what="std", whatLabels="std", layout="tree2", rotation=2,
edge.label.cex=0.65, node.label.cex=0.75, edge.color="black",
style="ram", residuals=TRUE, intercepts=FALSE, nCharNodes=10,
color=list(lat="#CECBF6", man="#9FE1CB"), mar=c(2,2,2,2))
## Warning in qgraph::qgraph(Edgelist, labels = nLab, bidirectional = Bidir, : The
## following arguments are not documented and likely not arguments of qgraph and
## thus ignored: node.label.cex
dev.off()
## png
## 2
cat("Path diagram disimpan: plot_path_diagram_SEM.png\n")
## Path diagram disimpan: plot_path_diagram_SEM.png
# --- 6.2 Plot Factor Loadings per konstruk -----------------------------------
cat("\nMembuat plot factor loadings...\n")
##
## Membuat plot factor loadings...
p_loadings <- ggplot(params_cfa,
aes(x = reorder(Indikator, Loading_Std), y = Loading_Std, fill = Konstruk)) +
geom_col(width = 0.7, color = "white", linewidth = 0.3) +
geom_hline(yintercept = 0.50, linetype = "dashed", color = "red", linewidth = 0.5) +
geom_text(aes(label = sprintf("%.2f", Loading_Std),
y = Loading_Std + 0.03), size = 2.8, hjust = 0) +
coord_flip() +
facet_wrap(~Konstruk, scales = "free_y", ncol = 2) +
scale_fill_brewer(palette = "Set2") +
scale_y_continuous(limits = c(0, 1.1)) +
labs(
title = "Factor Loadings (Standardized) per Konstruk",
subtitle = "Garis merah putus-putus = batas minimal 0.50",
x = NULL,
y = "Standardized Loading (λ)",
fill = "Konstruk"
) +
theme_minimal(base_size = 10) +
theme(
strip.text = element_text(face = "bold", size = 10),
plot.title = element_text(face = "bold", size = 13),
legend.position = "none",
panel.grid.minor = element_blank()
)
print(p_loadings)
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_col()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_text()`).
ggsave("plot_factor_loadings.png", p_loadings, width = 12, height = 9, dpi = 150)
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_col()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_text()`).
cat("Plot factor loadings disimpan: plot_factor_loadings.png\n")
## Plot factor loadings disimpan: plot_factor_loadings.png
# --- 6.3 Plot CR & AVE per konstruk ------------------------------------------
cat("\nMembuat plot CR dan AVE...\n")
##
## Membuat plot CR dan AVE...
cr_ave_long <- reliability_table %>%
select(Konstruk, CR, AVE) %>%
tidyr::pivot_longer(cols = c(CR, AVE), names_to = "Metrik", values_to = "Nilai")
cr_ave_long$Threshold <- ifelse(cr_ave_long$Metrik == "CR", 0.70, 0.50)
cr_ave_long$Status <- ifelse(cr_ave_long$Nilai >= cr_ave_long$Threshold, "Memenuhi", "Tidak")
p_crave <- ggplot(cr_ave_long,
aes(x = Konstruk, y = Nilai, fill = Status)) +
geom_col(width = 0.6, position = position_dodge(0.7), color = "white") +
geom_hline(data = data.frame(Metrik = c("CR","AVE"), Threshold = c(0.70, 0.50)),
aes(yintercept = Threshold), linetype = "dashed", color = "red",
linewidth = 0.6) +
geom_text(aes(label = sprintf("%.3f", Nilai), y = Nilai + 0.02),
position = position_dodge(0.7), size = 3) +
facet_wrap(~Metrik, ncol = 2) +
scale_fill_manual(values = c("Memenuhi" = "#1D9E75", "Tidak" = "#D85A30")) +
scale_y_continuous(limits = c(0, 1.1)) +
labs(
title = "Composite Reliability (CR) & Average Variance Extracted (AVE)",
subtitle = "CR ≥ 0.70 dan AVE ≥ 0.50 menunjukkan validitas & reliabilitas baik",
x = NULL, y = "Nilai",
fill = "Status"
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 13),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 30, hjust = 1),
panel.grid.minor = element_blank()
)
print(p_crave)
ggsave("plot_cr_ave.png", p_crave, width = 10, height = 5, dpi = 150)
cat("Plot CR & AVE disimpan: plot_cr_ave.png\n")
## Plot CR & AVE disimpan: plot_cr_ave.png
# --- 6.4 Plot path coefficients (struktural) ---------------------------------
cat("\nMembuat plot koefisien jalur...\n")
##
## Membuat plot koefisien jalur...
path_plot_data <- params_sem %>%
filter(Arah == "~") %>%
mutate(
Jalur = paste(Dari, "→", Ke),
Warna = ifelse(Beta_Std > 0, "Positif", "Negatif"),
Sig_Sym = case_when(
pvalue < 0.001 ~ "***",
pvalue < 0.01 ~ "**",
pvalue < 0.05 ~ "*",
TRUE ~ "ns"
)
)
p_paths <- ggplot(path_plot_data,
aes(x = reorder(Jalur, Beta_Std), y = Beta_Std, fill = Warna)) +
geom_col(width = 0.6, color = "white", linewidth = 0.3) +
geom_errorbar(aes(ymin = CI_Low, ymax = CI_High),
width = 0.25, linewidth = 0.6) +
geom_text(aes(label = Sig_Sym,
y = ifelse(Beta_Std >= 0, CI_High + 0.04, CI_Low - 0.04)),
size = 4, fontface = "bold") +
geom_hline(yintercept = 0, linewidth = 0.5, color = "gray40") +
coord_flip() +
scale_fill_manual(values = c("Positif" = "#1D9E75", "Negatif" = "#D85A30")) +
labs(
title = "Structural Path Coefficients (Standardized β)",
subtitle = "Error bars = 95% CI | *** p<.001 ** p<.01 * p<.05 ns = tidak signifikan",
x = NULL,
y = "Standardized β",
fill = "Arah Pengaruh"
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 13),
legend.position = "bottom",
panel.grid.minor = element_blank()
)
print(p_paths)
ggsave("plot_path_coefficients.png", p_paths, width = 9, height = 6, dpi = 150)
cat("Plot path coefficients disimpan: plot_path_coefficients.png\n")
## Plot path coefficients disimpan: plot_path_coefficients.png
# --- 6.5 Plot R-squared per konstruk endogen ---------------------------------
if (nrow(r2_df) > 0) {
cat("\nMembuat plot R-squared...\n")
p_r2 <- ggplot(r2_df, aes(x = reorder(Konstruk, R_squared), y = R_squared)) +
geom_col(fill = "#5B85C8", width = 0.6, color = "white") +
geom_text(aes(label = Pct, y = R_squared + 0.02),
size = 4, fontface = "bold") +
geom_hline(yintercept = 0.10, linetype = "dashed",
color = "#D85A30", linewidth = 0.5) +
coord_flip() +
scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
labs(
title = "R-squared: Varians Konstruk Endogen yang Dijelaskan",
subtitle = "Garis merah = ambang batas minimal 10%",
x = NULL,
y = "R²"
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 13),
panel.grid.minor = element_blank()
)
print(p_r2)
ggsave("plot_r_squared.png", p_r2, width = 8, height = 5, dpi = 150)
cat("Plot R-squared disimpan: plot_r_squared.png\n")
}
##
## Membuat plot R-squared...
## Plot R-squared disimpan: plot_r_squared.png
# --- 6.6 Plot perbandingan fit indices CFA vs SEM ----------------------------
cat("\nMembuat plot perbandingan fit indices...\n")
##
## Membuat plot perbandingan fit indices...
fit_comparison <- data.frame(
Model = c("CFA (Measurement)", "SEM (Structural)"),
CFI = c(safe_fit(fit_cfa, "cfi.scaled"), safe_fit(fit_sem, "cfi.scaled")),
TLI = c(safe_fit(fit_cfa, "tli.scaled"), safe_fit(fit_sem, "tli.scaled")),
RMSEA = c(safe_fit(fit_cfa, "rmsea.scaled"), safe_fit(fit_sem, "rmsea.scaled")),
SRMR = c(safe_fit(fit_cfa, "srmr"), safe_fit(fit_sem, "srmr"))
)
print(fit_comparison)
## Model CFI TLI RMSEA SRMR
## 1 CFA (Measurement) 0.661 0.640 0.063 0.083
## 2 SEM (Structural) 0.660 0.641 0.063 0.084
fit_long <- fit_comparison %>%
tidyr::pivot_longer(cols = c(CFI, TLI, RMSEA, SRMR),
names_to = "Indeks", values_to = "Nilai") %>%
mutate(
Threshold = case_when(
Indeks %in% c("CFI","TLI") ~ 0.90,
Indeks == "RMSEA" ~ 0.08,
Indeks == "SRMR" ~ 0.08
),
Status = case_when(
Indeks %in% c("CFI","TLI") & Nilai >= 0.95 ~ "Excellent",
Indeks %in% c("CFI","TLI") & Nilai >= 0.90 ~ "Acceptable",
Indeks %in% c("RMSEA","SRMR") & Nilai <= 0.06 ~ "Excellent",
Indeks %in% c("RMSEA","SRMR") & Nilai <= 0.08 ~ "Acceptable",
TRUE ~ "Poor"
)
)
p_fit <- ggplot(fit_long,
aes(x = Model, y = Nilai, fill = Status, group = Model)) +
geom_col(width = 0.5, color = "white") +
geom_hline(aes(yintercept = Threshold),
linetype = "dashed", color = "red", linewidth = 0.5) +
geom_text(aes(label = sprintf("%.3f", Nilai), y = Nilai + 0.01),
size = 3.5, fontface = "bold") +
facet_wrap(~Indeks, scales = "free_y", nrow = 1) +
scale_fill_manual(values = c("Excellent" = "#1D9E75",
"Acceptable" = "#F0A500",
"Poor" = "#D85A30")) +
labs(
title = "Perbandingan Fit Indices: CFA vs SEM",
subtitle = "Garis merah = ambang batas | CFI/TLI ≥ 0.90 | RMSEA/SRMR ≤ 0.08",
x = NULL, y = "Nilai Indeks",
fill = "Kategori"
) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(face = "bold", size = 13),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 20, hjust = 1),
panel.grid.minor = element_blank()
)
print(p_fit)
ggsave("plot_fit_comparison.png", p_fit, width = 10, height = 5, dpi = 150)
cat("Plot fit comparison disimpan: plot_fit_comparison.png\n")
## Plot fit comparison disimpan: plot_fit_comparison.png
# --- 6.7 Plot efek mediasi ---------------------------------------------------
if (nrow(indirect_boot) > 0) {
cat("\nMembuat plot efek mediasi...\n")
p_med <- ggplot(indirect_boot,
aes(x = reorder(Efek, Beta_Std), y = Beta_Std,
color = CI_Tidak_Nol)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = CI_Low, ymax = CI_High),
width = 0.15, linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed",
color = "gray50", linewidth = 0.5) +
coord_flip() +
scale_color_manual(values = c("✓ Mediasi" = "#1D9E75",
"✗ Tidak" = "#D85A30")) +
labs(
title = "Efek Mediasi (Indirect Effects)",
subtitle = "Error bars = 95% CI | Hijau = mediasi signifikan (CI tidak mencakup 0)",
x = NULL,
y = "Standardized Indirect Effect (β)",
color = "Kesimpulan"
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 13),
legend.position = "bottom",
panel.grid.minor = element_blank()
)
print(p_med)
ggsave("plot_mediasi.png", p_med, width = 8, height = 5, dpi = 150)
cat("Plot mediasi disimpan: plot_mediasi.png\n")
}
##
## Membuat plot efek mediasi...
## Plot mediasi disimpan: plot_mediasi.png
# --- 6.8 Tabel ringkasan hipotesis -------------------------------------------
cat("\n--- Ringkasan Hasil Hipotesis ---\n")
##
## --- Ringkasan Hasil Hipotesis ---
hipotesis <- data.frame(
Hipotesis = c("H1","H2","H3","H4","H5","H6","H7","H8 (mediasi)"),
Jalur = c(
"JOB_DEMANDS → WELLBEING",
"JOB_RESOURCES → WELLBEING",
"WORK_ENGAGE → WELLBEING",
"WELLBEING → ITL",
"JOB_DEMANDS → ITL (langsung)",
"WELLBEING → JOB_PERF",
"WORK_ENGAGE → JOB_PERF",
"JOB_DEMANDS → WELLBEING → ITL"
),
Arah_Prediksi = c("-","+","+","-","+","+","+","-")
)
print(hipotesis)
## Hipotesis Jalur Arah_Prediksi
## 1 H1 JOB_DEMANDS → WELLBEING -
## 2 H2 JOB_RESOURCES → WELLBEING +
## 3 H3 WORK_ENGAGE → WELLBEING +
## 4 H4 WELLBEING → ITL -
## 5 H5 JOB_DEMANDS → ITL (langsung) +
## 6 H6 WELLBEING → JOB_PERF +
## 7 H7 WORK_ENGAGE → JOB_PERF +
## 8 H8 (mediasi) JOB_DEMANDS → WELLBEING → ITL -
# --- 6.9 Simpan semua ke Excel -----------------------------------------------
addWorksheet(wb, "Fit_Comparison")
writeData(wb, "Fit_Comparison", fit_comparison)
addWorksheet(wb, "Hipotesis_Summary")
writeData(wb, "Hipotesis_Summary", hipotesis)
if (!is.null(fit_indices_sem)) {
addWorksheet(wb, "SEM_Fit_Indices")
writeData(wb, "SEM_Fit_Indices",
as.data.frame(t(round(fit_indices_sem, 4))))
}
if (nrow(r2_df) > 0) {
addWorksheet(wb, "SEM_R_Squared")
writeData(wb, "SEM_R_Squared", r2_df)
}
saveWorkbook(wb, "SEM_Hasil_Lengkap.xlsx", overwrite = TRUE)
cat("\nSemua hasil disimpan ke: SEM_Hasil_Lengkap.xlsx\n")
##
## Semua hasil disimpan ke: SEM_Hasil_Lengkap.xlsx
# --- 6.10 Ringkasan akhir ----------------------------------------------------
cat("\n=============================================================\n")
##
## =============================================================
cat(" RINGKASAN AKHIR ANALISIS SEM\n")
## RINGKASAN AKHIR ANALISIS SEM
cat("=============================================================\n")
## =============================================================
cat("Dataset :", nrow(df_clean), "responden,", length(all_items), "item\n")
## Dataset : 200 responden, 47 item
cat("Konstruk :", length(indicators), "konstruk laten\n")
## Konstruk : 6 konstruk laten
cat("Estimator : MLR (Robust Maximum Likelihood)\n")
## Estimator : MLR (Robust Maximum Likelihood)
cat("JP1-JP4 : Dinormalisasi dari skala 1-10 ke 1-5\n")
## JP1-JP4 : Dinormalisasi dari skala 1-10 ke 1-5
cat("\nFit Indeks CFA:\n")
##
## Fit Indeks CFA:
cat(" CFI =", safe_fit(fit_cfa, "cfi.scaled"), "\n")
## CFI = 0.661
cat(" RMSEA =", safe_fit(fit_cfa, "rmsea.scaled"), "\n")
## RMSEA = 0.063
cat(" SRMR =", safe_fit(fit_cfa, "srmr"), "\n")
## SRMR = 0.083
cat("\nFit Indeks SEM:\n")
##
## Fit Indeks SEM:
cat(" CFI =", safe_fit(fit_sem, "cfi.scaled"), "\n")
## CFI = 0.66
cat(" RMSEA =", safe_fit(fit_sem, "rmsea.scaled"), "\n")
## RMSEA = 0.063
cat(" SRMR =", safe_fit(fit_sem, "srmr"), "\n")
## SRMR = 0.084
cat("\nFile output:\n")
##
## File output:
cat(" SEM_Hasil_Lengkap.xlsx\n")
## SEM_Hasil_Lengkap.xlsx
cat(" plot_korelasi_konstruk.png\n")
## plot_korelasi_konstruk.png
cat(" plot_path_diagram_SEM.png\n")
## plot_path_diagram_SEM.png
cat(" plot_factor_loadings.png\n")
## plot_factor_loadings.png
cat(" plot_cr_ave.png\n")
## plot_cr_ave.png
cat(" plot_path_coefficients.png\n")
## plot_path_coefficients.png
cat(" plot_r_squared.png\n")
## plot_r_squared.png
cat(" plot_fit_comparison.png\n")
## plot_fit_comparison.png
cat(" plot_mediasi.png\n")
## plot_mediasi.png
cat("=============================================================\n")
## =============================================================