or Floor Data

Documented in ggm_cfe plot_censoring

#' @import mvtnorm
#' @import qgraph
#' @importFrom Matrix nearPD
#' @import psych
#' @importFrom graphics par barplot
#' @importFrom utils globalVariables
#' @importFrom stats cor cov sd qnorm pnorm dnorm pt p.adjust uniroot
NULL

utils::globalVariables(c("variable", "pct", "type"))
# 
# ggmcfe.R  User-facing function for GGM with Ceiling/Floor Effects
# Liu (2026). Psychological Networks with Ceiling/Floor Effects.
# 
# Install required packages if needed:
#   install.packages(c("mvtnorm","qgraph","Matrix","psych"))
# Usage:
#   source("ggmcfe.R")
#   fit <- ggm_cfe(data, floor=0, ceiling=6)
#   print(fit); plot(fit); summary(fit)
# 

#  Internal helpers (not exported)

.correct_marginals <- function(y, a, b) {
  n <- length(y); l <- sum(y<=a); r <- sum(y>=b)
  ym <- y[y>a & y<b]; m <- length(ym)
  if (m < 5) return(list(mu=mean(y), sigma=sd(y), n_floor=l, n_ceil=r,
                         floor_pct=l/n, ceil_pct=r/n))
  Mp <- mean(ym); sp <- sd(ym); eps <- 1e-8
  ah <- qnorm(l/n+eps); bh <- qnorm(1-r/n-eps)
  fa <- dnorm(ah); fb <- dnorm(bh); dP <- pnorm(bh)-pnorm(ah)
  denom <- max(1+(ah*fa-bh*fb)/dP-((fa-fb)/dP)^2, 0.05)
  sig <- sqrt(sp^2/denom)
  list(mu=Mp+sig*(fb-fa)/dP, sigma=sig, n_floor=l, n_ceil=r,
       alpha=ah, beta=bh, floor_pct=l/n, ceil_pct=r/n)
}

# Gauss-Hermite quadrature nodes and weights (50-point, computed once at load)
.gh50 <- local({
  gh <- function(n = 50) {
    i <- seq_len(n - 1); J <- matrix(0, n, n)
    off <- sqrt(i / 2)
    J[cbind(i, i+1)] <- off; J[cbind(i+1, i)] <- off
    ev <- eigen(J, symmetric=TRUE)
    ord <- order(ev$values)
    list(nodes=ev$values[ord], weights=sqrt(pi)*ev$vectors[1,ord]^2)
  }
  gh(50)
})

# Expected mean of a clipped standard normal: E[min(max(Z,a),b)]
.clip_mean_std <- function(a, b) {
  fp <- if (is.finite(a)) a * pnorm(a)               else 0
  cp <- if (is.finite(b)) b * pnorm(b, lower.tail=FALSE) else 0
  fp + (dnorm(a) - dnorm(b)) + cp
}

# Expected mean of clip(N(mu,sd^2), a, b) for conditional mean computation
.clip_norm_mean <- function(mu, sd, a, b) {
  if (!is.finite(a) && !is.finite(b)) return(mu)
  za <- (a - mu) / sd; zb <- (b - mu) / sd
  fp <- if (is.finite(a)) a * pnorm(za)               else 0
  cp <- if (is.finite(b)) b * pnorm(zb, lower.tail=FALSE) else 0
  fp + mu*(pnorm(zb)-pnorm(za)) + sd*(dnorm(za)-dnorm(zb)) + cp
}

# Expected censored covariance by moment matching (Gauss-Hermite quadrature)
# Computes E[Z1*(rho) * Z2*(rho)] - E[Z1*(rho)] * E[Z2*(rho)]
# where Zi* = clip(Zi, ai, bi) and (Z1,Z2) ~ BVN(0,0,1,1,rho)
.cens_cov_mom <- function(rho, a1, b1, a2, b2, gh=.gh50) {
  rho <- max(min(rho, 0.9999), -0.9999)
  x   <- sqrt(2) * gh$nodes
  w   <- gh$weights / sqrt(pi)
  cx  <- pmin(pmax(x, a1), b1)
  sd_cond <- sqrt(max(1 - rho^2, 1e-10))
  ez_clip <- .clip_norm_mean(rho * x, sd_cond, a2, b2)
  cross <- sum(w * cx * ez_clip)
  cross - .clip_mean_std(a1, b1) * .clip_mean_std(a2, b2)
}

