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

## ----population---------------------------------------------------------------
K <- 3
ipf <- 6
J <- K * ipf
cluster <- rep(seq_len(K), each = ipf)

set.seed(7)
mla <- matrix(0, J, K + 1)
mla[, 1] <- round(runif(J, .35, .75), 2)
for (k in seq_len(K)) {
  mla[cluster == k, k + 1] <- round(runif(ipf, .35, .70), 2)
}
while (max(rowSums(mla^2)) > .9) mla <- round(mla * .97, 3)

sim <- sim_fa(N = 800, mla = mla, phi = 0, rseed = 1)
Y <- sim$dat

## ----backbones----------------------------------------------------------------
Q0_AO <- matrix(-1L, J, 2)
for (k in 1:2) Q0_AO[which(cluster == k)[1:2], k] <- 1L

Q0_AZ <- Q0_AO
for (k in 1:2) {
  anchors <- which(cluster == k)[1:2]
  Q0_AZ[anchors, ] <- 0L
  Q0_AZ[anchors, k] <- 1L
}

Q0_show <- cbind(Q0_AO[1:8, ], Q0_AZ[1:8, ])
colnames(Q0_show) <- c("AO_1", "AO_2", "AZ_1", "AZ_2")
Q0_show

## ----step1--------------------------------------------------------------------
AO <- pefa(Q0_AO, Y, Kmin = 2, Kmax = 5, verbose = FALSE)
AZ <- pefa(Q0_AZ, Y, Kmin = 2, Kmax = 5, verbose = FALSE)

AO$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]
AZ$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]

## ----step1-transitions--------------------------------------------------------
AO$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
AZ$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]

## ----step1-persistence--------------------------------------------------------
round(AO$persistence$phi, 3)
round(AZ$persistence$phi, 3)
ssl(AO)[c("3", "4")]
ssl(AZ)[c("3", "4")]

## ----step1-reading------------------------------------------------------------
## Analysis-owned profiles. ELBO and BIC each use the full fitted window,
## a 20% gain cutoff, sustain one, a post-peak strict crossing, and no fallback.
profile_specs <- data.frame(
  profile = c(".85/r1", ".80/r2", ".70/r3"),
  phi_cut = c(.85, .80, .70), r = 1:3,
  practical_default = c(FALSE, TRUE, FALSE)
)
profile_sources <- 2:3

read_sweep <- function(x, phi_cut, r, sources) {
  K <- as.integer(x$sweep$K)
  converged <- x$sweep$converged

  read_count <- function(score) {
    usable <- length(K) >= 2L && all(is.finite(K)) && all(diff(K) == 1L) &&
      length(score) == length(K) && all(is.finite(score)) &&
      is.logical(converged) && all(converged %in% TRUE)
    if (!usable) return(list(usable = FALSE, Khat = NA_integer_))
    gain <- diff(score)
    gain_max <- max(gain)
    if (!(gain_max > 0)) return(list(usable = TRUE, Khat = NA_integer_))
    peak <- max(which(gain == gain_max))
    hit <- which(seq_along(gain) > peak & gain < .20 * gain_max)
    list(usable = TRUE,
         Khat = if (length(hit)) K[hit[1L]] else NA_integer_)
  }

  ELBO <- read_count(x$sweep$ELBO)
  BIC <- read_count(-x$sweep$BIC)
  C20 <- sort(unique(c(ELBO$Khat, BIC$Khat)[
    is.finite(c(ELBO$Khat, BIC$Khat))]))

  phi <- x$persistence$phi
  collision <- x$persistence$collision
  conv <- setNames(converged, K)
  state <- vapply(sources, function(k) {
    targets <- k + seq_len(r)
    edge <- vapply(targets, function(target) {
      from <- as.character(k); to <- as.character(target)
      if (!from %in% rownames(phi) || !to %in% colnames(phi) ||
          !isTRUE(conv[[from]]) || !isTRUE(conv[[to]]) ||
          is.na(collision[from, to])) return(NA_integer_)
      if (isTRUE(collision[from, to])) return(0L)
      value <- phi[from, to]
      if (!is.finite(value)) NA_integer_ else as.integer(value >= phi_cut)
    }, integer(1))
    if (any(edge %in% 0L)) 0L else
      if (length(edge) && all(edge %in% 1L)) 1L else NA_integer_
  }, integer(1))

  persistent <- which(state == 1L)
  if (!length(state)) {
    resolution <- "unresolved"
    Kp <- NA_integer_
  } else if (!length(persistent)) {
    resolution <- if (all(state %in% 0L)) "nonpersistent" else "unresolved"
    Kp <- NA_integer_
  } else {
    top <- max(persistent)
    higher <- if (top < length(state)) state[(top + 1L):length(state)] else 0L
    resolution <- if (all(higher %in% 0L)) "persistent" else "unresolved"
    Kp <- if (resolution == "persistent") sources[top] else NA_integer_
  }

  hit_ELBO <- is.finite(Kp) && is.finite(ELBO$Khat) && ELBO$Khat == Kp
  hit_BIC <- is.finite(Kp) && is.finite(BIC$Khat) && BIC$Khat == Kp
  support <- if (hit_ELBO && hit_BIC) "both" else if (hit_ELBO) "ELBO only" else
    if (hit_BIC) "BIC only" else "none"
  any_count_usable <- ELBO$usable || BIC$usable
  layer <- if (resolution == "persistent" && Kp %in% C20) "L1" else
    if (resolution == "persistent" && any_count_usable) "L2" else
    if (resolution == "nonpersistent") "L3" else
      "unclassified"
  gap <- if (resolution == "persistent" && length(C20))
    min(abs(C20 - Kp)) else NA_real_

  list(Khat_ELBO = ELBO$Khat, Khat_BIC = BIC$Khat, C20 = C20,
       count_usable = c(ELBO = ELBO$usable, BIC = BIC$usable),
       support = support, Kp = Kp, source_state = setNames(state, sources),
       resolution = resolution, layer = layer, g_CP = gap)
}

