# ============================================================================
# PSYC 434 — Lab 10: Measurement Invariance
# self-standing script — run from top to bottom
# ============================================================================

# --- packages ---------------------------------------------------------------

library(causalworkshop)
library(psych)
library(lavaan)
library(tidyverse)

# --- generate data ----------------------------------------------------------

d <- simulate_measurement_items(n = 2000, seed = 2026)
dim(d)
names(d)

# check true factor loadings and intercepts
attr(d, "true_loadings")
attr(d, "true_intercepts_group0")
attr(d, "true_intercepts_group1")

# --- exploratory factor analysis (EFA) --------------------------------------

# select items
items <- d |> select(item_1:item_6)

# factorability
psych::KMO(items)
psych::cortest.bartlett(cor(items), n = nrow(items))

# one-factor solution
fa_1 <- psych::fa(items, nfactors = 1, fm = "ml", rotate = "none")
print(fa_1$loadings, cutoff = 0.3)

# two-factor solution (for comparison)
fa_2 <- psych::fa(items, nfactors = 2, fm = "ml", rotate = "oblimin")
print(fa_2$loadings, cutoff = 0.3)

# --- confirmatory factor analysis (CFA) ------------------------------------

# specify one-factor model
model <- "
  distress =~ item_1 + item_2 + item_3 + item_4 + item_5 + item_6
"

# fit CFA on full sample
fit_cfa <- cfa(model, data = d)
summary(fit_cfa, fit.measures = TRUE, standardized = TRUE)

# extract key fit indices
fit_indices <- fitmeasures(fit_cfa, c("cfi", "rmsea", "srmr"))
print(round(fit_indices, 3))

# --- multigroup CFA: invariance testing -------------------------------------

# step 1: configural invariance
fit_configural <- cfa(model, data = d, group = "group")
summary(fit_configural, fit.measures = TRUE)

# step 2: metric invariance (equal loadings)
fit_metric <- cfa(model, data = d, group = "group",
                  group.equal = "loadings")
summary(fit_metric, fit.measures = TRUE)

# compare configural vs metric
lavTestLRT(fit_configural, fit_metric)

# step 3: scalar invariance (equal loadings + intercepts)
fit_scalar <- cfa(model, data = d, group = "group",
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, fit.measures = TRUE)

# compare metric vs scalar
lavTestLRT(fit_metric, fit_scalar)

# --- partial scalar invariance ----------------------------------------------

# free intercepts for items 3 and 5
model_partial <- "
  distress =~ item_1 + item_2 + item_3 + item_4 + item_5 + item_6
  item_3 ~ c(i3a, i3b) * 1
  item_5 ~ c(i5a, i5b) * 1
"

fit_partial <- cfa(model_partial, data = d, group = "group",
                   group.equal = c("loadings", "intercepts"))
summary(fit_partial, fit.measures = TRUE)

# compare partial scalar vs metric
lavTestLRT(fit_metric, fit_partial)

# --- compare all models -----------------------------------------------------

models <- list(
  Configural = fit_configural,
  Metric = fit_metric,
  Scalar = fit_scalar,
  "Partial Scalar" = fit_partial
)

fit_table <- map_dfr(names(models), function(name) {
  fm <- fitmeasures(models[[name]], c("cfi", "rmsea", "srmr", "chisq", "df"))
  tibble(
    model = name,
    cfi = round(fm["cfi"], 3),
    rmsea = round(fm["rmsea"], 3),
    srmr = round(fm["srmr"], 3),
    chisq = round(fm["chisq"], 1),
    df = fm["df"]
  )
})

print(fit_table)