# Pairwise MOM estimator: find rho such that model-implied censored cov = observed
.est_rho_mom <- function(yj, yk, mj, sj, mk, sk, aj, bj, ak, bk) {
  # Standardize thresholds
  aj_s <- (aj - mj) / sj; bj_s <- (bj - mj) / sj
  ak_s <- (ak - mk) / sk; bk_s <- (bk - mk) / sk
  # Observed censored covariance, scaled by corrected SDs
  target <- cov(yj, yk) / (sj * sk)
  f <- function(rho) .cens_cov_mom(rho, aj_s, bj_s, ak_s, bk_s) - target
  lo <- -0.999; hi <- 0.999
  flo <- f(lo); fhi <- f(hi)
  if (!is.finite(flo) || !is.finite(fhi)) return(cor(yj, yk))
  if (flo > 0) return(lo)
  if (fhi < 0) return(hi)
  uniroot(f, interval=c(lo, hi), tol=1e-7)$root
}

.pcor_from_omega <- function(Om) {
  p<-nrow(Om); PC<-matrix(0,p,p)
  for(j in 1:(p-1)) for(k in (j+1):p){
    v<- -Om[j,k]/sqrt(Om[j,j]*Om[k,k]); PC[j,k]<-PC[k,j]<-v}; PC
}

# 
#' Gaussian Graphical Model with Ceiling/Floor Effect Correction
#'
#' Estimates a psychological network (GGM) from data with ceiling and/or floor
#' effects using a two-step pairwise censored-normal correction. Step 1 corrects
#' item means and variances via truncated-normal moment matching (Liu & Wang, 2021).
#' Step 2 estimates pairwise latent correlations by solving the moment equation
#' E[cov*(rho)] = observed censored covariance, evaluated using 50-point
#' Gauss-Hermite quadrature. The corrected covariance matrix is supplied as input
#' to EBICglasso or significance-based edge selection, preserving the standard
#' qgraph/bootnet workflow.
#'
#' @param data  A data frame or numeric matrix (n x p). Rows = observations,
#'              columns = variables (items/nodes).
#' @param floor  Floor threshold(s). Either a single value applied to all
#'               variables, or a named numeric vector with one value per column.
#'               Use NULL to indicate no floor effect. Default: NULL.
#' @param ceiling  Ceiling threshold(s). Same format as \code{floor}.
#'                 Default: NULL.
#' @param method  Estimation method(s): "EBICglasso" (default), "FDR",
#'                "Bonferroni", or "all" to run all three.
#' @param gamma   EBIC hyperparameter (0 = BIC, 0.5 = default). Only used
#'                when method includes "EBICglasso".
#' @param nlambda Number of tuning parameters searched by EBICglasso.
#'                Default: 100.
#' @param lambda.min.ratio Smallest lambda searched, as a fraction of the
#'                largest lambda. Passed to \code{qgraph::EBICglasso}.
#'                Default: 0.01.
#' @param threshold Logical. Passed to \code{qgraph::EBICglasso}; \code{TRUE}
#'                enforces higher specificity at the cost of sensitivity.
#'                Default: FALSE.
#' @param alpha   Significance level for p-value methods. Default: 0.05.
#' @param verbose Logical. Print progress and diagnostics. Default: TRUE.
#'
#' @return An object of class \code{"ggm_cfe"} containing:
#'   \item{network}{Partial correlation matrix from the corrected method
#'                  (primary estimator: EBICglasso or first specified method).}
#'  \item{network_corrected}{Partial correlation matrix from naive EBICglasso
#'                        (ignoring ceiling/floor effects).}
#'   \item{network_naive}{Partial correlation matrix estimated from the raw
#'                  (uncorrected) sample covariance using the primary method.}
#'   \item{networks}{Named list of partial correlation matrices for all
#'                   requested methods.}
#'  \item{networks_corrected}{Partial correlation matrix from naive EBICglasso
#'                        (ignoring ceiling/floor effects).}
#'  \item{networks_naive}{Named list of naive partial correlation matrices
#'                  (one per requested method), estimated without
#'                  ceiling/floor correction.}
#'   \item{Sigma_corrected}{The corrected covariance matrix (p x p).}
#'   \item{Sigma_naive}{The raw sample covariance matrix (p x p).}
#'   \item{diagnostics}{List with per-variable censoring statistics, nearPD
#'                      correction magnitude, and floor/ceiling thresholds.
#'                      \describe{
#'       \item{censoring}{Data frame (p rows) with per-variable censoring
#'                  statistics: variable, \code{n_floor}, \code{n_ceil},
#'                  \code{pct_floor}, \code{pct_ceil}, naive and corrected
#'                  means (\code{mu_naive}, \code{mu_corrected}), and
#'                  standard deviations (\code{sd_naive},
#'                  \code{sd_corrected}).}
#'       \item{npd_correction}{Relative Frobenius-norm magnitude of the
#'                  nearest positive-definite projection applied to the
#'                  corrected covariance matrix. Values above 0.05 suggest
#'                  instability and should be interpreted cautiously.}
#'       \item{a_vec}{Length-p numeric vector of effective floor thresholds
#'                  used internally (with \code{NULL} inputs replaced by
#'                  data-driven lower bounds).}
#'       \item{b_vec}{Length-p numeric vector of effective ceiling thresholds
#'                  (with \code{NULL} inputs replaced by data-driven upper
#'                  bounds).}
#'     }
#' }
#'   \item{method}{Character vector of estimation method(s) run.}      
#'   \item{gamma}{EBIC hyperparameter value used.}   
#'   \item{nlambda}{Number of regularization parameters searched.}  
#'   \item{lambda.min.ratio}{Smallest lambda as a fraction of the largest.} 
#'   \item{threshold}{Logical; whether EBICglasso hard-thresholding was applied.}     
#'   \item{alpha}{Significance level for p-value methods.}    
#'   \item{n}{The number of observations.}     
#'   \item{p}{The number of variables.}     
#'   \item{varnames}{The character vector of variable (column) names.}   
#'   \item{data}{Original data matrix.}
#'   \item{call}{Matched call.}
#'
#' @examples
#' # Simulate data with ceiling effects
#' set.seed(42)
#' Y <- MASS::mvrnorm(200, rep(0,5), diag(5) + 0.3)
#' Y_cens <- pmin(Y, 1.0)   # ceiling at 1.0 SD
#' fit <- ggm_cfe(Y_cens, floor=NULL, ceiling=1.0)
#' print(fit)
#' plot(fit)
#'
#' # With psych::bfi personality data
#' # library(psych)
#' # fit <- ggm_cfe(bfi[,1:25], floor=1, ceiling=6)