profile_row <- function(x, backbone, spec) {
  z <- read_sweep(x, spec$phi_cut, spec$r, profile_sources)
  data.frame(backbone = backbone, profile = spec$profile,
             practical_default = spec$practical_default,
             Khat_ELBO = z$Khat_ELBO, Khat_BIC = z$Khat_BIC,
             C20 = if (length(z$C20)) paste0("{", paste(z$C20,
               collapse = ","), "}") else "{}",
             support = z$support, Kp = z$Kp, layer = z$layer,
             g_CP = z$g_CP, row.names = NULL)
}

sweeps <- list(AO = AO, AZ = AZ)
profiles <- do.call(rbind, lapply(names(sweeps), function(backbone)
  do.call(rbind, lapply(seq_len(nrow(profile_specs)), function(i)
    profile_row(sweeps[[backbone]], backbone, profile_specs[i, ])))))
profiles

## ----step2-design-------------------------------------------------------------
default_AO <- profiles[profiles$backbone == "AO" &
                       profiles$profile == ".80/r2", , drop = FALSE]
if (nrow(default_AO) != 1L || !default_AO$layer %in% c("L1", "L2"))
  stop("The declared profile did not deliver a Step-2 proposal.")
K_step2 <- as.integer(default_AO$Kp)
K0      <- ncol(AO$Q0)
Lam     <- AO$loadings[[as.character(K_step2)]]
d     <- unique(colSums(AO$Q0 == 1L))
if (length(d) != 1L || d[1] < 1L) {
  stop("Choose one marker depth before building Q2.")
}

Q2 <- cbind(AO$Q0, matrix(-1L, J, K_step2 - K0))
colnames(Q2) <- paste0("F", seq_len(K_step2))
if (K_step2 > K0) {
  for (k in (K0 + 1L):K_step2) {
    top <- head(order(-abs(Lam[, k]), seq_len(J)), d)
    Q2[top, k] <- 1L
  }
}

c(K_step2 = K_step2, K0 = K0, depth = d)
Q2[rowSums(Q2 == 1L) > 0L, , drop = FALSE]

## ----step2-fit----------------------------------------------------------------
f_obl <- vbfa(Y, Q2)
f_bif <- vbfa(Y, Q2, bifactor = TRUE)

idx <- c("RMSEA", "SRMR", "CFI", "TLI", "AIC", "BIC")
round(rbind(
  oblique  = fit_stats(f_obl)[idx],
  bifactor = fit_stats(f_bif)[idx]
), 3)

## ----step2-effects------------------------------------------------------------
special_effects(f_bif)

## ----bifactor-sweep-----------------------------------------------------------
BI <- pefa(Q0_AO, Y, Kmin = 2, Kmax = 5, bifactor = TRUE, verbose = FALSE)
BI$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]
BI$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
round(BI$persistence$phi, 3)
round(BI$persistence$rmsd, 3)
ssl(BI)