## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(plavaan)
library(lavaan)
data(HolzingerSwineford1939)

## -----------------------------------------------------------------------------
# Select the 9 cognitive test variables
hs_data <- HolzingerSwineford1939[, c("school", "x1", "x2", "x3", "x4", 
                                       "x5", "x6", "x7", "x8", "x9")]

# Convert to ordinal with 3 points
for (i in 2:10) {
    hs_data[[i]] <- cut(
        hs_data[[i]], 
        breaks = 3, 
        labels = FALSE,
        include.lowest = TRUE
    )
    hs_data[[i]] <- ordered(hs_data[[i]])
}

head(hs_data)

## -----------------------------------------------------------------------------
mod_base <- "
  visual =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed =~ x7 + x8 + x9
  x1 ~~ 1 * x1
  x2 ~~ 1 * x2
  x3 ~~ 1 * x3
  x4 ~~ 1 * x4
  x5 ~~ 1 * x5
  x6 ~~ 1 * x6
  x7 ~~ 1 * x7
  x8 ~~ 1 * x8
  x9 ~~ 1 * x9
"
fit_mg <- cfa(mod_base, data = hs_data, ordered = TRUE, std.lv = TRUE,
              parameterization = "theta", group = "school")
summary(fit_mg, fit.measures = TRUE)

## -----------------------------------------------------------------------------
fit_mg@optim$fx

## -----------------------------------------------------------------------------
# Strict invariance: constrain loadings, thresholds, and residual variances
fit_strict <- cfa(mod_base, data = hs_data, ordered = TRUE, std.lv = TRUE,
                  parameterization = "theta", group = "school",
                  group.equal = c("loadings", "thresholds", "residuals"))

# Score test
lavTestScore(fit_strict)

## -----------------------------------------------------------------------------
mod_un <- "
  visual =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed =~ x7 + x8 + x9
  visual ~~ c(1, NA) * visual
  textual ~~ c(1, NA) * textual
  speed ~~ c(1, NA) * speed
  visual ~ c(0, NA) * 1
  textual ~ c(0, NA) * 1
  speed ~ c(0, NA) * 1
  x1 ~~ 1 * x1
  x2 ~~ 1 * x2
  x3 ~~ 1 * x3
  x4 ~~ 1 * x4
  x5 ~~ 1 * x5
  x6 ~~ 1 * x6
  x7 ~~ 1 * x7
  x8 ~~ 1 * x8
  x9 ~~ 1 * x9
"
fit_mg_nofit <- cfa(mod_base, data = hs_data, ordered = TRUE, std.lv = TRUE,
                    auto.fix.first = FALSE,
                    parameterization = "theta", group = "school", do.fit = FALSE)

## -----------------------------------------------------------------------------
pt <- parTable(fit_mg_nofit)
# Show loadings
pt[pt$op == "=~", c("lhs", "op", "rhs", "group", "free")]

## -----------------------------------------------------------------------------
# Show thresholds
pt[pt$op == "|", c("lhs", "op", "rhs", "group", "free")]

## -----------------------------------------------------------------------------
# Loadings: group 1 (Pasteur) and group 2 (Grant-White)
load_g1 <- pt$free[pt$op == "=~" & pt$group == 1 & pt$free > 0]
load_g2 <- pt$free[pt$op == "=~" & pt$group == 2 & pt$free > 0]

# Thresholds: group 1 and group 2
thresh_g1 <- pt$free[pt$op == "|" & pt$group == 1 & pt$free > 0]
thresh_g2 <- pt$free[pt$op == "|" & pt$group == 2 & pt$free > 0]

print(list(
    loadings_g1 = load_g1,
    loadings_g2 = load_g2,
    thresholds_g1 = thresh_g1,
    thresholds_g2 = thresh_g2
))

## -----------------------------------------------------------------------------
fit_pen_mg <- penalized_est(
    fit_mg_nofit,
    w = 0.03,
    pen_diff_id = list(
        loadings = rbind(load_g1, load_g2),
        thresholds = rbind(thresh_g1, thresh_g2)
    )
)
summary(fit_pen_mg)

## -----------------------------------------------------------------------------
# Loadings
load_ests_g1 <- as.numeric(coef(fit_pen_mg)[load_g1])
load_ests_g2 <- as.numeric(coef(fit_pen_mg)[load_g2])
load_mat <- rbind(load_ests_g1, load_ests_g2)
colnames(load_mat) <- names(coef(fit_pen_mg))[load_g1]
eff_load_diff <- composite_pair_loss(load_mat, fun = l0a)

# Thresholds
thresh_ests_g1 <- as.numeric(coef(fit_pen_mg)[thresh_g1])
thresh_ests_g2 <- as.numeric(coef(fit_pen_mg)[thresh_g2])
thresh_mat <- rbind(thresh_ests_g1, thresh_ests_g2)
colnames(thresh_mat) <- names(coef(fit_pen_mg))[thresh_g1]
eff_thresh_diff <- composite_pair_loss(thresh_mat, fun = l0a)

cat("Penalized Loading Estimates:\n")
print(load_mat, digits = 3)

cat("\nPenalized Threshold Estimates:\n")
print(thresh_mat, digits = 3)

cat("Effective number of non-invariant loadings:", eff_load_diff, "\n")
cat("Effective number of non-invariant thresholds:", eff_thresh_diff, "\n")