R/umxSummary.MxModel.R

In umx: Structural Equation Modeling and Twin Modeling in R

#' Shows a compact, publication-style, summary of a RAM model
#'
#' Report the fit of a model in a compact form suitable for a journal. 
#' It reports parameters in a markdown or html table (optionally standardized), and fit indices
#' RMSEA (an absolute fit index, comparing the model to a perfect model) and CFI and TLI (incremental fit indices comparing a model with the worst fit).
#' 
#' `umxSummary` alerts you when model fit is worse than accepted criterion (TLI >= .95 and RMSEA <= .06; (Hu & Bentler, 1999; Yu, 2002).
#' 
#' Note: For some (multi-group) models, you will need to fall back on [summary()]
#'
#' **Robust inference (`uncertainty = "MLR"`)**
#'
#' This option means **engine-appropriate robust inference**, not one sandwich formula for every fit type:
#' \itemize{
#'   \item **ML / FIML (raw continuous):** casewise score sandwich SEs via [OpenMx::imxRobustSE()], plus robust CFI/TLI/RMSEA via [xmu_robust_ML_fit()] (Yuan-Bentler / Brosseau-Liard & Savalei scaling).
#'   \item **WLS / DWLS / ULS:** do **not** use the ML casewise sandwich. Parameter SEs are the usual **WLS/GMM moment sandwich** already returned by OpenMx (\eqn{W}, \eqn{\Delta}, \eqn{\Gamma}). Robust CFI/TLI/RMSEA are applied when \code{output$implied_jacobian} is present via [xmu_robust_WLS_fit()] (Satorra-Bentler 2010 for continuous WLS; Savalei 2021 catML for ordinal WLS). A statistical note states that SEs are WLS asymptotic (moment sandwich).
#' }
#'
#' **Robust Fit Statistics for ML**
#' Standard Maximum Likelihood (ML) is powerful and efficient, with chi2, CFI, TLI, and RMSEA measures of fit. Non-normality and heteroscedasticity distort the sampling distribution of the scores and the ML test statistic itself, necessitating changes to SEs and fit statistics. For ML, `uncertainty = "MLR"` implements the sandwich estimator for SEs and MLR-style fit scaling (Yuan-Bentler / Savalei 2018-style robust indices).
#'
#' **Robust Fit Statistics for WLS/DWLS (The GenomicMx Ecosystem)**
#' 
#' Simulation work shows CFI/TLI behave differently under DWLS/WLS than under ML
#' (Shi et al., 2020): incremental indices use an independence baseline under the
#' weight matrix, so conventional cutoffs (e.g. CFI > 0.95) break down for continuous WLS
#' and Genomic SEM. Prefer **SRMR** and residual inspection for absolute fit; use
#' [umxCompare()] for nested tests. See `?umxCompare` for the full guidance.
#'
#' *Implementation:* when a Jacobian is available, `umxSummary` routes WLS models
#' through [xmu_robust_WLS_fit()] (independence baseline + trace scaling). Parameter SEs
#' remain the OpenMx WLS moment sandwich (see `uncertainty = "MLR"` above). Robust
#' CFI/TLI/RMSEA are **reported** for continuous WLS but should be treated as descriptive
#' (no Hu-Bentler cutoffs). Genomic SEM has the same rule, with even less trust in RMSEA
#' under bookkeeping N.
#'
#' **Two statistic families (by design).** WLS output separates omnibus testing from fit indices:
#' \itemize{
#'   \item **Displayed \eqn{\chi^2} and \eqn{p}** - Satorra-Bentler (2010) scaled WLS test when robust fit succeeds (continuous and ordinal).
#'   \item **Displayed CFI, TLI, RMSEA** - robust/scaled indices; interpretation depends on continuous vs ordinal (below).
#' }
#'
#' **Continuous WLS** (all manifests numeric): SB (2010) scaled pipeline
#' (\eqn{\chi^2 = F/c} from `output$fit`, same Chi as [umxCompare()] when Jacobians
#' are present). Saturated residual df (ChiDoF = 0): Chi = 0, p = NA, CFI = TLI = 1,
#' RMSEA = 0 (not NA). Printed note: no conventional cutoffs; prefer SRMR; nested
#' Strict SB via [umxCompare()].
#' **Ordinal / categorical WLS** (at least one `ordered`/`factor` on raw data): robust CFI/TLI/RMSEA use **Savalei (2021)** catML (Categorical Maximum Likelihood) corrections; **Hu & Bentler (1999) cutoffs (CFI > .95, RMSEA < .06) apply to these robust indices**. Displayed \eqn{\chi^2}/\eqn{p} remain SB-scaled WLS omnibus tests.
#' \itemize{
#'   \item **catML (Categorical Maximum Likelihood):** Evaluates parameter estimates inside a Maximum Likelihood correlation scaffold to obtain non-inflated fit indices for ordinal data.
#'   \item **\eqn{c_{\text{model}}} and \eqn{c_{\text{null}}}:** Mean-and-variance scaling correction multipliers for model and baseline discrepancies (values near 1.0 indicate sample sampling variability matches standard ML expectations; values > 1.0 adjust for categorical overdispersion).
#'   \item **Suggested Manuscript Copy:** *"Model fit was evaluated using DWLS with robust categorical ML fit indices (Savalei, 2021; CFI = X.XX, RMSEA = X.XX)."*
#' }
#'
#' Requires `implied_jacobian` after `mxRun` on WLS. If robust computation fails, falls back to OpenMx `summary()` with a note.
#' 
#' **CIs and Identification**
#' This function uses the standard errors reported by OpenMx to produce the CIs you see in umxSummary.
#' These are used to derive confidence intervals based on the formula 95%CI = estimate +/- 1.96*SE.
#' 
#' Sometimes SEs appear NA. This may reflect a model which is not identified (see <http://davidakenny.net/cm/identify.htm>).
#' This can include empirical under-identification - for instance two factors
#' that are essentially identical in structure. Use [OpenMx::mxCheckIdentification()] to check identification.
#' 
#' Solutions: If there are paths estimated at or close to zero, this suggests that fixing one or two of 
#' these to zero may fix the standard error calculation.
#' 
#' If factor loadings can flip sign and provide identical fit, this creates another form of 
#' under-identification and can break confidence interval estimation.
#' *Solution*: Fixing a factor loading to 1 and estimating factor variances can help here.
#'
#' @aliases umxSummary.MxModel umxSummary.MxRAMModel
#' @param model The [OpenMx::mxModel()] whose fit will be reported
#' @param std If TRUE, model is standardized (Default FALSE, NULL means "don't show").
#' @param digits How many decimal places to report (Default 2)
#' @param report If "html", then show results in browser (default = "markdown")
#' @param uncertainty What type of uncertainty/inference packaging to report: `"SE"` (default SEs from the fit), `"MLR"` (engine-appropriate robust inference: ML casewise sandwich SEs + robust AFIs; for WLS, moment-sandwich SEs + robust WLS AFIs when available), `"CI"` (profile likelihood CIs), or `"none"`.
#' @param means Whether to include means in the summary (TRUE)
#' @param residuals Whether to include residuals in the summary (TRUE)
#' @param filter whether to show significant paths (SIG) or NS paths (NS) or all paths (ALL)
#' @param RMSEA_CI Whether to compute the CI on RMSEA (Defaults to FALSE)
#' @param refModels Saturated models if needed for fit indices (see example below:
#'  If NULL will be computed on demand. If FALSE will not be computed.
#' @param matrixAddresses Whether to show "matrix address" columns (Default = FALSE)
#' @param SE Deprecated. Please use `uncertainty` instead.
#' @param ... Other parameters to control model summary
#' @family Summary functions
#' @seealso - [umxRAM()], [xmu_robust_WLS_fit()]
#' @references
#' * Hu, L., & Bentler, P. M. (1999). Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. *Structural Equation Modeling*, **6**, 1-55. 
#' * Satorra, A., & Bentler, P. M. (2010). Ensuring positiveness of the scaled difference chi-square test statistic. *Psychometrika*, **75**(2), 243-248.
#' * Satorra, A. and Bentler P. M. (1994). Corrections to test statistics and standard errors in covariance structure analysis. Latent variables analysis: Applications for developmental research. A. von Eye and C. C. Clogg, Sage: 399-419.
#' * Savalei, V. (2018). On the computation of the RMSEA and CFI from the mean-and-variance corrected test statistic with nonnormal data in SEM. *Multivariate Behavioral Research*, **53**(3), 419--429.
#' * Savalei, V. (2021). Improving fit indices in SEM with categorical data. *Multivariate Behavioral Research*, **56**(3), 390--407.
#' * Brosseau-Liard, P. E., & Savalei, V. (2012). Adjusting incremental fit indices for nonnormality. *Multivariate Behavioral Research*, **47**(5), 647--677.
#' * Yu, C.Y. (2002). Evaluating cutoff criteria of model fit indices for latent variable models with binary and continuous outcomes. University of California, Los Angeles. Retrieved from <https://www.statmodel.com/download/Yudissertation.pdf>
#' * Yuan, K.-H. and P. M. Bentler (2000). "5. Three Likelihood-Based Methods for Mean and Covariance Structure Analysis with Nonnormal Missing Data." Sociological Methodology 30(1): 165-200.
#' @export
#' @import OpenMx
#' @return - parameterTable returned invisibly, if estimates requested
#' @examples
#' \dontrun{
#' require(umx)
#' data(demoOneFactor)
#' manifests = names(demoOneFactor)
#' m1 = umxRAM("One Factor", data = demoOneFactor, type = "cov",
#'  umxPath("G", to = manifests),
#'  umxPath(var = manifests),
#'  umxPath(var = "G", fixedAt = 1)
#' )
#' umxSummary(m1, std = TRUE)
#' # output as latex
#' umx_set_table_format("latex")
#' umxSummary(m1, std = TRUE)
#' umx_set_table_format("markdown")
#' # output as raw
#' umxSummary(m1, std = FALSE)
#' 
#' # switch to a raw data model
#' m1 = umxRAM("One Factor", data = demoOneFactor[1:100, ],
#'  umxPath("G", to = manifests),
#'  umxPath(v.m. = manifests),
#'  umxPath(v1m0 = "G")
#' )
#' umxSummary(m1, std = TRUE, filter = "NS")
#' }
#'
umxSummary.MxModel <- function(model, refModels = NULL, std = FALSE, digits = 2, report = c("markdown", "html"), means= TRUE, residuals= TRUE, uncertainty = c("SE", "MLR", "CI", "none"), filter = c("ALL", "NS", "SIG"), RMSEA_CI = FALSE, SE = TRUE, matrixAddresses = FALSE, ...){
	# TODO make table take lists of models...
	if (!is.null(model$fitfunction) && is(model$fitfunction, "MxFitFunctionGLM")) {
		return(umxSummary.MxModelGLM(model, refModels = refModels, std = std, digits = digits, report = report, means = means, residuals = residuals, uncertainty = uncertainty, filter = filter, RMSEA_CI = RMSEA_CI, SE = SE, matrixAddresses = matrixAddresses, ...))
	}
	commaSep = paste0(umx_set_separator(silent = TRUE), " ")
	report   = match.arg(report)
	filter   = match.arg(filter)
	
	# Handle deprecated SE parameter and map it to uncertainty
	if (missing(uncertainty)) {
		if (!missing(SE)) {
			if (is.logical(SE)) {
				uncertainty = ifelse(SE, "SE", "none")
			} else {
				uncertainty = SE
			}
		} else {
			uncertainty = "SE"
		}
	} else {
		uncertainty = match.arg(uncertainty)
	}

	# RobustSE: engine-appropriate robust inference (not one sandwich for all fit types)
	isWLS = xmu_is_wls(model)
	if (uncertainty %in% c("MLR")) {
		if (isWLS) {
			# Plan A: keep OpenMx WLS/GMM moment-sandwich SEs (already in output$vcov).
			# Do not call imxRobustSE (casewise ML sandwich is wrong for WLS).
			if (is.null(model$output$vcov)) {
				warning("WLS model has no parameter covariance matrix in output; SEs unavailable. Run the model with standard errors enabled.", call. = FALSE)
			}
		} else {
			if (is.null(model$data) || is.null(model$data$observed) || identical(model$data$type, "cov")) {
				stop("Robust standard errors for ML require raw data. The model has covariance or missing data.")
			}
			resRobust = imxRobustSE(model, details = TRUE)
			model@output$vcov = resRobust$cov
			model@output$standardErrors = resRobust$SE
		}
	}

	SE = (uncertainty %in% c("SE", "MLR", "CI")) && !is.null(model$output$vcov)

	reportedWLS = FALSE
	
	# If the filter is not default, user must want something: Assume it's what would have been the default...
	if( filter != "ALL" & is.null(std) ) {
		std = FALSE
	}else if(!is.null(std)){
		 if(SE == FALSE && !is.null(model$output$vcov)){
			 # message("SE must be TRUE to show std, overriding to set SE = TRUE")
			 SE = TRUE
		 }
	}
	umx_has_been_run(model, stop = TRUE)

	if(is.null(refModels)) {
		# SaturatedModels not passed in from outside, so get them from the model
		modelSummary = summary(model)
		if(is.na(modelSummary$SaturatedLikelihood)){
			# no SaturatedLikelihood, compute refModels
			refModels = tryCatch({
				xmu_mxRefModels(model, run = TRUE, beginMessage = FALSE)
			}, error = function(x) {
			    print("Error calling mxRefModels: mxRefModels can't handle all designs https://github.com/OpenMx/OpenMx/issues/184")
			    NULL
			})

			if(!inherits(refModels, "list")){
				modelSummary = summary(model)
			} else {
				modelSummary = summary(model, refModels = refModels)
			}
		}
		if(is.null(model$data)){
			# TODO model with no data - no saturated solution?
			# message("Top model doesn't contain data. You may get extra information from summary() rather than umxSummary()")
		}
	} else if (is.logical(refModels) && !refModels){
		modelSummary = summary(model) # Don't use or generate refModels		
	}else{
		modelSummary = summary(model, refModels = refModels) # Using refModels supplied by user
	}

	if (isWLS) {
	    # GenomicMx / modern OpenMx: implied_jacobian present -> SB/robust WLS AFIs when computable
	    if (xmu_has_WLS_jacobian(model)) {
	        robustFit = NULL
	        tryCatch({
	            robustFit = xmu_robust_WLS_fit(model)
	        }, error = function(e) {
	            # Avoid raw R errors in user-facing notes
	            msg = conditionMessage(e)
	            if (grepl("replacement has length zero|Could not locate observed covariance|dimnames do not cover|Multigroup WLS robust", msg)) {
	            	# Multigroup or incomplete stack: quiet note once, not a stack of debug lines
	            	if (grepl("Multigroup WLS robust", msg)) {
	            		message("umxSummary Note: Multigroup WLS robust CFI/TLI/RMSEA not available yet (need stacked Jacobian across groups). Estimates and SEs are unaffected.")
	            	} else {
	            		message("umxSummary Note: robust CFI/TLI/RMSEA could not be computed; reporting unadjusted chi-square. Estimates and SEs are unaffected.")
	            	}
	            } else {
	            	message("umxSummary Note: robust CFI/TLI/RMSEA skipped: ", msg)
	            }
	        })
        
	        if (!is.null(robustFit)) {
	            # Patch incremental indices
	            modelSummary$CFI   = robustFit$CFI
	            modelSummary$TLI   = robustFit$TLI
	            modelSummary$RMSEA = robustFit$RMSEA
            
	            # Patch the base Chi-square statistics
	            modelSummary$Chi    = robustFit$Chi
	            modelSummary$ChiDoF = robustFit$ChiDoF
	            modelSummary$p      = robustFit$p

	            attr(modelSummary, "robustScalingFactors") = list(
	            	cModel = attr(robustFit, "c_model"),
	            	cNull = attr(robustFit, "c_null"),
	            	correction = attr(robustFit, "correction")
	            )
	        }
        
	    } else {
	        # Legacy catch: Warn the user (once per session), leave modelSummary unadjusted
	        if (is.null(getOption("umx_warned_legacy_wls")) || !getOption("umx_warned_legacy_wls")) {
	            warning(
	            	"Legacy OpenMx WLS engine detected (missing Jacobian). Fit statistics are unadjusted and unreliable.\n",
	            	xmu_openmx_install_message("WLS robust CFI/TLI/RMSEA and SB chi-square"),
	            	call. = FALSE
	            )
	            options(umx_warned_legacy_wls = TRUE)
	        }
	    }
	} else if (uncertainty %in% c("MLR")) {
		robustFit = NULL
		if (any(sapply(model$data$observed, is.ordered))) {
			message("Polite note: Robust SEs applied, but MLR scaling is not defined for ordinal models, so overall fit CFI/RMSEA.Chi are not scaled")
		} else {
			tryCatch({
				robustFit = xmu_robust_ML_fit(model, refModels = refModels)
			}, error = function(e) {
				message(paste("umxSummary Note: Frontend robust ML calculation skipped:", e$message))
			})
		}

		if (!is.null(robustFit)) {
			modelSummary$CFI   = robustFit$CFI
			modelSummary$TLI   = robustFit$TLI
			modelSummary$RMSEA = robustFit$RMSEA

			modelSummary$Chi    = robustFit$Chi
			modelSummary$ChiDoF = robustFit$ChiDoF
			modelSummary$p      = robustFit$p
			
			attr(modelSummary, "robustScalingFactors") = list(
				cModel = robustFit$scalingFactor,
				cNull = robustFit$scalingFactorNull
			)
		}
	}
	
	# DisplayColumns show
	if(!is.null(std)){
		# n.b: mxStandardizeRAMpaths returns the raw paths as well, so two birds, one stone.
		parameterTable = mxStandardizeRAMpaths(model, SE = SE) # Compute standard errors
		
		# Helper to recursively flatten lists of dataframes returned by mxStandardizeRAMpaths
		flatten_pTable <- function(x) {
			if (is.data.frame(x)) {
				return(x)
			} else if (is.list(x)) {
				if (length(x) == 0) return(NULL)
				dfs = lapply(x, flatten_pTable)
				dfs = dfs[vapply(dfs, is.data.frame, logical(1))]
				if (length(dfs) > 0) {
					return(do.call(rbind, dfs))
				}
			}
			return(NULL)
		}
		
		parameterTable = flatten_pTable(parameterTable)

		if (!is.null(parameterTable) && nrow(parameterTable) > 0) {
			#          name    label  matrix   row         col    Raw.Value  Raw.SE   Std.Value    Std.SE
			# 1  Dep.A[6,1]    age    A        mean_sdrr   age   -0.37       0.0284   -0.372350    .028
			# Raw.SE is new
			names(parameterTable) = c("label", "name", "matrix", "row", "col", "Estimate", "SE", "Std.Estimate", "Std.SE")

			if(matrixAddresses){
				naming = c("name", "matrix", "row", "col")
			} else {
				naming = c("name")
			}
			# TODO: umxSummary add p value, perhaps CI?
			
			if(std == TRUE){
				if (uncertainty == "none") {
					namesToShow = c(naming, "Std.Estimate")
				} else if (uncertainty == "CI") {
					namesToShow = c(naming, "Std.Estimate", "CI")
				} else {
					namesToShow = c(naming, "Std.Estimate", "Std.SE", "CI")
				}
			}else{ # must be raw
				if (uncertainty == "none") {
					namesToShow = c(naming, "Estimate")
				} else if (uncertainty == "CI") {
					namesToShow = c(naming, "Estimate", "CI")
				} else {
					namesToShow = c(naming, "Estimate", "SE", "CI")
				}
			}

			parameterTable$sig = TRUE
			if("CI" %in% namesToShow){
				parameterTable$CI  = ""
				cis = model$output$confidenceIntervals
				hasCIs = !is.null(cis) && (nrow(cis) > 0)
				for(i in 1:dim(parameterTable)[1]) {
					pName = parameterTable[i, "name"]
					est   = parameterTable[i, ifelse(std, "Std.Estimate", "Estimate")]
					
					if (uncertainty == "CI" && hasCIs && (pName %in% rownames(cis))) {
						lbound = cis[pName, "lbound"]
						ubound = cis[pName, "ubound"]
						parameterTable[i, "CI"] = paste0(round(est, digits), " [", round(lbound, digits), commaSep, round(ubound, digits), "]")
						if (!is.na(lbound) && !is.na(ubound)) {
							if (lbound <= 0 && ubound >= 0) {
								parameterTable[i, "sig"] = FALSE
							}
						}
					} else {
						seVal = suppressWarnings(as.numeric(parameterTable[i, ifelse(std, "Std.SE", "SE")]))
						CI95  = seVal * 1.96
						bounds = c(est - CI95, est + CI95)

						if(any(is.na(bounds))) {
							# protect cases with SE == NA from evaluation for significance
						} else {
							if (any(bounds <= 0) & any(bounds >= 0)){
								parameterTable[i, "sig"] = FALSE
							}
							parameterTable[i, "CI"] = paste0(round(est, digits), " [", round(est - CI95, digits), commaSep, round(est + CI95, digits), "]")
						}
					}
				}
			}
			if(filter == "NS") {
				toShow = parameterTable[parameterTable$sig == FALSE, namesToShow]
			} else if(filter == "SIG") {
				toShow = parameterTable[parameterTable$sig == TRUE, namesToShow]
			} else {
				toShow = parameterTable[, namesToShow]
			}
			toShow = xmu_summary_RAM_group_parameters(model, toShow,  means= means, residuals = residuals)
			toShow = unique.data.frame(toShow[,c(namesToShow, "type")])
			umx_print(toShow, digits = digits, report = report, caption = paste0("Parameter loadings for model ", omxQuotes(model$name)), na.print = "", zero.print = "0", justify = "none")
		}
	}
	if (!reportedWLS) {
		hasFitIndices = !is.null(modelSummary$TLI) && !is.null(modelSummary$CFI) && !is.null(modelSummary$RMSEA) && !is.null(modelSummary$Chi)
		if(hasFitIndices) {
			with(modelSummary, {
				if(!is.finite(TLI)){
					TLI_OK = "OK"
				} else {
					if(TLI > .95) {
						TLI_OK = "OK"
					} else {
						TLI_OK = "bad"
					}
				}
				if(!is.finite(RMSEA)) {
					RMSEA_OK = "OK"
				} else {
					if(RMSEA < .06){
					RMSEA_OK = "OK"
					} else {
						RMSEA_OK = "bad"
					}
				}
				if(report == "table"){
					x = data.frame(cbind(model$name, round(Chi,2), formatC(p, format="g"), round(CFI,3), round(TLI,3), round(RMSEA, 3)))
					names(x) = c("model","\u03C7","p","CFI", "TLI","RMSEA") # \u03A7 is unicode for chi
					print(x)
				} else {
					if(RMSEA_CI){
						RMSEA_CI = RMSEA(modelSummary)$txt
					} else {
						RMSEA_CI = paste0("RMSEA = ", round(RMSEA, 3))
					}
					fitMsg = paste0("\nModel Fit: \u03C7\u00B2(", ChiDoF, ") = ", round(Chi, 2), # was A7
						# "Chi2(", ChiDoF, ") = ", round(Chi, 2), # was A7
						", p "      , umx_APA_pval(p, .001, 3, addComparison = TRUE),
						"; CFI = "  , round(CFI, 3),
						"; TLI = "  , round(TLI, 3),
						"; ", RMSEA_CI
					)
					message(fitMsg)
					if (xmu_is_wls(model)) {
						robustScaling = attr(modelSummary, "robustScalingFactors")
						if (umx_is_GSEM(model)) {
							cat("\n*Statistical Note*: Genomic SEM - absolute fit: prefer SRMR (roughly < 0.10) and residual inspection. Nested models: use umxCompare (DWLS chi-square difference). De-emphasize CFI/TLI/RMSEA; do not apply Hu-Bentler cutoffs. See ?umxCompare.\n")
						} else if (!is.null(robustScaling) && identical(robustScaling$correction, "Savalei2021")) {
							cat(sprintf("\n*Statistical Note*: Ordinal WLS robust fit (Savalei (2021) categorical ML scaling (c_model = %.3f, c_null = %.3f). Unlike unscaled WLSMV which lacks this calibration, Standard Hu & Bentler (1999) cutoffs apply (CFI > .95, RMSEA < .06).\n", robustScaling$cModel, robustScaling$cNull))
						} else {
							cat("\n*Statistical Note*: Continuous WLS/DWLS - conventional CFI/TLI/RMSEA cutoffs do not apply. Prefer SRMR for absolute fit; nested comparisons use Strict Satorra-Bentler (2010) Delta chi-square via umxCompare. See ?umxCompare.\n")
						}
						if (uncertainty %in% c("MLR")) {
							cat("*SEs*: WLS robust SEs (sandwich estimator using asymptotic covariance of summary statistics).\n")
						}
					} else {
						robustScaling = attr(modelSummary, "robustScalingFactors")
						if (!is.null(robustScaling)) {
							isFiml = FALSE
							if (!is.null(model$data) && model$data$type == "raw" && !is.null(model$data$observed)) {
								manifests = model$manifestVars
								if (is.null(manifests) || length(manifests) == 0) {
									manifests = colnames(model$data$observed)
								}
								isFiml = any(is.na(model$data$observed[, manifests, drop = FALSE]))
							}
							if (isFiml) {
								cat(sprintf("\n*Statistical Note*: FIML robust indices use Yuan-Bentler trace scaling + Brosseau-Liard/Savalei adjustment (c_model = %.3f, c_null = %.3f).\n", robustScaling$cModel, robustScaling$cNull))
							} else {
								cat(sprintf("\n*Statistical Note*: Robust ML indices use sandwich-based scaling + Brosseau-Liard & Savalei (2012) adjustments to CFI/TLI/RMSEA (c_model = %.3f, c_null = %.3f).\n", robustScaling$cModel, robustScaling$cNull))
							}
						} else {
							if(TLI_OK   != "OK"){ message("TLI is worse than desired (>.95)") }
							if(RMSEA_OK != "OK"){ message("RMSEA is worse than desired (<.06)")}
						}
					}
				}
			})
		} else {
			# Fallback if fit indices are not available
			minus2LL = modelSummary$Minus2LogLikelihood
			df = modelSummary$degreesOfFreedom
			estimatedParameters = modelSummary$estimatedParameters
			
			# If estimatedParameters or minus2LL is not NULL/NA, calculate AIC
			if (!is.null(minus2LL) && !is.na(minus2LL) && !is.null(estimatedParameters)) {
				aic = minus2LL + 2 * estimatedParameters
				fitMsg = paste0("\nModel Fit: -2LL = ", round(minus2LL, 2), 
				                ", df = ", df, 
				                ", AIC = ", round(aic, 2))
			} else {
				fitMsg = "\nModel Fit: Fit statistics not available (model may not have run successfully)."
			}
			message(fitMsg)
		}
	}
	# TODO: umxSummary.MxRAMModel integrate interval printing into summary table
	if(!is.null(model$output$confidenceIntervals)){
		print(model$output$confidenceIntervals)
	}
	
	xmu_print_algebras(model)
	if(!is.null(std) && !is.null(parameterTable) && nrow(parameterTable) > 0){ # return these as  invisible for the user to filter, sort etc.
		robustScaling = attr(modelSummary, "robustScalingFactors")
		if (!is.null(robustScaling)) {
			attr(parameterTable, "robustScalingFactors") = robustScaling
			attr(parameterTable, "c_model") = robustScaling$cModel
			attr(parameterTable, "c_null") = robustScaling$cNull
		}
		if(filter == "NS"){
			invisible(parameterTable[parameterTable$sig == FALSE, namesToShow])
		}else if(filter == "SIG"){
			invisible(parameterTable[parameterTable$sig == TRUE, namesToShow])
		}else{
			invisible(parameterTable[,namesToShow])
		}
	} else {
		invisible(NULL)
	}
}

