R/umxCompare.R

In umx: Structural Equation Modeling and Twin Modeling in R

#' Print a comparison table of two or more [OpenMx::mxModel()]s, nicely formatted.
#'
#' @description
#' umxCompare compares two or more [OpenMx::mxModel()]s. It has several nice features:
#' 
#' 1. Direct control of precision via digits=. p-values use APA style.
#' 2. Report = directs publication to your chosen format: markdown, html, latex.
#' 3. Columns arranged to make for easy visual comparison.
#' 4. report = 'inline', includes a summary sentence suitable for your report.
#' 5. report = "html" opens a web table in your browser to paste into a word processor.
#' 
#'
#' **Interpreting fit under WLS/DWLS and Genomic SEM**
#'
#' Hu & Bentler (1999) style cutoffs (e.g. CFI > .95, RMSEA < .06) were developed for
#' maximum likelihood with continuous data. They do **not** transfer to WLS/DWLS
#' (Shi et al., 2020), and still less to Genomic SEM (DWLS on LDSC genetic
#' covariances with an estimated sampling covariance, often with bookkeeping `numObs = 1`).
#'
#' * **Incremental indices (CFI, TLI, NFI, ...)** compare the model to an independence
#'   baseline under the same weight matrix. That baseline behaves differently than under ML.
#'   Report CFI/TLI only **descriptively**; do **not** apply conventional cutoffs.
#' * **Absolute fit: SRMR** (and residual matrices) is preferred. SRMR lives on a
#'   standardized residual metric and is less hostage to weight-matrix / N scaling.
#'   Rough guide: SRMR < 0.10 suggests acceptable residual structure (not a hard law).
#' * **RMSEA** folds in chi-square scaling and N; it is often misleading for WLS and
#'   especially for GSEM. Do not use the ML cutoff of < .06.
#' * **Nested model building:** prefer scaled chi-square difference tests and DeltaSRMR
#'   over DeltaCFI alone (see Details for what `umxCompare` computes).
#'
#' @details
#' **What `umxCompare` does for WLS**
#'
#' * Refuses mixed WLS vs ML comparisons (same engine required).
#' * **Continuous WLS** with cached `implied_jacobian`:
#'   - Table **Chi** is the same Satorra-Bentler (2010) scaled omnibus as [umxSummary()]
#'     (\eqn{F/c} from `output$fit`), **not** OpenMx Browne residual `output$chi`
#'     (those two differ under DWLS). Saturated residual df reports Chi = 0, CFI = 1, RMSEA = 0.
#'   - Nested **diffFit** is SB-2010 scaled \eqn{\Delta F = F_{nested} - F_{base}}
#'     (requires both models at comparable minima of the same WLS objective: same moments,
#'     `useWeight`, and \eqn{N}). If \eqn{F_{nested} < F_{base}} (often multigroup
#'     optimization), diffFit is **NA** with a warning - not a negative chi-square with p = 1.
#' * **Genomic SEM** (`MxModelGSEM`): nested difference on the **GSEM DWLS chi-square
#'   scale** (same discrepancy as the structural fit; not a substitute for ML LRT).
#'   If LDSC matrices were `nearPD`-smoothed, a fiduciary warning notes that difference
#'   tests may look artificially precise.
#' * Table also reports **SRMR / DeltaSRMR** (preferred absolute residual summary) and
#'   **CFI / DeltaCFI** (descriptive only). AIC may appear but is not a primary GSEM decision
#'   rule under asymptotic / bookkeeping N.
#'
#' **Ordinal WLS** is a separate track: robust CFI/TLI/RMSEA use Savalei (2021) catML
#' corrections in [umxSummary()]; Hu-Bentler cutoffs can apply to **those** robust indices.
#' That exception does **not** apply to continuous WLS or Genomic SEM.
#'
#' Best practice: report estimates with SEs, SRMR (and residuals), nested SB/GSEM
#' difference tests, and avoid single-number "good fit" claims from CFI or RMSEA.
#'
#' **Constraints and the EP column**
#'
#' EP is OpenMx's count of free parameters, `summary(model)$estimatedParameters`.
#' Equality [OpenMx::mxConstraint()]s are not removed from that count. OpenMx adds
#' their rows to `observedStatistics`, and degrees of freedom already include them.
#' A common-pathway model that sets each factor with `A + C + E = 1` can therefore
#' show a larger EP than an independent-pathway model while \eqn{\Delta}{Delta} df is 0. Both
#' numbers are the OpenMx counts. When `sum(summary(model)$constraints)` is above 0,
#' `umxCompare` adds one line for that model:
#'
#' `Note: df for 'CP' reflects 3 constraints: 51 estimated parameters (48 effective after constraints).`
#'
#' With `silent = TRUE`, the lines are `attr(table, "constraintNotes")`, one element
#' per model. A comparison in which no model has an equality constraint does not
#' set that attribute. Inequality constraints contribute 0 and do not produce a line.
#'
#' @param base The base [OpenMx::mxModel()] for comparison
#' @param comparison The model (or list of models) which will be compared for fit with the base model (can be empty)
#' @param all Whether to make all possible comparisons if there is more than one base model (defaults to T)
#' @param digits rounding for p-values etc.
#' @param report "markdown" (default), "inline" (a sentence suitable for inclusion in a paper), or "html".
#' create a web table and open your default browser.
#' (handy for getting tables into Word, and other text systems!)
#' @param file file to write html too if report = "html" (defaults to "tmp.html")
#' @param compareWeightedAIC Show the Wagenmakers AIC weighted comparison (default = FALSE)
#' @param silent (don't print, just return the table as a dataframe (default = FALSE)
#' @param uncertainty What type of parameter uncertainty to report: "SE" (standard ML standard errors), "MLR" (robust standard errors and robust fit), "CI" (profile likelihood confidence intervals), or "none" (none).
#' @family Model Summary and Comparison
#' @seealso - [umxSummary()], [umxRAM()],[umxCompare()]
#' @references <https://github.com/tbates/umx>
#' @export
#' @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)
#' )
#'
#' m2 = umxModify(m1, update = "G_to_x2", name = "drop_path_2_x2")
#' umxCompare(m1, m2)
#' umxCompare(m1, m2, report = "inline") # Add English-sentence descriptions
#' umxCompare(m1, m2, report = "html") # Open table in browser
#'
#' # Two comparison models
#' m3 = umxModify(m2, update = "G_to_x3", name = "drop_path_2_x2_and_3")
#' 
#' umxCompare(m1, c(m2, m3))
#' umxCompare(m1, c(m2, m3), compareWeightedAIC = TRUE)
#' umxCompare(c(m1, m2), c(m2, m3), all = TRUE)
#'
#' manifests = names(demoOneFactor)
#' m1 = umxRAM("WLS", data = demoOneFactor, type = "DWLS",
#' 	umxPath("G", to = manifests),
#' 	umxPath(var = manifests),
#' 	umxPath(var = "G", fixedAt = 1)
#' )
#'
#' m2 = umxModify(m1, update = "G_to_x2", name = "drop_path_2_x2")
#' umxCompare(m1, m2)
#' umxCompare(m1, m2, report = "inline") # Add English-sentence descriptions
#' umxCompare(m1, m2, report = "html") # Open table in browser
#' }
umxCompare <- function(base = NULL, comparison = NULL, all = TRUE, digits = 3, report = c("markdown", "html", "inline"), compareWeightedAIC = FALSE, silent = FALSE, file = "tmp.html", uncertainty = c("none", "SE", "MLR", "CI")) {
	report = match.arg(report)
	uncertainty = match.arg(uncertainty)
	if(umx_is_MxModel(all)){
		stop("Provide all comparison models as a c() (You provided a model as input to 'all', and I'm guessing that's a mistake)")
	}
	if(is.null(comparison)){
		comparison = base
	} else if (is.null(base)) {
		stop("You must provide at least a base model for umxCompare")
	}
	if(length(base) == 1) {
		if(typeof(base) == "list"){
			base = base[[1]]
		}
		if(!umx_has_been_run(base)){
			warning("Base model not run yet!")		
		}
	}
	# A comparison list is already the models to walk. umx_is_MxModel()
	# warns, and returns FALSE, when handed that list.
	if (!is.list(comparison) && umx_is_MxModel(comparison)) {
		comparison = list(comparison)
	}
	if (!is.list(comparison)) {
		stop("comparison must be a list or an MxModel")
	}
	for (comp in comparison) {
		if (!umx_has_been_run(comp)) {
			stop("Comparison model has not been run!")
		}
	}

	# A base list longer than one is the models themselves. xmu_is_wls()
	# calls umx_is_MxModel(), which warns and returns FALSE on that list.
	baseModels = base
	if (!is.list(baseModels)) {
		baseModels = list(base)
	}
	nWlsBase = 0
	for (b in baseModels) {
		if (xmu_is_wls(b)) {
			nWlsBase = nWlsBase + 1
		}
	}
	if (nWlsBase > 0 && nWlsBase != length(baseModels)) {
		stop("Engine Mismatch: Cannot compare a WLS model with an ML model.")
	}
	baseIsWLS = nWlsBase > 0
	baseHasJac = xmu_has_WLS_jacobian(baseModels[[1]])

	if (baseIsWLS) {
		# 1. Enforce Homogeneity Constraint across all base and comparison models
		for (b in baseModels) {
			if (baseHasJac != xmu_has_WLS_jacobian(b)) {
				stop("Engine Mismatch: Cannot compare a legacy OpenMx WLS model with a GenomicMx WLS model. Both models must use the same engine.")
			}
		}
		for (comp in comparison) {
			if (!xmu_is_wls(comp)) {
				stop("Engine Mismatch: Cannot compare a WLS model with an ML model.")
			}
			if (baseHasJac != xmu_has_WLS_jacobian(comp)) {
				stop("Engine Mismatch: Cannot compare a legacy OpenMx WLS model with a GenomicMx WLS model. Both models must use the same engine.")
			}
		}

		# 1. Initialize an empty data frame to hold the master table
		finalTable = data.frame()

		# One base: every comparison against that base (previous behavior).
		# Several bases: all=TRUE is every pair; all=FALSE cycles, as mxCompare does.
		pairBase = list()
		pairComp = list()
		if (length(baseModels) == 1 || all) {
			for (b in baseModels) {
				for (comp in comparison) {
					pairBase[[length(pairBase) + 1]] = b
					pairComp[[length(pairComp) + 1]] = comp
				}
			}
		} else {
			maxLength = max(length(baseModels), length(comparison))
			for (i in seq_len(maxLength)) {
				pairBase[[i]] = baseModels[[(i - 1) %% length(baseModels) + 1]]
				pairComp[[i]] = comparison[[(i - 1) %% length(comparison) + 1]]
			}
		}

		constraintNotes = xmu_compare_constraint_note(c(baseModels, comparison))

		# 2. Iterate explicitly through each base vs. comparison pair
		for (i in seq_along(pairComp)) {
			compTable = xmu_compare_WLS(baseModel = pairBase[[i]], comparisonModel = pairComp[[i]])
			if (nrow(finalTable) == 0) {
				finalTable = compTable
			} else {
				finalTable = rbind(finalTable, compTable)
			}
		}

		# TODO match the output as much as possible with what we do for ML (see further down)
		# names(tablePub) = c("Model", "EP", "\u0394 Fit" , "\u0394 df" , "p", "AIC", "\u0394 AIC", "Compare with Model", "Fit units")
		if (!silent) {
			umx_print(finalTable, digits = digits, zero.print = "0", caption = "Table of Model Comparisons", report = report)
			units_str = summary(baseModels[[1]])$fitUnits
			if (is.null(units_str) || length(units_str) == 0) {
				units_str = "r'wr"
			}
			isGenomic = FALSE
			for (b in baseModels) {
				if (umx_is_GSEM(b)) {
					isGenomic = TRUE
				}
			}
			for (comp in comparison) {
				if (umx_is_GSEM(comp)) {
					isGenomic = TRUE
				}
			}

			if (isGenomic) {
				cat("\n*Note* EP: Estimated parameters; Chi: Browne (1984) residual chi-square; SRMR: Standardized Root Mean Square Residual; CFI: Comparative Fit Index; AIC: Akaike Information Criterion; diffFit: Browne (1984) residual chi-square difference test. See ?umxCompare for details.\n")
			} else {
				cat("\n*Note* EP: Estimated parameters; Chi: Satorra-Bentler (2010) scaled WLS chi-square; SRMR: Standardized Root Mean Square Residual; CFI: Comparative Fit Index; AIC: Akaike Information Criterion; diffFit: Satorra-Bentler (2010) scaled chi-square difference test. See ?umxCompare for details.\n")
			}
			for (constraintNote in constraintNotes) {
				cat("\n", constraintNote, "\n", sep = "")
			}
		}
		if (length(constraintNotes) > 0) {
			attr(finalTable, "constraintNotes") = constraintNotes
		}
		return(invisible(finalTable))
	} else {
		# Base is ML. Ensure no comparison models are WLS to prevent reverse-mismatch.
		for (comp in comparison) {
			if (xmu_is_wls(comp)) {
				stop("Engine Mismatch: Cannot compare an ML model with a WLS model.")
			}
		}
	}

	robustScalingFactors = list()
	tableOut = mxCompare(base = base, comparison = comparison, all = all)
	tableOut = as.data.frame(tableOut)

	tablePub = tableOut[, c("comparison", "ep", "diffFit", "diffdf", "p", "AIC", "base", "fitUnits")]

	# Subtract row-1 AIC from all values and place the resulting deltaAIC column after AIC 
	tablePub$deltaAIC = tablePub[, "AIC"] - tablePub[1, "AIC"]
	tablePub = tablePub[,c("comparison", "ep", "diffFit", "diffdf", "p", "AIC", "deltaAIC", "base", "fitUnits")]

	# Compute and patch robust Satorra-Bentler difference statistics if requested
	if (uncertainty %in% c("MLR")) {
		for (comp in comparison) {
			resSb = xmu_compare_robust_ML(base, comp)
			if (!is.null(resSb)) {
				compName = comp$name
				rowIdx = which(tablePub$comparison == compName)
				if (length(rowIdx) == 1) {
					tablePub[rowIdx, "diffFit"] = resSb$diffFit
					tablePub[rowIdx, "p"] = resSb$p
					robustScalingFactors[[compName]] = resSb$scalingFactor
				}
			}
		}
	}

	# c("1: Comparison", "2: Base", "3: EP", "4: AIC", "5: &Delta; -2LL", "6: &Delta; df", "7: p")
	# U+2206 = math delta
	# Fix problem where base model has compare set to its own name, and name set to NA
	nRows = dim(tablePub)[1]
	for (i in 1:nRows) {
		if(is.na(tablePub[i, "comparison"])){
			tablePub[i, "comparison"] = tablePub[i, "base"]
			tablePub[i, "base"] = NA
		}
	}
	tablePub[, "p"] = umx_APA_pval(tablePub[, "p"], min = (1/ 10^3), digits = digits, addComparison = NA)
	if(report == "inline"){
		n_rows = dim(tablePub)[1]
		for (i in 1:n_rows) {
			thisPValue = tableOut[i, "p"]
			if(!is.na(thisPValue) && !is.nan(thisPValue)){
				if(tableOut[i, "p"] < .05){
					this = ". This caused a significant loss of fit "
				} else {
					this = ". This did not lower fit significantly "
				}
				inlineMsg = paste0("The hypothesis that ", omxQuotes(tablePub[i, "comparison"]), 
				" was tested by dropping ", tablePub[i, "comparison"],
				" from ", omxQuotes(tablePub[i, "base"]), 
				this, "(\u03C7\u00B2(", tablePub[i, "diffdf"], ") = ", round(tablePub[i, "diffFit"], 2), # \u03A7 = Chi \u00B2 = superscript 2
				", p = ", tablePub[i, "p"], ": AIC = ", round(tablePub[i, "AIC"], digits), " change in AIC = ", round(tablePub[i, "deltaAIC"], digits), ")."
				)
				if(!silent){
					cat(inlineMsg)
				}
			}
		}
	}
	
	# Rename for printing
	names(tablePub) = c("Model", "EP", "\u0394 Fit" , "\u0394 df" , "p", "AIC", "\u0394 AIC", "Compare with Model", "Fit units")

	if(report == "inline"){ report= "markdown"}
	
	# Check if any model in the comparison is a WLS model
	anyWLS = FALSE
	checkModels = list()
	if (is.list(base)) {
		checkModels = c(checkModels, base)
	} else if (!is.null(base)) {
		checkModels = c(checkModels, list(base))
	}
	if (is.list(comparison)) {
		checkModels = c(checkModels, comparison)
	} else if (!is.null(comparison)) {
		checkModels = c(checkModels, list(comparison))
	}
	
	for (m in checkModels) {
		if (umx_is_MxModel(m) && xmu_is_wls(m)) {
			anyWLS = TRUE
			break
		}
	}

	noteModels = comparison
	if (is.list(base)) {
		noteModels = c(base, noteModels)
	} else if (!is.null(base)) {
		noteModels = c(list(base), noteModels)
	}
	constraintNotes = xmu_compare_constraint_note(noteModels)

	if(!silent){
		umx_print(tablePub, digits = digits, zero.print = "0", caption = "Table of Model Comparisons", report = report)
		if (anyWLS) {
			units_val = unique(tablePub$`Fit units`[!is.na(tablePub$`Fit units`)])
			if (length(units_val) == 0) {
				units_str = "WLS discrepancy function"
			} else {
				units_str = paste(units_val, collapse = ", ")
			}
			cat(paste0("\n*Note*: EP = Estimated (i.e. free) parameters; \u0394 Fit = change in fit (units: ", units_str, "); \u0394 df = Change in degrees of freedom with respect to the comparison model; \u0394 AIC = Change in Akaike Information Criterion; 'Compared to' = The baseline model for this comparison.\n"))
			cat("\n*Statistical Note*: WLS/GSEM - conventional CFI/TLI/RMSEA cutoffs do not apply. Prefer SRMR for absolute fit; nested comparisons use scaled Delta chi-square (diffFit). De-emphasize incremental indices. See ?umxCompare.\n")

		} else {
			if (length(robustScalingFactors) > 0) {
				factorsStr = paste(sapply(names(robustScalingFactors), function(n) sprintf("%s (c_d = %.3f)", n, robustScalingFactors[[n]])), collapse = ", ")
				cat(sprintf("\n*Statistical Note*: Robust ML difference tests use Satorra-Bentler (2001) scaled difference formula: %s.\n", factorsStr))
			} else {
				cat("\n*Note*: EP = Estimated (i.e. free) parameters; \u0394-2LL = change in -2 \u00D7 Log-Likelihood of the model; \u0394 df = Change in degrees of freedom with respect to the comparison model; \u0394 AIC = Change in Akaike Information Criterion; 'Compared to' = The baseline model for this comparison.\n")
			}
		}
		for (constraintNote in constraintNotes) {
			cat("\n", constraintNote, "\n", sep = "")
		}
	}

	if(compareWeightedAIC){
		modelList = c(base, comparison)
		# get list of AICs
		AIClist = c()
		for (i in modelList) {
			AIClist = c(AIClist, AIC(i))
		}
		whichBest = which.min(AIClist)
		bestModel = modelList[[whichBest]]
		aic.weights = round(aic_weights(AIClist), 2)
		if(!silent){
			cat("The ", omxQuotes(bestModel$name), " model is the best fitting model according to AIC.")
			cat("AIC weight-based  {Wagenmakers, 2004, 192-196} conditional probabilities of being the best model for ", 
				omxQuotes(namez(modelList)), " respectively are: ", 
				omxQuotes(aic.weights), " Using MuMIn::Weights(AIC()).")	
		}
		
	}
	if (length(robustScalingFactors) > 0) {
		attr(tablePub, "robustScalingFactors") = robustScalingFactors
		attr(tablePub, "c_d") = robustScalingFactors
	}
	if (length(constraintNotes) > 0) {
		attr(tablePub, "constraintNotes") = constraintNotes
	}
	invisible(tablePub)
}