#' @export
ggm_cfe <- function(data,
                    floor   = NULL,
                    ceiling = NULL,
                    method  = "EBICglasso",
                    gamma   = 0.5,
                    nlambda = 100,
                    lambda.min.ratio = 0.01,
                    threshold = FALSE,
                    alpha   = 0.05,
                    verbose = TRUE) {
  cl <- match.call()
  
  # Input validation
  if (is.data.frame(data)) data <- as.matrix(data)
  if (!is.numeric(data)) stop("'data' must be a numeric matrix or data frame.")
  
  cc <- complete.cases(data)
  if (!all(cc)) {
    warning(sprintf(
      "Removed %d row(s) with missing values before estimation.",
      sum(!cc)
    ), call. = FALSE)
    data <- data[cc, , drop = FALSE]
  }
  if (any(!is.finite(data))) {
    stop("'data' must not contain infinite values.")
  }
  
  n <- nrow(data); p <- ncol(data)
  if (n < p + 3) warning("n is small relative to p; results may be unstable.")
  if (any(apply(data, 2, sd) == 0)) {
    stop("'data' contains one or more constant columns.")
  }
  
  vnames <- if (!is.null(colnames(data))) colnames(data) else paste0("V", seq_len(p))
  colnames(data) <- vnames
  
  method <- match.arg(method, c("EBICglasso","FDR","Bonferroni","all"),
                      several.ok=FALSE)
  run_methods <- if (method=="all") c("EBICglasso","FDR","Bonferroni") else method
  
  # Expand floor/ceiling to length-p vectors
  expand_thresh <- function(th, name) {
    if (is.null(th)) return(rep(-Inf, p))   # effectively no threshold
    if (length(th)==1) return(rep(as.numeric(th), p))
    if (length(th)==p) return(as.numeric(th))
    stop(sprintf("'%s' must be NULL, a single value, or a vector of length p=%d.", name, p))
  }
  a_vec <- expand_thresh(floor,   "floor")
  b_vec <- expand_thresh(ceiling, "ceiling")
  # Replace +/- Inf with extreme data-driven values for computation
  a_vec[!is.finite(a_vec)] <- apply(data, 2, min)[!is.finite(a_vec)] - 1e6
  b_vec[!is.finite(b_vec)] <- apply(data, 2, max)[!is.finite(b_vec)] + 1e6
  
  #  Step 1: Marginal corrections 
  if (verbose) cat("Step 1: Correcting marginal distributions...\n")
  mg <- lapply(seq_len(p), function(j) .correct_marginals(data[,j], a_vec[j], b_vec[j]))
  mu_t  <- sapply(mg, `[[`, "mu")
  sig_t <- sapply(mg, `[[`, "sigma")
  
  cens_df <- data.frame(
    variable   = vnames,
    n_floor    = sapply(mg, `[[`, "n_floor"),
    n_ceil     = sapply(mg, `[[`, "n_ceil"),
    pct_floor  = round(sapply(mg, `[[`, "n_floor") / n * 100, 1),
    pct_ceil   = round(sapply(mg, `[[`, "n_ceil")  / n * 100, 1),
    mu_naive   = colMeans(data),
    mu_corrected = mu_t,
    sd_naive   = apply(data, 2, sd),
    sd_corrected = sig_t
  )
  if (verbose) {
    cat(sprintf("  Variables: %d  |  n: %d\n", p, n))
    cat(sprintf("  Mean floor%%: %.1f  |  Mean ceiling%%: %.1f\n",
                mean(cens_df$pct_floor), mean(cens_df$pct_ceil)))
  }
  
  # Step 2: Pairwise rho estimation (MOM via Gauss-Hermite quadrature) 
  if (verbose) cat("Step 2: Estimating pairwise correlations (MOM)...\n")
  rho_mat <- diag(p); total_pairs <- p*(p-1)/2; cnt <- 0
  for (j in seq_len(p-1)) for (k in (j+1):p) {
    cnt <- cnt+1
    r <- tryCatch(
      .est_rho_mom(data[,j],data[,k],mu_t[j],sig_t[j],mu_t[k],sig_t[k],
                   a_vec[j],b_vec[j],a_vec[k],b_vec[k]),
      error=function(e) cor(data[,j],data[,k]))
    rho_mat[j,k] <- rho_mat[k,j] <- r
    if (verbose && cnt %% max(1,floor(total_pairs/5)) == 0)
      cat(sprintf("  Pair %d / %d\r", cnt, total_pairs))
  }
  if (verbose) cat(sprintf("  Pair %d / %d\n", total_pairs, total_pairs))
  
  # Step 3: Corrected covariance + nearPD 
  Sig_raw <- outer(sig_t, sig_t) * rho_mat
  rownames(Sig_raw) <- colnames(Sig_raw) <- vnames
  pd_obj  <- tryCatch(Matrix::nearPD(Sig_raw, corr=FALSE, keepDiag=TRUE),
                      error=function(e) list(mat=Sig_raw+diag(1e-4,p), converged=FALSE))
  Sig_pd  <- as.matrix(pd_obj$mat)
  npd_rel <- sqrt(sum((Sig_pd-Sig_raw)^2)) / max(sqrt(sum(Sig_raw^2)), 1e-10)
  if (verbose) {
    cat(sprintf("Step 3: nearPD correction: %.2f%%", 100*npd_rel))
    if (is.finite(npd_rel) && npd_rel > 0.05) cat("  *** WARNING: >5%%  interpret with caution")
    cat("\n")
  }
  
  # Step 4: Network estimation 
  if (verbose) cat("Step 4: Estimating network...\n")
  S_naive <- cov(data)
  .run_gl  <- function(S, n_) tryCatch(
    EBICglasso(S, n=n_, gamma=gamma, nlambda=nlambda,
               lambda.min.ratio=lambda.min.ratio, threshold=threshold,
               returnAllResults=FALSE),
    error=function(e){ warning("EBICglasso failed; returning empty network."); matrix(0,p,p) })
  .pval_net <- function(S, n_eff, method_) {
    Om <- tryCatch(solve(S), error=function(e) solve(S+diag(1e-4,p)))
    D  <- sqrt(diag(Om)); PC <- -Om/outer(D,D); diag(PC) <- 0
    df <- max(floor(n_eff)-p-2, 1); ut <- upper.tri(PC)
    rv <- PC[ut]; tv <- rv*sqrt(df)/sqrt(pmax(1-rv^2,1e-8))
    pv <- 2*pt(-abs(tv), df=df); nt <- sum(ut)
    sig <- if(method_=="Bonferroni") pv < alpha/nt else p.adjust(pv,"BH") < alpha
    M <- matrix(0,p,p); M[ut] <- PC[ut]*sig; M+t(M)
  }
  
  nets_prop <- list(); nets_naive <- list()
  for (meth in run_methods) {
    nets_prop[[meth]] <- switch(meth,
                                EBICglasso = .run_gl(Sig_pd, n),
                                FDR        = .pval_net(Sig_pd, n, "FDR"),
                                Bonferroni = .pval_net(Sig_pd, n, "Bonferroni")
    )
    nets_naive[[meth]] <- switch(meth,
                                 EBICglasso = .run_gl(S_naive, n),
                                 FDR        = .pval_net(S_naive, n, "FDR"),
                                 Bonferroni = .pval_net(S_naive, n, "Bonferroni")
    )
    rownames(nets_prop[[meth]]) <- colnames(nets_prop[[meth]]) <- vnames
    rownames(nets_naive[[meth]])<- colnames(nets_naive[[meth]])<- vnames
  }
  
  structure(list(
    network          = nets_prop[[run_methods[1]]],
    network_corrected = nets_prop[[run_methods[1]]],
    network_naive    = nets_naive[[run_methods[1]]],
    networks         = nets_prop,
    networks_corrected = nets_prop,
    networks_naive   = nets_naive,
    Sigma_corrected  = Sig_pd,
    Sigma_naive      = S_naive,
    diagnostics      = list(
      censoring      = cens_df,
      npd_correction = npd_rel,
      a_vec          = a_vec,
      b_vec          = b_vec
    ),
    method    = run_methods,
    gamma     = gamma,
    nlambda   = nlambda,
    lambda.min.ratio = lambda.min.ratio,
    threshold = threshold,
    alpha     = alpha,
    n         = n, p = p,
    varnames  = vnames,
    data      = data,
    call      = cl
  ), class = "ggm_cfe")
}

