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

## ----fit----------------------------------------------------------------------
sim <- sim_fa(N = 500, K = 3, ipf = 6, lam = .7, lac = .3, rseed = 1)
Y <- sim$dat

## an AZ (anchor-zero) design: two anchors per factor, each specified (1) on
## its own factor and fixed to zero on the other two; every non-anchor row is
## left entirely to the data
Q <- matrix(-1L, ncol(Y), 3)
for (k in 1:3) {
  a <- which(rep(1:3, each = 6) == k)[1:2]
  Q[a, ] <- 0L
  Q[a, k] <- 1L
}

fit <- vbfa(Y, Q)
fit

## ----access-------------------------------------------------------------------
round(fit$Lam[1:6, ], 2)   # posterior mean loadings
round(fit$pi[1:6, ], 2)    # PIPs of the unspecified entries
round(fit$Phi, 2)          # factor correlations (oblique by default)

## ----recovery-----------------------------------------------------------------
active <- (Q == 1) | (Q == -1 & fit$pi >= .5)
table(truth = sim$MLA != 0, active = active)

## ----vbfit--------------------------------------------------------------------
round(fit_stats(fit), 3)

## ----rank-adjustment, eval=FALSE----------------------------------------------
# fit_stats(fit, rank_adjust = TRUE, rank_max_J = 100)

## ----ld-----------------------------------------------------------------------
simLD <- sim_fa(N = 500, K = 3, ipf = 6, lam = .7, lac = .3, ecr = .3,
                rseed = 2)
fLD <- vbfa(simLD$dat, Q, ld = TRUE, max_it = 300, tolVal = 1e-3)

## the largest recovered residual edges, vs the planted pairs
Poff <- abs(fLD$Psi); Poff[lower.tri(Poff, diag = TRUE)] <- 0
which(Poff >= sort(Poff, decreasing = TRUE)[3], arr.ind = TRUE)
simLD$ofd_ind

## ----ld-restricted------------------------------------------------------------
J <- ncol(simLD$dat)
Qe <- matrix(-1L, J, J)
Qe[1:3, 1:3] <- 1L    # a known testlet: freely estimated among these 3 items
Qe[4:6, 7:9] <- 0L    # a block known to be residually independent
Qe[7:9, 4:6] <- 0L
isSymmetric(unname(Qe))

fRestricted <- vbfa(simLD$dat, Q, ld = TRUE, Qe = Qe, max_it = 300,
                    tolVal = 1e-3)

round(fRestricted$Psi[1:3, 1:3], 3)   # freely estimated: off-diagonal is not forced
round(fRestricted$Psi[4:6, 7:9], 3)   # fixed absent: driven to (numerical) zero

## ----ld-compare---------------------------------------------------------------
## fully exploratory loadings isolate the comparison to the residual side
Qexp <- matrix(-1L, J, 3)

fDiag <- vbfa(simLD$dat, Qexp)
fLDc  <- vbfa(simLD$dat, Qexp, ld = TRUE, max_it = 300, tolVal = 1e-3)

round(rbind(diagonal = fit_stats(fDiag)[c("BIC", "RMSEA")],
            ld       = fit_stats(fLDc)[c("BIC", "RMSEA")]), 3)

## ----ld-control---------------------------------------------------------------
fCtrl <- vbfa(simLD$dat, Q, ld = TRUE, max_it = 300, tolVal = 1e-3,
              ld_control = list(xi0 = c(0.1, 0.5, 1)))
fCtrl$converged   # TRUE if the final v0 stage met the tolerance

## ----mar-sim------------------------------------------------------------------
simM <- sim_fa(N = 400, K = 3, ipf = 6, lam = .7, lac = .3, rseed = 1)
Ym0  <- simM$dat

## items 5 and 11 go missing depending on the OBSERVED value of anchor items
## 1 and 7 (higher values make missingness more likely); the anchors
## themselves stay fully observed. This is MAR, not MCAR: the probability of
## missingness varies systematically with observed, not missing, data.
set.seed(42)
Ymar <- Ym0
p5  <- ifelse(Ym0[, 1] > stats::median(Ym0[, 1]), .40, .05)
p11 <- ifelse(Ym0[, 7] > stats::median(Ym0[, 7]), .40, .05)
Ymar[stats::rbinom(nrow(Ym0), 1, p5)  == 1, 5]  <- NA
Ymar[stats::rbinom(nrow(Ym0), 1, p11) == 1, 11] <- NA
sum(is.na(Ymar))

fMar <- vbfa(Ymar, Q)
fMar$preprocess$n_missing
round(fMar$Lam[c(1, 5, 7, 11), ], 2)   # loadings recovered despite the missingness

## ----mar-sim-clean------------------------------------------------------------
fClean <- vbfa(Ym0, Q)
max(abs(fMar$Lam - fClean$Lam))    # largest discrepancy
mean(abs(fMar$Lam - fClean$Lam))   # typical discrepancy is much smaller

## ----nlsy---------------------------------------------------------------------
data(nlsy27)
Yn <- as.matrix(nlsy27$dat)
dim(Yn)
sum(is.na(Yn))            # incomplete cells, handled in-loop

fn <- vbfa(Yn, nlsy27$Q)
fn
fn$preprocess$n_missing   # recorded on the fit
round(fn$Lam, 2)

## ----nlsy-cc------------------------------------------------------------------
nrow(Yn) - sum(stats::complete.cases(Yn))   # respondents listwise deletion drops

## ----softq--------------------------------------------------------------------
## a toy graded membership matrix for the simulated items
set.seed(3)
S <- matrix(runif(18 * 3, 0, .2), 18, 3)          # baseline noise
S[cbind(1:18, rep(1:3, each = 6))] <- runif(18, .55, .95)  # true memberships

Qsoft <- matrix(0L, 18, 3)          # S < .20 stays 0
Qsoft[S >= .20 & S < .80] <- -1L
Qsoft[S >= .80] <- 1L
table(Qsoft)

fsoft <- vbfa(Y, Qsoft)
fsoft$converged