# Analysis Clinic Case 002: a controlled lme4 singular-fit example
#
# Purpose: reproduce lme4's boundary (singular) message, distinguish
# singularity from optimizer nonconvergence, and show how isSingular(),
# VarCorr(), and rePCA() diagnose an unsupported random-slope dimension.
# The data are simulated with random intercept variation but no random slope
# variation, so removing the random slope is justified by the known design of
# this teaching example. Real analyses require design- and theory-based review.

suppressPackageStartupMessages(library(lme4))

set.seed(1)

n_groups <- 30
n_per_group <- 8
group <- factor(rep(seq_len(n_groups), each = n_per_group))
x <- rep(as.numeric(scale(seq_len(n_per_group), scale = FALSE)), n_groups)
random_intercept <- rnorm(n_groups, mean = 0, sd = 1)
y <- 5 + 0.6 * x + random_intercept[group] + rnorm(n_groups * n_per_group)
dat <- data.frame(y, x, group)

# This fitted model asks the data to estimate both random-intercept and
# random-slope variation even though the generating process contains no random
# slope variation. lme4 reports the singularity as a message, not a warning.
fit_singular_messages <- character()
fit_singular <- withCallingHandlers(
  lmer(y ~ x + (1 + x | group), data = dat),
  message = function(m) {
    fit_singular_messages <<- c(fit_singular_messages, conditionMessage(m))
    invokeRestart("muffleMessage")
  }
)

# In this controlled example, a random-intercept model matches the known data-
# generating structure. This is not a rule to delete random slopes whenever a
# singular fit appears in an empirical analysis.
fit_supported <- lmer(y ~ x + (1 | group), data = dat)

singular_vc <- as.data.frame(VarCorr(fit_singular))
slope_variance <- singular_vc$vcov[
  singular_vc$grp == "group" & singular_vc$var1 == "x" & is.na(singular_vc$var2)
]
intercept_slope_correlation <- singular_vc$sdcor[
  singular_vc$grp == "group" &
    singular_vc$var1 == "(Intercept)" & !is.na(singular_vc$var2) &
    singular_vc$var2 == "x"
]
pca_sd <- rePCA(fit_singular)$group$sd

cat("R version:", R.version.string, "\n")
cat("lme4 version:", as.character(packageVersion("lme4")), "\n")
cat("seed: 1\n")
cat("groups:", n_groups, "\n")
cat("observations per group:", n_per_group, "\n")
cat("total observations:", nrow(dat), "\n\n")

cat("RANDOM-INTERCEPT + RANDOM-SLOPE FIT\n")
cat("message:", paste(trimws(fit_singular_messages), collapse = " | "), "\n")
cat("optimizer convergence code:", fit_singular@optinfo$conv$opt, "\n")
cat("lme4 post-fit messages:",
    if (is.null(fit_singular@optinfo$conv$lme4$messages)) "none" else
      paste(fit_singular@optinfo$conv$lme4$messages, collapse = " | "), "\n")
cat("isSingular:", isSingular(fit_singular), "\n")
cat("random-slope variance:", format(slope_variance, scientific = TRUE), "\n")
cat("intercept-slope correlation:", round(intercept_slope_correlation, 6), "\n")
cat("rePCA standard deviations:", paste(format(pca_sd, scientific = TRUE), collapse = ", "), "\n")
cat("fixed slope estimate:", round(unname(fixef(fit_singular)["x"]), 6), "\n\n")

cat("RANDOM-INTERCEPT FIT (matches the known generating structure)\n")
cat("isSingular:", isSingular(fit_supported), "\n")
cat("fixed slope estimate:", round(unname(fixef(fit_supported)["x"]), 6), "\n\n")

cat("INTERPRETATION\n")
cat("The optimizer completed successfully, but the random-effects covariance matrix is rank deficient.\n")
cat("The near-zero PCA dimension and boundary correlation show that the data do not distinguish every requested random-effect dimension.\n")
cat("The random-intercept model is justified here only because the simulation is known to contain no random-slope variation.\n")
cat("An empirical analysis must consider the design, grouping structure, focal estimand, and planned random effects before simplifying.\n")

stopifnot(
  any(grepl("boundary \\(singular\\) fit", fit_singular_messages)),
  fit_singular@optinfo$conv$opt == 0,
  any(grepl("boundary \\(singular\\) fit", fit_singular@optinfo$conv$lme4$messages)),
  isSingular(fit_singular),
  min(pca_sd) < 1e-4,
  slope_variance < 0.01,
  abs(intercept_slope_correlation) > 0.999,
  !isSingular(fit_supported)
)