# S3 methods 
#' @exportS3Method print ggm_cfe
print.ggm_cfe <- function(x, ...) {
  cat("\n=== GGM with Ceiling/Floor Correction (ggm_cfe) ===\n")
  cat(sprintf("Variables: %d  |  Observations: %d\n", x$p, x$n))
  cat(sprintf("Method(s): %s\n", paste(x$method, collapse=", ")))
  d <- x$diagnostics$censoring
  cat(sprintf("\nCensoring summary (floor / ceiling %%):\n"))
  cat(sprintf("  Mean:  %.1f%% / %.1f%%\n", mean(d$pct_floor), mean(d$pct_ceil)))
  cat(sprintf("  Range: %.1f-%.1f%% / %.1f-%.1f%%\n",
              min(d$pct_floor), max(d$pct_floor), min(d$pct_ceil), max(d$pct_ceil)))
  cat(sprintf("\nnearPD correction: %.2f%%", 100*x$diagnostics$npd_correction))
  if (is.finite(x$diagnostics$npd_correction) && x$diagnostics$npd_correction > 0.05)
    cat("  *** >5%% threshold check stability")
  cat("\n\n")
  
  net <- x$network
  n_edges <- sum(net[upper.tri(net)] != 0)
  cat(sprintf("Estimated edges (corrected): %d  |  Density: %.1f%%\n",
              n_edges, 100*n_edges/(x$p*(x$p-1)/2)))
  n_edges_naive <- sum(x$network_naive[upper.tri(x$network_naive)] != 0)
  cat(sprintf("Estimated edges (naive):    %d  |  Density: %.1f%%\n",
              n_edges_naive, 100*n_edges_naive/(x$p*(x$p-1)/2)))
  invisible(x)
}

