If following approach for testing SEM and Moderation Suitable?
The main doubt is, should moderator variable be added in the structural model part, considering it is being tested in Parts 3 and 4 of the code?
If yes why, if not why?
# ==============================================================================
# PART 1: SETUP & LIBRARIES
# ==============================================================================
library(lavaan)
library(semTools)
library(readxl)
# Load Data
Fulldata <- read_excel("your_data_path.xlsx")
# --- DATA CLEANING ---
Fulldata <- na.omit(Fulldata)
# Convert Demographics to Numeric (Generic Labels)
Fulldata$C1 <- as.numeric(as.factor(Fulldata$Gender)) - 1
Fulldata$C2 <- as.numeric(as.factor(Fulldata$Occupation))
Fulldata$C3 <- as.numeric(as.factor(Fulldata$Education))
Fulldata$C4 <- as.numeric(as.factor(Fulldata$Income))
Fulldata$C5 <- as.numeric(Fulldata$Age)
# ==============================================================================
# PART 2: THE MAIN EFFECTS MODEL
# ==============================================================================
Model_Main <- '
# Measurement Model
IV_1 =~ AB1 + AB2 + AB3
IV_2 =~ CD1 + CD2 + CD3
IV_3 =~ EF1 + EF2 + EF3
IV_4 =~ GH1 + GH2 + GH3
IV_5 =~ IJ1 + IJ2 + IJ3
IV_6 =~ KL1 + KL2 + KL3
Moderator =~ MN1 + MN2 + MN3
Dependent_Var =~ OP1 + OP2 + OP3
# Structural Model
Dependent_Var ~ IV_1 + IV_2 + IV_3 + IV_4 +
IV_5 + IV_6 + Moderator +
C1 + C2 + C3 + C4 + C5
'
# Using MLM for robust estimation
fit_main <- sem(Model_Main, estimator="MLM", data=Fulldata, meanstructure=TRUE)
print("--- MAIN EFFECTS MODEL RESULTS ---")
summary(fit_main, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE)
print(fitmeasures(fit_main, c("chisq.scaled","df.scaled","pvalue.scaled","tli.robust","cfi.robust","rmsea.robust","srmr")))
# ==============================================================================
# PART 3: AUTOMATED MODERATION ANALYSIS
# ==============================================================================
# 1. VARIABLE DICTIONARY
item_dict <- list(
IV_1 = paste0("AB", 1:3),
IV_2 = paste0("CD", 1:3),
IV_3 = paste0("EF", 1:3),
IV_4 = paste0("GH", 1:3),
IV_5 = paste0("IJ", 1:3),
IV_6 = paste0("KL", 1:3),
Moderator = paste0("MN", 1:3),
Dependent_Var = paste0("OP", 1:3)
)
control_vars <- "C1 + C2 + C3 + C4 + C5"
# 2. MODERATION FUNCTION
run_moderation_robust <- function(df, iv_name, mod_name, dv_name) {
iv_items <- item_dict[[iv_name]]
mod_items <- item_dict[[mod_name]]
dv_items <- item_dict[[dv_name]]
# Calculate Mean Scores and Center for Interaction
iv_score <- rowMeans(df[, iv_items], na.rm = TRUE)
mod_score <- rowMeans(df[, mod_items], na.rm = TRUE)
iv_centered <- scale(iv_score, center=TRUE, scale=FALSE)
mod_centered <- scale(mod_score, center=TRUE, scale=FALSE)
# Add Interaction Term to df
df$Interaction_Observed <- as.numeric(iv_centered * mod_centered)
model_syntax <- paste0("
# Measurement Models
", iv_name, " =~ ", paste(iv_items, collapse=" + "), "
", mod_name, " =~ ", paste(mod_items, collapse=" + "), "
", dv_name, " =~ ", paste(dv_items, collapse=" + "), "
# Structural Regression
", dv_name, " ~ ", iv_name, " + ", mod_name, " + Interaction_Observed + ", control_vars, "
# Covariances
Interaction_Observed ~~ ", iv_name, "
Interaction_Observed ~~ ", mod_name, "
")
fit <- sem(model_syntax, data = df,
std.lv = TRUE, estimator="MLM")
return(fit)
}
# ==============================================================================
# PART 4: EXECUTION LOOP
# ==============================================================================
IV_List <- c("IV_1", "IV_2", "IV_3", "IV_4", "IV_5", "IV_6")
Mod_List <- c("Moderator")
DV_List <- c("Dependent_Var")
results_table <- data.frame()
for (mod in Mod_List) {
for (dv in DV_List) {
for (iv in IV_List) {
tryCatch({
fit <- run_moderation_robust(Fulldata, iv, mod, dv)
if (inspect(fit, "converged")) {
params <- parameterEstimates(fit, standardized = TRUE)
interaction_row <- subset(params, op == "~" & rhs == "Interaction_Observed")
if(nrow(interaction_row) > 0) {
new_row <- data.frame(
Independent_Var = iv,
Moderator = mod,
Dependent_Var = dv,
Interaction_Est = round(interaction_row$est, 3),
P_Value = round(interaction_row$pvalue, 3),
Std_Beta = round(interaction_row$std.all, 3),
Significant = ifelse(interaction_row$pvalue < 0.05, "YES", "NO")
)
results_table <- rbind(results_table, new_row)
}
}
}, error = function(e) {
cat(paste("Error in model for IV:", iv, "\n"))
})
}
}
}
# ==============================================================================
# PART 5: VIEW RESULTS
# ==============================================================================
print("--- FINAL MODERATION RESULTS ---")
print(results_table)