#' Lines reconciling raw free parameters with equality constraints.
#'
#' OpenMx stores equality [OpenMx::mxConstraint()] rows in `summary(model)$constraints`
#' and adds that total to `observedStatistics`. `estimatedParameters` stays the raw
#' free-parameter count, which is also the EP column of [umxCompare()].
#'
#' @param models A list of [OpenMx::mxModel()] objects.
#' @return A character vector, one line per model whose equality-constraint row count
#' is not 0. An empty character vector when no model qualifies. Unrun models are skipped.
#' @family xmu internal not for end user
xmu_compare_constraint_note <- function(models) {
	# Equality constraints are extra observed statistics, not fewer free parameters.
	# sum(constraints) is that row count. Inequalities contribute 0.
	notes = character()
	seen = character()
	for (m in models) {
		if (!umx_is_MxModel(m)) {
			next
		}
		if (!umx_has_been_run(m)) {
			next
		}
		modelName = m$name
		if (is.null(modelName) || is.na(modelName) || modelName %in% seen) {
			next
		}
		seen = c(seen, modelName)
		modelSummary = summary(m)
		constraintRows = modelSummary$constraints
		if (is.null(constraintRows) || length(constraintRows) == 0) {
			next
		}
		nConstraints = sum(constraintRows)
		epLive = modelSummary$estimatedParameters
		if (!is.finite(nConstraints) || !is.finite(epLive) || nConstraints == 0) {
			next
		}
		if (abs(nConstraints) == 1) {
			constraintWord = "constraint"
		} else {
			constraintWord = "constraints"
		}
		epEffective = epLive - nConstraints
		notes = c(notes, paste0(
			"Note: df for '", modelName, "' reflects ", nConstraints, " ", constraintWord, ": ",
			epLive, " estimated parameters (", epEffective, " effective after constraints)."
		))
	}
	notes
}