#' @exportS3Method summary ggm_cfe
summary.ggm_cfe <- function(object, ...) {
  cat("\n=== Censoring Statistics per Variable ===\n")
  print(object$diagnostics$censoring[, c("variable","pct_floor","pct_ceil",
                                         "mu_naive","mu_corrected",
                                         "sd_naive","sd_corrected")],
        row.names=FALSE, digits=3)
  cat(sprintf("\nnearPD correction: %.3f%%\n", 100*object$diagnostics$npd_correction))
  invisible(object)
}

#' @exportS3Method plot ggm_cfe
plot.ggm_cfe <- function(x, which="both",
                         layout="spring", title_prefix="", ...) {
  # which: "corrected", "naive", or "both" (side-by-side).
  # "proposed" is accepted for backward compatibility.
  if (which == "proposed") which <- "corrected"
  if (which == "both") {
    par_old <- par(mfrow=c(1,2))
    on.exit(par(par_old))
  }
  if (which %in% c("corrected","both"))
    qgraph(x$network,
           layout     = layout,
           title      = paste0(title_prefix, "Corrected"),
           posCol     = "#2980B9", negCol = "#E74C3C",
           maximum    = 0.5, labels = x$varnames, ...)
  if (which %in% c("naive","both"))
    qgraph(x$network_naive,
           layout     = layout,
           title      = paste0(title_prefix, "Naive (uncorrected)"),
           posCol     = "#2980B9", negCol = "#E74C3C",
           maximum    = 0.5, labels = x$varnames, ...)
  invisible(x)
}