#' Summary for a RAM model with mxFitFunctionGLM
#'
#' Parameter table and -2LL / AIC / BIC. No CFI/TLI/RMSEA (no saturated GLM
#' pair yet). Nested tests: [umxCompare()] / AIC. Called from [umxSummary()]
#' when the fit function is GLM; not a new S4 class.
#'
#' @param model an [OpenMx::mxModel()] with `OpenMx::mxFitFunctionGLM()`
#' @param refModels ignored (no GLM saturated model yet)
#' @param std ignored (the box is eta; standardization is not defined)
#' @param digits rounding
#' @param report "markdown" or "html"
#' @param means unused (kept for umxSummary signature)
#' @param residuals unused
#' @param uncertainty "SE" (default) or "none"
#' @param filter unused
#' @param RMSEA_CI unused
#' @param SE deprecated; use uncertainty
#' @param matrixAddresses unused
#' @param ... unused
#' @return parameter table (invisibly)
#' @exportS3Method umxSummary MxModelGLM
#' @md
umxSummary.MxModelGLM <- function(model, refModels = NULL, std = FALSE, digits = 2, report = c("markdown", "html"), means = TRUE, residuals = TRUE, uncertainty = c("SE", "MLR", "CI", "none"), filter = c("ALL", "NS", "SIG"), RMSEA_CI = FALSE, SE = TRUE, matrixAddresses = FALSE, ...) {
	report = match.arg(report)
	if (missing(uncertainty)) {
		if (!missing(SE) && is.logical(SE) && !SE) {
			uncertainty = "none"
		} else {
			uncertainty = "SE"
		}
	} else {
		uncertainty = match.arg(uncertainty)
	}
	umx_has_been_run(model, stop = TRUE)
	if (isTRUE(std)) {
		message("Polite note: standardized paths are not defined for GLM (the box is eta). Showing raw estimates.")
	}
	if (!is.null(refModels) && !isFALSE(refModels)) {
		message("Polite note: GLM has no saturated/independence pair yet; ignoring refModels. Use umxCompare or AIC.")
	}
	modelSummary = summary(model)
	pars = modelSummary$parameters
	if (!is.null(pars) && nrow(pars) > 0) {
		tab = data.frame(
			name = pars$name,
			Estimate = pars$Estimate,
			SE = if ("Std.Error" %in% names(pars)) pars$Std.Error else NA_real_,
			stringsAsFactors = FALSE
		)
		if (uncertainty == "none") {
			umx_print(tab[, c("name", "Estimate")], digits = digits, report = report, caption = paste0("Parameters for GLM model ", omxQuotes(model$name)), na.print = "")
		} else {
			umx_print(tab, digits = digits, report = report, caption = paste0("Parameters for GLM model ", omxQuotes(model$name)), na.print = "")
		}
	}
	minus2LL = modelSummary$Minus2LogLikelihood
	nPar = modelSummary$estimatedParameters
	aic = modelSummary$AIC
	bic = modelSummary$BIC
	nObs = NA
	if (!is.null(model$data) && !is.null(model$data$numObs)) nObs = model$data$numObs
	status = NA
	if (!is.null(model$output$status$code)) status = model$output$status$code
	fitBits = c()
	if (!is.null(minus2LL) && is.finite(minus2LL)) fitBits = c(fitBits, paste0("-2LL = ", round(minus2LL, digits)))
	if (!is.null(nPar)) fitBits = c(fitBits, paste0("k = ", nPar))
	if (!is.null(aic) && is.finite(aic)) fitBits = c(fitBits, paste0("AIC = ", round(aic, digits)))
	if (!is.null(bic) && is.finite(bic)) fitBits = c(fitBits, paste0("BIC = ", round(bic, digits)))
	if (is.finite(nObs)) fitBits = c(fitBits, paste0("n = ", nObs))
	if (length(status) == 1 && !is.na(status)) fitBits = c(fitBits, paste0("status = ", status))
	if (length(fitBits) > 0) {
		message(paste0("\nGLM fit: ", paste(fitBits, collapse = "; ")))
	}
	predNames = c()
	if (!is.null(model$S) && !is.null(dimnames(model$S$values))) {
		sNames = dimnames(model$S$values)[[1]]
		famN = c()
		if (is(model$fitfunction, "MxFitFunctionGLM") && length(model$fitfunction$family) > 0) {
			famN = names(model$fitfunction$family)
		}
		for (nm in sNames) {
			if (nm %in% famN) next
			if (!isTRUE(model$S$free[nm, nm]) && abs(model$S$values[nm, nm]) < 1e-12) {
				predNames = c(predNames, nm)
			}
		}
	}
	if (length(predNames) > 0) {
		message(paste0("Observed predictors (no mean/variance; arrows use the data column): ", paste(predNames, collapse = ", ")))
	}
	message("No CFI/TLI/RMSEA for GLM yet. Nested models: umxCompare / AIC.")
	if (!is.null(pars) && nrow(pars) > 0) {
		invisible(pars)
	} else {
		invisible(NULL)
	}
}

#' @export
umxSummary.MxRAMModel <- umxSummary.MxModel

#' @export
umxSummary.MxModelGSEM = umxSummary.MxModel