# Convenience: censoring profile plot 
#' Censoring Profile Plot for ggm_cfe Objects
#'
#' Plots the floor and ceiling proportions for each variable in a
#' \code{ggm_cfe} object. Uses \pkg{ggplot2} if available, otherwise
#' falls back to base graphics.
#'
#' @param x An object of class \code{"ggm_cfe"} as returned by
#'   \code{\link{ggm_cfe}}.
#'
#' @return Invisibly returns \code{x}.
#'
#' @examples
#' set.seed(42)
#' Y <- MASS::mvrnorm(200, rep(0, 5), diag(5) + 0.3)
#' Y_cens <- pmin(Y, 1.0)
#' fit <- ggm_cfe(Y_cens, floor = NULL, ceiling = 1.0)
#' plot_censoring(fit)
#'
#' @export
plot_censoring <- function(x) {
  d <- x$diagnostics$censoring
  d_long <- rbind(
    data.frame(variable=d$variable, type="Floor", pct=d$pct_floor),
    data.frame(variable=d$variable, type="Ceiling", pct=d$pct_ceil)
  )
  d_long$variable <- factor(d_long$variable, levels=rev(x$varnames))
  if (!requireNamespace("ggplot2", quietly=TRUE)) {
    barplot(rbind(d$pct_floor, d$pct_ceil), beside=TRUE, names.arg=d$variable,
            col=c("#3498DB","#E74C3C"), legend.text=c("Floor","Ceiling"),
            main="Censoring proportions (%)", las=2, ylab="%")
    return(invisible(x))
  }
  p <- ggplot2::ggplot(d_long, ggplot2::aes(x=variable, y=pct, fill=type)) +
    ggplot2::geom_col(position="dodge", alpha=0.85) +
    ggplot2::coord_flip() +
    ggplot2::scale_fill_manual(values=c(Floor="#3498DB", Ceiling="#E74C3C")) +
    ggplot2::geom_hline(yintercept=c(10,20,30), linetype="dashed",
                        color="grey60", linewidth=0.4) +
    ggplot2::labs(title="Ceiling and Floor Proportions by Variable",
                  x=NULL, y="Proportion (%)", fill="Effect") +
    ggplot2::theme_bw(base_size=11) +
    ggplot2::theme(legend.position="bottom")
  print(p); invisible(x)
}