R/umx_fit_umxFitIndices.R
In tbates/umx: Structural Equation Modeling and Twin Modeling in R
Defines functions
print.umxFitIndices
xmuResidualIndices
xmuResiduals
xmuModelMatrices
umxFitIndices
Documented in
umxFitIndices
#' Get additional fit-indices for a model with umxFitIndices
#'
#' Computes a variety of fit indices. Originated in this thread: <http://openmx.ssri.psu.edu/thread/765>
#'
#' Note: This function is currently not robust across multi-group designs or
#' definition variables. It is designed to provide residual-based fit indices
#' (SRMR, CRMR, SMAR, CMAR, etc.) and less-often reported fit indices where
#' Reviewer 2 wants something other than CFA/TLI/RMSEA.
#'
#'
#' Fit information reported includes:
#'
#' **Model characteristics**:
#' numObs, estimated parameters,
#' observed statistics, observed summary statistics,
#' -2*log(Likelihood), degrees of freedom
#'
#' **Chi-squared test**:
#' Chi, ChiDoF, p (of Chi), ChiPerDoF,
#'
#' **Noncentrality-based indices**:
#' RMSEA, RMSEACI, RMSEANull, RMSEAClose (p value), independenceRMSEA,
#' NCP, NCPCI, F0, F0CI, Mc (aka NCI, MFI)
#'
#' **Comparative fit indices**:
#' TLI (aka NNFI), CFI, IFI, PRATIO, PCFI
#'
#' **Residual-based indices**:
#' RMR, SRMR, SRMR_mplus, CRMR,
#' MAR, SMAR, SMAR_mplus, CMAR
#'
#' **Information-theory criteria** (computed using chi-square or -2LL; df or parameters penalties)
#' AIC, AICc, BIC, SABIC, CAIC, BCC
#' ECVI, ECVICI, MECVI, MECVICI
#'
#' **LISREL and other early fit indices** (we recommend not reporting these)
#' GFI, AGFI, PGFI, GH, NFI, PNFI, RFI
#'
#'
#' Want more? *Open an Issue* at [GitHub](https://github.com/tbates/umx/issues).
#'
#' @param model The \code{\link{mxModel}} for which you want fit indices.
#' @param ... Additional parameters passed to \code{\link{summary.MxModel}}.
#' @return List of fit statistics
#' @export
#' @family Reporting functions
#' @author Brenton M. Wiernik, Athanassios Protopapas, Paolo Ghisletta, Markus Brauer
#' @md
#' @examples
#' \dontrun{
#' library(umx)
#' data(demoOneFactor)
#' latents = c("G")
#' manifests = names(demoOneFactor)
#' m1 = umxRAM("One Factor",
#' data = mxData(cov(demoOneFactor), type = "cov", numObs = 500),
#' umxPath(latents, to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = latents, fixedAt = 1)
#' )
#' umxFitIndices(m1)
#'
#' # And with raw data
#' m2 = umxRAM("m1", data = demoOneFactor,
#' umxPath(latents, to = manifests),
#' umxPath(v.m. = manifests),
#' umxPath(v1m0 = latents)
#' )
#' umxFitIndices(m1, refModels = mxRefModels(m2, run = TRUE))
#' }
#'
umxFitIndices = function(model, ...) {
if (is.list(model)) {
stop("Lists of models not supported in `umxFitIndices()`.\nSpecify one model at a time.", call. = FALSE)
}
if (!umx_is_MxModel(model)) {
stop("`model` must be an MxModel.", call. = FALSE)
}
if (!umx_is_RAM(model)) {
message("Not fully tested against non-RAM models...")
}
modSumm = summary(model, ...)
numObs = modSumm$numObs
nParam = modSumm$estimatedParameters
deviance = modSumm$Minus2LogLikelihood
df = modSumm$degreesOfFreedom
devInd = modSumm$IndependenceLikelihood
devSat = modSumm$SaturatedLikelihood
dfInd = modSumm$independenceDoF
dfSat = modSumm$saturatedDoF
Chi = modSumm$Chi
dfChi = modSumm$ChiDoF
pChi = modSumm$p
ChiPerDf = Chi / dfChi
ChiInd = devInd - devSat
dfChiInd = dfInd - dfSat
if (model$data$type == "raw") {
if (!is.null(model$manifestVars)) {
manifests = model$manifestVars
nManifest = length(manifests)
} else if (!(length(model$expectation$dims) == 1 && is.na(model$expectation$dims))) {
manifests = model$expectation$dims[model$expectation$dims %in% colnames(model$data$observed)]
nManifest = length(manifests)
} else {
# TODO: Extract dimnames from expectation slot objects according to the MxExpectation class
warning(paste("Models with raw data and expectation type:", class(model$expectation),
"only supported if dimnames are specified in the expectation.",
"\nFit indices computed using the number of observed variables will be `NA`."))
manifests = NA_character_
nManifest = NA_integer_
}
} else {
nManifest = NCOL(model$data$observed)
}
q = (nManifest * (nManifest + 1)) / 2
if (model$data$type == "raw" || ! is.na(model$data$means)) {
q_m = q + nManifest
} else {
q_m = NA_integer_
}
RMSEASquared = modSumm$RMSEASquared
ciRMSEA = modSumm$RMSEACI
NCP = RMSEASquared * (numObs * dfChi) # Chi - dfChi
ciNCP = ciRMSEA^2 * (modSumm$numObs * modSumm$ChiDoF)
F0 = NCP / numObs
ciF0 = ciNCP / numObs
NFI = (devInd - deviance) / (devInd - devSat)
RFI = 1 - (Chi / dfChi) / (ChiInd / dfChiInd)
IFI = (ChiInd - Chi) / (ChiInd - dfChi)
CFI = modSumm$CFI
PRATIO = dfChi / dfChiInd
PNFI = PRATIO * NFI
PCFI = PRATIO * CFI
# https://doi.org/10.1080/10705519409539970 gives a different formula for PCFI:
# PCFI = CFI * df / ((nManifest * (nManifest - 1)) / 2)
# But most users of PCFI use the form given here
GH = nManifest / (nManifest + 2 * F0)
Mc = exp(-0.5 * NCP / numObs)
RMSEASquaredInd = (ChiInd/(dfChiInd) - 1)/numObs
RMSEAInd = ifelse(RMSEASquaredInd < 0, 0, sqrt(RMSEASquaredInd))
# Information-based indices
infoCrit = modSumm$informationCriteria
Fvalue = infoCrit["AIC:", "par"] - 2 * nParam
infoCrit = rbind(modSumm$informationCriteria,
`CAIC:` = c(Fvalue - (log(numObs) + 1) * df,
Fvalue + (log(numObs) + 1) * nParam,
NA),
`BCC:` = c(Fvalue - (2 * numObs / (numObs - nManifest - 1)) * df,
Fvalue + (2 * numObs / (numObs - nManifest - 1)) * nParam,
NA)
)
if (Fvalue == deviance) {
infoCritDev = infoCrit
infoCritChi = infoCrit - deviance + Chi
} else {
infoCritChi = infoCrit
infoCritDev = infoCrit - Chi + deviance
}
informationCriteria = list(Chi = infoCritChi,
deviance = infoCritDev)
ECVI = (Chi + 2 * nParam) / numObs
MECVI = (Chi + 2 * nParam * numObs / (numObs - nManifest - 1)) / numObs
ciECVI = (ciNCP + dfChi + 2 * nParam) / numObs
ciMECVI = (ciNCP + dfChi + 2 * nParam * numObs / (numObs - nManifest - 1)) / numObs
modelMatrices = xmuModelMatrices(model, manifests)
IdManifest = diag(nManifest)
# Residual-based indices
if (!inherits(model$fitfunction, "MxFitFunctionML")) {
warning("Bias correction and confidence intervals for residual-based\nfit indices currently only designed to work with maximum likelihood estimation.")
}
residualIndices = with(modelMatrices,
xmuResidualIndices(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans))
# LISREL indices
if (model$data$type == "cor") {
invSigS = solve(modelMatrices$estCor) %*% modelMatrices$obsCor
} else {
invSigS = solve(modelMatrices$estCov) %*% modelMatrices$obsCov
}
GFI = 1 - tr(invSigS - IdManifest)^2 / tr(invSigS)^2
AGFI = 1 - (q / dfChi) * (1 - GFI)
PGFI = GFI * dfChi / q
indices =
list(
numObs = numObs,
estimatedParameters = nParam,
observedStatistics = modSumm$observedStatistics,
observedSummaryStatistics = ifelse(model$data$type == "raw" || !(length(model$data$means) == 1 && is.na(model$data$means)), q_m, q),
Minus2LogLikelihood = deviance,
degreesOfFreedom = df,
Chi = Chi,
ChiDoF = dfChi,
p = pChi,
ChiPerDoF = ChiPerDf,
RMSEA = modSumm$RMSEA,
RMSEACI = ciRMSEA,
RMSEANull = modSumm$RMSEANull,
RMSEAClose = modSumm$RMSEAClose,
independenceRMSEA = RMSEAInd,
NCP = NCP, NCPCI = ciNCP,
F0 = F0, F0CI = ciF0,
Mc = Mc,
residualIndices = residualIndices,
GFI = GFI, AGFI = AGFI, PGFI = PGFI, GH = GH,
NFI = NFI, RFI = RFI, IFI = IFI,
TLI = modSumm$TLI, CFI = CFI,
PRATIO = PRATIO, PNFI = PNFI, PCFI = PCFI,
informationCriteria = informationCriteria,
ECVI = ECVI, ECVICI = ciECVI, MECVI = MECVI, MECVICI = ciMECVI
)
class(indices) = c("umxFitIndices", "list")
attr(indices, "name") <- attr(model, "name")
attr(indices, "call") <- deparse(sys.call())
return(indices)
}
xmuModelMatrices = function(model, manifests) {
if (model$data$type == "raw"){
if (any(is.na(model$data$observed))) {
obsFiml = psych::corFiml(model$data$observed[,manifests], covar = TRUE)
obsCov = obsFiml$cov
obsCor = obsFiml$cor
obsMeans = obsFiml$mean
} else {
obsCov = cov(model$data$observed[,manifests])
obsCor = cov2cor(obsCov)
obsMeans = colMeans(model$data$observed)
}
} else if (model$data$type == "cor") {
obsCov = NA_real_
obsCor = model$data$observed
obsMeans = model$data$means
} else {
obsCov = model$data$observed
obsCor = cov2cor(obsCov)
obsMeans = model$data$means
}
estCov = mxGetExpected(model, "covariance")
estCor = cov2cor(estCov)
if(!(length(obsMeans == 1) && is.na(obsMeans))) {
estMeans = mxGetExpected(model, "means")
} else {
estMeans = NA_real_
}
return(list(
obsCov = obsCov,
obsCor = obsCor,
obsMeans = obsMeans,
estCov = estCov,
estCor = estCor,
estMeans = estMeans
))
}
xmuResiduals = function(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans,
which = c("raw", "cor.bentler", "cor.bollen", "cor.mplus")) {
which = match.arg(which, c("raw", "normalized", "standardized",
"standardized.mplus", "cor.bentler",
"cor.bollen", "cor.mplus"), several.ok = TRUE)
# Raw
raw = list(
Cov = obsCov - estCov,
Means = obsMeans - estMeans
)
# TODO: Compute normalized and standardized residuals
# # Normalized
# normalized = c(
# Cov = obsCov - estCov
# Means = obsMeans - estMeans
# )
#
# # Standardized
# standardized = c(
# Cov = obsCov - estCov
# Means = obsMeans - estMeans
# )
#
# # Standardized - MPLus-style
# standardized.mplus = c(
# Cov = obsCov - estCov
# Means = obsMeans - estMeans
# )
# Correlation - EQS/Benter-style
if (length(obsCov) == 1 && is.na(obsCov)) {
cor.bentler = list(
Cov = matrix(NA_real_, nrow = nrow(obsCor), ncol = ncol(obsCov)),
Means = rep(NA_real_, length(obsMeans))
)
} else {
cor.bentler = list(
Cov = (obsCov - estCov) / sqrt(diag(obsCov) %*% t(diag(obsCov))),
Means = ifelse(length(obsMeans) == 1 && is.na(obsMeans), NA_real_, (obsMeans - estMeans) / sqrt(diag(obsCov)))
)
}
# Correlation - Bollen-style
cor.bollen = list(
Cov = obsCor - estCor,
Means = ifelse(length(obsMeans) == 1 && is.na(obsMeans), NA_real_, obsMeans / sqrt(diag(obsCov)) - estMeans / sqrt(diag(estCov)))
)
# Correlation - MPlus-style (for Information = Observed; http://www.statmodel.com/download/SRMR.pdf)
if (length(obsCov) == 1 && is.na(obsCov)) {
cor.mplus = cor.bentler
} else {
cor.mplus = cor.bollen
diag(cor.mplus$Cov) = diag(cor.bentler$Cov)
}
return(list(
raw = raw,
# normalized = normalized,
# standardized = standardized,
# standardized.mplus = standardized.mplus,
cor.bentler = cor.bentler,
cor.bollen = cor.bollen,
cor.mplus = cor.mplus
)[which])
}
xmuResidualIndices = function(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans) {
# TODO: Add unbiased stdandardized estimators
# TODO: Add confidence intervals
resMat = xmuResiduals(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans,
which = c("raw", "cor.bentler", "cor.bollen", "cor.mplus"))
# SRMR/SMAR = residuals standardized using the EQS/Bentler approach
# SRMR_mplus/SMAR_mplus = residuals standardized using the MPlus approach (Information=Observed; http://www.statmodel.com/download/SRMR.pdf)
# CRMR/CMAR = residuals standardized using the Bollen approach (residual correlations)
if(length(obsMeans) == 1 && is.na(obsMeans)) {
if (!(length(obsCov) == 1 && is.na(obsCov))) {
RMR = sqrt(mean( c(vech(resMat$raw$Cov))^2))
SRMR = sqrt(mean( c(vech(resMat$cor.bentler$Cov))^2))
SRMR_mplus = sqrt(mean( c(vech(resMat$cor.mplus$Cov))^2))
MAR = mean(abs(c(vech(resMat$raw$Cov))))
SMAR = mean(abs(c(vech(resMat$cor.bentler$Cov))))
SMAR_mplus = mean(abs(c(vech(resMat$cor.mplus$Cov))))
} else {
RMR = SRMR = SRMR_mplus = MAR = SMAR = SMAR_mplus = NA
}
CRMR = sqrt(mean( c(vech(resMat$cor.bollen$Cov))^2))
CMAR = mean(abs(c(vech(resMat$cor.bollen$Cov))))
} else {
if (!(length(obsCov) == 1 && is.na(obsCov))) {
RMR = sqrt(mean( c(vech(resMat$raw$Cov), resMat$raw$Means)^2))
SRMR = sqrt(mean( c(vech(resMat$cor.bentler$Cov), resMat$cor.bentler$Means)^2))
SRMR_mplus = sqrt(mean( c(vech(resMat$cor.mplus$Cov), resMat$cor.mplus$Means)^2))
MAR = mean(abs(c(vech(resMat$raw$Cov), resMat$raw$Means)))
SMAR = mean(abs(c(vech(resMat$cor.bentler$Cov), resMat$cor.bentler$Means)))
SMAR_mplus = mean(abs(c(vech(resMat$cor.mplus$Cov), resMat$cor.mplus$Means)))
} else {
RMR = SRMR = SRMR_mplus = MAR = SMAR = SMAR_mplus = NA
}
CRMR = sqrt(mean( c(vech(resMat$cor.bollen$Cov), resMat$cor.bollen$Means)^2))
CMAR = mean(abs(c(vech(resMat$cor.bollen$Cov), resMat$cor.bollen$Means)))
}
# Return valus structure in preparation for confidence intervals and SEs
return(
list(
estimate = c(RMR = RMR, SRMR = SRMR, SRMR_mplus = SRMR_mplus, CRMR = CRMR,
MAR = MAR, SMAR = SMAR, SMAR_mplus = SMAR_mplus, CMAR = CMAR)
)
)
}
print.umxFitIndices = function(x, digits = max(1L, getOption("digits") - 3L), ...) {
if (!inherits(x, "umxFitIndices")) {
stop(gettextf("'x' must inherit from class %s", dQuote("umxFitIndices")),
domain = NA)
}
cat('\nFit indices for', attr(x, "name"), "\n\n")
cat('Model characteristics\n')
cat('=====================\n')
cat('Number of observations (sample size): ', round(x$numObs, digits), '\n')
cat('Number of estimated parameters: ', round(x$estimatedParameters, digits), '\n')
cat('Number of observed statistics: ', round(x$observedStatistics, digits), '\n')
cat('Number of observed summary statistics: ', round(x$observedSummaryStatistics, digits), '\n')
cat(' (means, variances, covariances)\n')
cat('Deviance (-2 * LogLikelihood): ', round(x$Minus2LogLikelihood, digits), '\n')
cat('\n\n')
if (!is.na(x$Chi)) {
cat('Chi squared test\n')
cat('================\n')
if (x$Chi <= 0) {
cat("chi-square: ", "\u03C7\u00B2 (df=", x$ChiDoF, ") = ",
round(x$Chi, digits),
", p = ", format.pval(x$p, digits, eps = 0), "\n", sep = "")
} else {
cat("chi-square: ", "\u03C7\u00B2 (df=", x$ChiDoF, ") = ",
format(round(x$Chi, max(0, digits - log10(x$Chi)))),
", p = ", format.pval(x$p, digits, eps = 0), "\n", sep = "")
}
cat('\n')
cat('chi-square per df (\u03C7\u00B2 / df):', round(x$ChiPerDoF, digits), '\n')
cat('\n\n')
cat('Noncentrality-based indices\n')
cat('===========================\n')
if (length(x$RMSEASquared) == 1 && !is.na(x$RMSEASquared) && x$RMSEASquared < 0) {
cat('** (Non-centrality parameter is negative)')
}
cat('RMSEA (root mean squared error of approximation): ', round(x$RMSEA, digits),
' [95% CI (', round(x$RMSEACI[1], digits), ', ',
round(x$RMSEACI[2], digits), ')]', '\n', sep = "")
cat(' Prob(RMSEA <= ', round(x$RMSEANull, digits), "): ",
format.pval(x$RMSEAClose, digits, eps = 0), "\n", sep = "")
cat('\n')
cat('RMSEA of independence model: ',
rep(' ', nchar(round(x$NCP, digits)) - nchar(round(x$independenceRMSEA, digits))),
round(x$independenceRMSEA, digits), '\n', sep = "")
cat('Noncentrality parameter (NCP or d): ', round(x$NCP, digits),
' [95% CI (', round(x$NCPCI[1], digits), ', ',
round(x$NCPCI[2], digits), ')]', '\n', sep = "")
cat('Rescaled NCP (F0 or t): ',
rep(' ', nchar(round(x$NCP, digits)) - nchar(round(x$F0, digits))),
round(x$F0, digits),
' [95% CI (', round(x$F0CI[1], digits), ', ',
round(x$F0CI[2], digits), ')]', '\n', sep = "")
cat('Mc (McDonald centrality index): ',
rep(' ', nchar(round(x$NCP, digits)) - nchar(round(x$Mc, digits))),
round(x$Mc, digits),
' [Also called MFI (McDonald fit index)\n',
rep(' ', 46), 'or NCI (Noncentrality index)]\n', sep = "")
cat('\n\n')
cat('Incremental fit indices\n')
cat('=======================\n')
cat('CFI (Comparative fit index):', round(x$CFI, digits), ' [Also called RNI (Relative noncentrality index)]\n')
cat('TLI (Tucker-Lewis index): ', round(x$TLI, digits), ' [Also called NNFI (Non-normed fit index)]\n')
cat('IFI (Incremental fit index):', round(x$IFI, digits), ' [Also called BL89 (Bollen, 1989, fit index)]\n')
cat('\n')
cat('Parsimony ratio (df / df_indepedence):', round(x$PRATIO, digits), '\n')
cat('PCFI (Parsimonious CFI): ', round(x$PCFI, digits), '\n')
if (!is.na(x$independenceRMSEA) && x$independenceRMSEA <= .158) {
message(paste0('\nNote: Independence (Null) model has RMSEA = ',
round(x$independenceRMSEA, digits),
'.\nIf the model shows good fit (RMSEA <= .05), TLI has a maximum value <= .90.',
'\nInterpret incremental fit indices (TLI, CFI, etc.) with caution.\n'))
}
cat('\n\n')
}
squaredResisualIndices = data.frame(
Estimate = round(x$residualIndices$estimate[c('RMR', 'SRMR', 'SRMR_mplus', 'CRMR')], digits)
)
rownames(squaredResisualIndices) =
c('RMR (Root mean squared residual):',
'SRMR (Standardized root mean squared residual):',
' (MPlus method):',
'CRMR (Correlation root mean squared residual):'
)
absoluteResisualIndices = data.frame(
Estimate = round(x$residualIndices$estimate[c('MAR', 'SMAR', 'SMAR_mplus', 'CMAR')], digits)
)
rownames(absoluteResisualIndices) =
c('MAR (Mean absolute residual):',
'SMAR (Standardized mean absolute residual):',
' (MPlus method): ',
'CMAR (Correlation mean absolute residual):'
)
cat('Residuals-based indices\n')
cat('=======================\n')
if (is.na(x$residualIndices$estimate['SRMR'])) {
cat('Observed values are correlations',
' Use CRMR instead of RMR or SRMR.',
' Use CMAR instead of MAR or SMAR.',
'\n')
print(rbind(squaredResisualIndices[4,], ` ` = ""))
cat('\n')
cat('CRMR, CMAR:\n Residuals standaridzed using the Bollen method.\n')
} else {
print(rbind(squaredResisualIndices, ` ` = "", absoluteResisualIndices))
cat('\n')
cat(
' SRMR, SMAR: Residuals standardized using the Bentler method (used by lavaan, sem, EQS, AMOS).\n',
' (MPlus): Residuals standardized using the MPlus "observed information" method.\n',
' CRMR, CMAR: Residuals standaridzed using the Bollen method.\n',
sep = ""
)
}
cat('\n\n')
colnames(x$informationCriteria$Chi) =
colnames(x$informationCriteria$deviance) =
c(" | df Penalty", " | Parameters Penalty", " | Sample-Size Adjusted")
rownames(x$informationCriteria$Chi) =
rownames(x$informationCriteria$deviance) =
c('AIC (Akaike information criterion):',
'BIC (Bayesian information criterion):',
'CAIC (Consistent AIC):',
'BCC (Browne-Cudeck criterion):'
)
if (!is.na(x$Chi)) {
cat('Information-based fit indices (computed using \u03C7\u00B2)\n')
cat('=================================================\n')
print(round(x$informationCriteria$Chi, digits))
cat("\n")
cat('ECVI (Expected cross-validation index): ', round(x$ECVI, digits),
' [95% CI (', round(x$ECVICI[1], digits), ', ',
round(x$ECVICI[2], digits), ')]', '\n', sep = "")
cat('MECVI (Modified ECVI): ', round(x$ECVI, digits),
' [95% CI (', round(x$MECVICI[1], digits), ', ',
round(x$MECVICI[2], digits), ')]', '\n', sep = "")
cat('\n\n')
}
cat('Information-based fit indices (computed using -2lnL)\n')
cat('====================================================\n')
print(round(x$informationCriteria$deviance, digits))
cat("\n\n")
cat('LISREL and other early fit indices\n')
cat('==================================\n')
cat('** These indices are strongly discouraged due to numerous problems.',
'\n We recommend that you do not report them.\n')
cat('GFI (Goodness of fit index):', round(x$GFI, digits), '\n')
cat('AGFI (Adjusted GFI): ', round(x$AGFI, digits), '\n')
cat('PGFI (Parsimonious GFI): ', round(x$PGFI, digits), '\n')
if (!is.na(x$Chi)) {
cat('NFI (Normed fit index): ', round(x$NFI, digits), '\n')
cat('PNFI (Parsimonious NFI): ', round(x$PNFI, digits), ' [Also called PFI (Parsimonious fit index)]\n')
cat('RFI (Relative fit index): ', round(x$RFI, digits), '\n')
cat('GH (Gamma hat): ', round(x$GH, digits), ' [Estimated population GFI]\n')
}
objectName = regmatches(attr(x, "call"),
regexec("umxFitIndices\\((?:model *= *)?([^,)]+)",
attr(x, "call"))
)[[1]][[2]]
if (is.na(x$Chi)) {
cat('\n\nFor additional fit indices, specify reference models using:',
'\n umxFitIndices(', objectName,
', refModels = mxRefModels(', objectName, ', run = TRUE))', sep="")
}
invisible(x)
}
tbates/umx documentation built on Sept. 19, 2024, 1:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.umxFitIndices xmuResidualIndices xmuResiduals xmuModelMatrices umxFitIndices
Documented in umxFitIndices
#' Get additional fit-indices for a model with umxFitIndices
#'
#' Computes a variety of fit indices. Originated in this thread: <http://openmx.ssri.psu.edu/thread/765>
#'
#' Note: This function is currently not robust across multi-group designs or
#' definition variables. It is designed to provide residual-based fit indices
#' (SRMR, CRMR, SMAR, CMAR, etc.) and less-often reported fit indices where
#' Reviewer 2 wants something other than CFA/TLI/RMSEA.
#'
#'
#' Fit information reported includes:
#'
#' **Model characteristics**:
#' numObs, estimated parameters,
#' observed statistics, observed summary statistics,
#' -2*log(Likelihood), degrees of freedom
#'
#' **Chi-squared test**:
#' Chi, ChiDoF, p (of Chi), ChiPerDoF,
#'
#' **Noncentrality-based indices**:
#' RMSEA, RMSEACI, RMSEANull, RMSEAClose (p value), independenceRMSEA,
#' NCP, NCPCI, F0, F0CI, Mc (aka NCI, MFI)
#'
#' **Comparative fit indices**:
#' TLI (aka NNFI), CFI, IFI, PRATIO, PCFI
#'
#' **Residual-based indices**:
#' RMR, SRMR, SRMR_mplus, CRMR,
#' MAR, SMAR, SMAR_mplus, CMAR
#'
#' **Information-theory criteria** (computed using chi-square or -2LL; df or parameters penalties)
#' AIC, AICc, BIC, SABIC, CAIC, BCC
#' ECVI, ECVICI, MECVI, MECVICI
#'
#' **LISREL and other early fit indices** (we recommend not reporting these)
#' GFI, AGFI, PGFI, GH, NFI, PNFI, RFI
#'
#'
#' Want more? *Open an Issue* at [GitHub](https://github.com/tbates/umx/issues).
#'
#' @param model The \code{\link{mxModel}} for which you want fit indices.
#' @param ... Additional parameters passed to \code{\link{summary.MxModel}}.
#' @return List of fit statistics
#' @export
#' @family Reporting functions
#' @author Brenton M. Wiernik, Athanassios Protopapas, Paolo Ghisletta, Markus Brauer
#' @md
#' @examples
#' \dontrun{
#' library(umx)
#' data(demoOneFactor)
#' latents = c("G")
#' manifests = names(demoOneFactor)
#' m1 = umxRAM("One Factor",
#' data = mxData(cov(demoOneFactor), type = "cov", numObs = 500),
#' umxPath(latents, to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = latents, fixedAt = 1)
#' )
#' umxFitIndices(m1)
#'
#' # And with raw data
#' m2 = umxRAM("m1", data = demoOneFactor,
#' umxPath(latents, to = manifests),
#' umxPath(v.m. = manifests),
#' umxPath(v1m0 = latents)
#' )
#' umxFitIndices(m1, refModels = mxRefModels(m2, run = TRUE))
#' }
#'
umxFitIndices = function(model, ...) {
if (is.list(model)) {
stop("Lists of models not supported in `umxFitIndices()`.\nSpecify one model at a time.", call. = FALSE)
}
if (!umx_is_MxModel(model)) {
stop("`model` must be an MxModel.", call. = FALSE)
}
if (!umx_is_RAM(model)) {
message("Not fully tested against non-RAM models...")
}
modSumm = summary(model, ...)
numObs = modSumm$numObs
nParam = modSumm$estimatedParameters
deviance = modSumm$Minus2LogLikelihood
df = modSumm$degreesOfFreedom
devInd = modSumm$IndependenceLikelihood
devSat = modSumm$SaturatedLikelihood
dfInd = modSumm$independenceDoF
dfSat = modSumm$saturatedDoF
Chi = modSumm$Chi
dfChi = modSumm$ChiDoF
pChi = modSumm$p
ChiPerDf = Chi / dfChi
ChiInd = devInd - devSat
dfChiInd = dfInd - dfSat
if (model$data$type == "raw") {
if (!is.null(model$manifestVars)) {
manifests = model$manifestVars
nManifest = length(manifests)
} else if (!(length(model$expectation$dims) == 1 && is.na(model$expectation$dims))) {
manifests = model$expectation$dims[model$expectation$dims %in% colnames(model$data$observed)]
nManifest = length(manifests)
} else {
# TODO: Extract dimnames from expectation slot objects according to the MxExpectation class
warning(paste("Models with raw data and expectation type:", class(model$expectation),
"only supported if dimnames are specified in the expectation.",
"\nFit indices computed using the number of observed variables will be `NA`."))
manifests = NA_character_
nManifest = NA_integer_
}
} else {
nManifest = NCOL(model$data$observed)
}
q = (nManifest * (nManifest + 1)) / 2
if (model$data$type == "raw" || ! is.na(model$data$means)) {
q_m = q + nManifest
} else {
q_m = NA_integer_
}
RMSEASquared = modSumm$RMSEASquared
ciRMSEA = modSumm$RMSEACI
NCP = RMSEASquared * (numObs * dfChi) # Chi - dfChi
ciNCP = ciRMSEA^2 * (modSumm$numObs * modSumm$ChiDoF)
F0 = NCP / numObs
ciF0 = ciNCP / numObs
NFI = (devInd - deviance) / (devInd - devSat)
RFI = 1 - (Chi / dfChi) / (ChiInd / dfChiInd)
IFI = (ChiInd - Chi) / (ChiInd - dfChi)
CFI = modSumm$CFI
PRATIO = dfChi / dfChiInd
PNFI = PRATIO * NFI
PCFI = PRATIO * CFI
# https://doi.org/10.1080/10705519409539970 gives a different formula for PCFI:
# PCFI = CFI * df / ((nManifest * (nManifest - 1)) / 2)
# But most users of PCFI use the form given here
GH = nManifest / (nManifest + 2 * F0)
Mc = exp(-0.5 * NCP / numObs)
RMSEASquaredInd = (ChiInd/(dfChiInd) - 1)/numObs
RMSEAInd = ifelse(RMSEASquaredInd < 0, 0, sqrt(RMSEASquaredInd))
# Information-based indices
infoCrit = modSumm$informationCriteria
Fvalue = infoCrit["AIC:", "par"] - 2 * nParam
infoCrit = rbind(modSumm$informationCriteria,
`CAIC:` = c(Fvalue - (log(numObs) + 1) * df,
Fvalue + (log(numObs) + 1) * nParam,
NA),
`BCC:` = c(Fvalue - (2 * numObs / (numObs - nManifest - 1)) * df,
Fvalue + (2 * numObs / (numObs - nManifest - 1)) * nParam,
NA)
)
if (Fvalue == deviance) {
infoCritDev = infoCrit
infoCritChi = infoCrit - deviance + Chi
} else {
infoCritChi = infoCrit
infoCritDev = infoCrit - Chi + deviance
}
informationCriteria = list(Chi = infoCritChi,
deviance = infoCritDev)
ECVI = (Chi + 2 * nParam) / numObs
MECVI = (Chi + 2 * nParam * numObs / (numObs - nManifest - 1)) / numObs
ciECVI = (ciNCP + dfChi + 2 * nParam) / numObs
ciMECVI = (ciNCP + dfChi + 2 * nParam * numObs / (numObs - nManifest - 1)) / numObs
modelMatrices = xmuModelMatrices(model, manifests)
IdManifest = diag(nManifest)
# Residual-based indices
if (!inherits(model$fitfunction, "MxFitFunctionML")) {
warning("Bias correction and confidence intervals for residual-based\nfit indices currently only designed to work with maximum likelihood estimation.")
}
residualIndices = with(modelMatrices,
xmuResidualIndices(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans))
# LISREL indices
if (model$data$type == "cor") {
invSigS = solve(modelMatrices$estCor) %*% modelMatrices$obsCor
} else {
invSigS = solve(modelMatrices$estCov) %*% modelMatrices$obsCov
}
GFI = 1 - tr(invSigS - IdManifest)^2 / tr(invSigS)^2
AGFI = 1 - (q / dfChi) * (1 - GFI)
PGFI = GFI * dfChi / q
indices =
list(
numObs = numObs,
estimatedParameters = nParam,
observedStatistics = modSumm$observedStatistics,
observedSummaryStatistics = ifelse(model$data$type == "raw" || !(length(model$data$means) == 1 && is.na(model$data$means)), q_m, q),
Minus2LogLikelihood = deviance,
degreesOfFreedom = df,
Chi = Chi,
ChiDoF = dfChi,
p = pChi,
ChiPerDoF = ChiPerDf,
RMSEA = modSumm$RMSEA,
RMSEACI = ciRMSEA,
RMSEANull = modSumm$RMSEANull,
RMSEAClose = modSumm$RMSEAClose,
independenceRMSEA = RMSEAInd,
NCP = NCP, NCPCI = ciNCP,
F0 = F0, F0CI = ciF0,
Mc = Mc,
residualIndices = residualIndices,
GFI = GFI, AGFI = AGFI, PGFI = PGFI, GH = GH,
NFI = NFI, RFI = RFI, IFI = IFI,
TLI = modSumm$TLI, CFI = CFI,
PRATIO = PRATIO, PNFI = PNFI, PCFI = PCFI,
informationCriteria = informationCriteria,
ECVI = ECVI, ECVICI = ciECVI, MECVI = MECVI, MECVICI = ciMECVI
)
class(indices) = c("umxFitIndices", "list")
attr(indices, "name") <- attr(model, "name")
attr(indices, "call") <- deparse(sys.call())
return(indices)
}
xmuModelMatrices = function(model, manifests) {
if (model$data$type == "raw"){
if (any(is.na(model$data$observed))) {
obsFiml = psych::corFiml(model$data$observed[,manifests], covar = TRUE)
obsCov = obsFiml$cov
obsCor = obsFiml$cor
obsMeans = obsFiml$mean
} else {
obsCov = cov(model$data$observed[,manifests])
obsCor = cov2cor(obsCov)
obsMeans = colMeans(model$data$observed)
}
} else if (model$data$type == "cor") {
obsCov = NA_real_
obsCor = model$data$observed
obsMeans = model$data$means
} else {
obsCov = model$data$observed
obsCor = cov2cor(obsCov)
obsMeans = model$data$means
}
estCov = mxGetExpected(model, "covariance")
estCor = cov2cor(estCov)
if(!(length(obsMeans == 1) && is.na(obsMeans))) {
estMeans = mxGetExpected(model, "means")
} else {
estMeans = NA_real_
}
return(list(
obsCov = obsCov,
obsCor = obsCor,
obsMeans = obsMeans,
estCov = estCov,
estCor = estCor,
estMeans = estMeans
))
}
xmuResiduals = function(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans,
which = c("raw", "cor.bentler", "cor.bollen", "cor.mplus")) {
which = match.arg(which, c("raw", "normalized", "standardized",
"standardized.mplus", "cor.bentler",
"cor.bollen", "cor.mplus"), several.ok = TRUE)
# Raw
raw = list(
Cov = obsCov - estCov,
Means = obsMeans - estMeans
)
# TODO: Compute normalized and standardized residuals
# # Normalized
# normalized = c(
# Cov = obsCov - estCov
# Means = obsMeans - estMeans
# )
#
# # Standardized
# standardized = c(
# Cov = obsCov - estCov
# Means = obsMeans - estMeans
# )
#
# # Standardized - MPLus-style
# standardized.mplus = c(
# Cov = obsCov - estCov
# Means = obsMeans - estMeans
# )
# Correlation - EQS/Benter-style
if (length(obsCov) == 1 && is.na(obsCov)) {
cor.bentler = list(
Cov = matrix(NA_real_, nrow = nrow(obsCor), ncol = ncol(obsCov)),
Means = rep(NA_real_, length(obsMeans))
)
} else {
cor.bentler = list(
Cov = (obsCov - estCov) / sqrt(diag(obsCov) %*% t(diag(obsCov))),
Means = ifelse(length(obsMeans) == 1 && is.na(obsMeans), NA_real_, (obsMeans - estMeans) / sqrt(diag(obsCov)))
)
}
# Correlation - Bollen-style
cor.bollen = list(
Cov = obsCor - estCor,
Means = ifelse(length(obsMeans) == 1 && is.na(obsMeans), NA_real_, obsMeans / sqrt(diag(obsCov)) - estMeans / sqrt(diag(estCov)))
)
# Correlation - MPlus-style (for Information = Observed; http://www.statmodel.com/download/SRMR.pdf)
if (length(obsCov) == 1 && is.na(obsCov)) {
cor.mplus = cor.bentler
} else {
cor.mplus = cor.bollen
diag(cor.mplus$Cov) = diag(cor.bentler$Cov)
}
return(list(
raw = raw,
# normalized = normalized,
# standardized = standardized,
# standardized.mplus = standardized.mplus,
cor.bentler = cor.bentler,
cor.bollen = cor.bollen,
cor.mplus = cor.mplus
)[which])
}
xmuResidualIndices = function(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans) {
# TODO: Add unbiased stdandardized estimators
# TODO: Add confidence intervals
resMat = xmuResiduals(obsCov, obsCor, obsMeans,
estCov, estCor, estMeans,
which = c("raw", "cor.bentler", "cor.bollen", "cor.mplus"))
# SRMR/SMAR = residuals standardized using the EQS/Bentler approach
# SRMR_mplus/SMAR_mplus = residuals standardized using the MPlus approach (Information=Observed; http://www.statmodel.com/download/SRMR.pdf)
# CRMR/CMAR = residuals standardized using the Bollen approach (residual correlations)
if(length(obsMeans) == 1 && is.na(obsMeans)) {
if (!(length(obsCov) == 1 && is.na(obsCov))) {
RMR = sqrt(mean( c(vech(resMat$raw$Cov))^2))
SRMR = sqrt(mean( c(vech(resMat$cor.bentler$Cov))^2))
SRMR_mplus = sqrt(mean( c(vech(resMat$cor.mplus$Cov))^2))
MAR = mean(abs(c(vech(resMat$raw$Cov))))
SMAR = mean(abs(c(vech(resMat$cor.bentler$Cov))))
SMAR_mplus = mean(abs(c(vech(resMat$cor.mplus$Cov))))
} else {
RMR = SRMR = SRMR_mplus = MAR = SMAR = SMAR_mplus = NA
}
CRMR = sqrt(mean( c(vech(resMat$cor.bollen$Cov))^2))
CMAR = mean(abs(c(vech(resMat$cor.bollen$Cov))))
} else {
if (!(length(obsCov) == 1 && is.na(obsCov))) {
RMR = sqrt(mean( c(vech(resMat$raw$Cov), resMat$raw$Means)^2))
SRMR = sqrt(mean( c(vech(resMat$cor.bentler$Cov), resMat$cor.bentler$Means)^2))
SRMR_mplus = sqrt(mean( c(vech(resMat$cor.mplus$Cov), resMat$cor.mplus$Means)^2))
MAR = mean(abs(c(vech(resMat$raw$Cov), resMat$raw$Means)))
SMAR = mean(abs(c(vech(resMat$cor.bentler$Cov), resMat$cor.bentler$Means)))
SMAR_mplus = mean(abs(c(vech(resMat$cor.mplus$Cov), resMat$cor.mplus$Means)))
} else {
RMR = SRMR = SRMR_mplus = MAR = SMAR = SMAR_mplus = NA
}
CRMR = sqrt(mean( c(vech(resMat$cor.bollen$Cov), resMat$cor.bollen$Means)^2))
CMAR = mean(abs(c(vech(resMat$cor.bollen$Cov), resMat$cor.bollen$Means)))
}
# Return valus structure in preparation for confidence intervals and SEs
return(
list(
estimate = c(RMR = RMR, SRMR = SRMR, SRMR_mplus = SRMR_mplus, CRMR = CRMR,
MAR = MAR, SMAR = SMAR, SMAR_mplus = SMAR_mplus, CMAR = CMAR)
)
)
}
print.umxFitIndices = function(x, digits = max(1L, getOption("digits") - 3L), ...) {
if (!inherits(x, "umxFitIndices")) {
stop(gettextf("'x' must inherit from class %s", dQuote("umxFitIndices")),
domain = NA)
}
cat('\nFit indices for', attr(x, "name"), "\n\n")
cat('Model characteristics\n')
cat('=====================\n')
cat('Number of observations (sample size): ', round(x$numObs, digits), '\n')
cat('Number of estimated parameters: ', round(x$estimatedParameters, digits), '\n')
cat('Number of observed statistics: ', round(x$observedStatistics, digits), '\n')
cat('Number of observed summary statistics: ', round(x$observedSummaryStatistics, digits), '\n')
cat(' (means, variances, covariances)\n')
cat('Deviance (-2 * LogLikelihood): ', round(x$Minus2LogLikelihood, digits), '\n')
cat('\n\n')
if (!is.na(x$Chi)) {
cat('Chi squared test\n')
cat('================\n')
if (x$Chi <= 0) {
cat("chi-square: ", "\u03C7\u00B2 (df=", x$ChiDoF, ") = ",
round(x$Chi, digits),
", p = ", format.pval(x$p, digits, eps = 0), "\n", sep = "")
} else {
cat("chi-square: ", "\u03C7\u00B2 (df=", x$ChiDoF, ") = ",
format(round(x$Chi, max(0, digits - log10(x$Chi)))),
", p = ", format.pval(x$p, digits, eps = 0), "\n", sep = "")
}
cat('\n')
cat('chi-square per df (\u03C7\u00B2 / df):', round(x$ChiPerDoF, digits), '\n')
cat('\n\n')
cat('Noncentrality-based indices\n')
cat('===========================\n')
if (length(x$RMSEASquared) == 1 && !is.na(x$RMSEASquared) && x$RMSEASquared < 0) {
cat('** (Non-centrality parameter is negative)')
}
cat('RMSEA (root mean squared error of approximation): ', round(x$RMSEA, digits),
' [95% CI (', round(x$RMSEACI[1], digits), ', ',
round(x$RMSEACI[2], digits), ')]', '\n', sep = "")
cat(' Prob(RMSEA <= ', round(x$RMSEANull, digits), "): ",
format.pval(x$RMSEAClose, digits, eps = 0), "\n", sep = "")
cat('\n')
cat('RMSEA of independence model: ',
rep(' ', nchar(round(x$NCP, digits)) - nchar(round(x$independenceRMSEA, digits))),
round(x$independenceRMSEA, digits), '\n', sep = "")
cat('Noncentrality parameter (NCP or d): ', round(x$NCP, digits),
' [95% CI (', round(x$NCPCI[1], digits), ', ',
round(x$NCPCI[2], digits), ')]', '\n', sep = "")
cat('Rescaled NCP (F0 or t): ',
rep(' ', nchar(round(x$NCP, digits)) - nchar(round(x$F0, digits))),
round(x$F0, digits),
' [95% CI (', round(x$F0CI[1], digits), ', ',
round(x$F0CI[2], digits), ')]', '\n', sep = "")
cat('Mc (McDonald centrality index): ',
rep(' ', nchar(round(x$NCP, digits)) - nchar(round(x$Mc, digits))),
round(x$Mc, digits),
' [Also called MFI (McDonald fit index)\n',
rep(' ', 46), 'or NCI (Noncentrality index)]\n', sep = "")
cat('\n\n')
cat('Incremental fit indices\n')
cat('=======================\n')
cat('CFI (Comparative fit index):', round(x$CFI, digits), ' [Also called RNI (Relative noncentrality index)]\n')
cat('TLI (Tucker-Lewis index): ', round(x$TLI, digits), ' [Also called NNFI (Non-normed fit index)]\n')
cat('IFI (Incremental fit index):', round(x$IFI, digits), ' [Also called BL89 (Bollen, 1989, fit index)]\n')
cat('\n')
cat('Parsimony ratio (df / df_indepedence):', round(x$PRATIO, digits), '\n')
cat('PCFI (Parsimonious CFI): ', round(x$PCFI, digits), '\n')
if (!is.na(x$independenceRMSEA) && x$independenceRMSEA <= .158) {
message(paste0('\nNote: Independence (Null) model has RMSEA = ',
round(x$independenceRMSEA, digits),
'.\nIf the model shows good fit (RMSEA <= .05), TLI has a maximum value <= .90.',
'\nInterpret incremental fit indices (TLI, CFI, etc.) with caution.\n'))
}
cat('\n\n')
}
squaredResisualIndices = data.frame(
Estimate = round(x$residualIndices$estimate[c('RMR', 'SRMR', 'SRMR_mplus', 'CRMR')], digits)
)
rownames(squaredResisualIndices) =
c('RMR (Root mean squared residual):',
'SRMR (Standardized root mean squared residual):',
' (MPlus method):',
'CRMR (Correlation root mean squared residual):'
)
absoluteResisualIndices = data.frame(
Estimate = round(x$residualIndices$estimate[c('MAR', 'SMAR', 'SMAR_mplus', 'CMAR')], digits)
)
rownames(absoluteResisualIndices) =
c('MAR (Mean absolute residual):',
'SMAR (Standardized mean absolute residual):',
' (MPlus method): ',
'CMAR (Correlation mean absolute residual):'
)
cat('Residuals-based indices\n')
cat('=======================\n')
if (is.na(x$residualIndices$estimate['SRMR'])) {
cat('Observed values are correlations',
' Use CRMR instead of RMR or SRMR.',
' Use CMAR instead of MAR or SMAR.',
'\n')
print(rbind(squaredResisualIndices[4,], ` ` = ""))
cat('\n')
cat('CRMR, CMAR:\n Residuals standaridzed using the Bollen method.\n')
} else {
print(rbind(squaredResisualIndices, ` ` = "", absoluteResisualIndices))
cat('\n')
cat(
' SRMR, SMAR: Residuals standardized using the Bentler method (used by lavaan, sem, EQS, AMOS).\n',
' (MPlus): Residuals standardized using the MPlus "observed information" method.\n',
' CRMR, CMAR: Residuals standaridzed using the Bollen method.\n',
sep = ""
)
}
cat('\n\n')
colnames(x$informationCriteria$Chi) =
colnames(x$informationCriteria$deviance) =
c(" | df Penalty", " | Parameters Penalty", " | Sample-Size Adjusted")
rownames(x$informationCriteria$Chi) =
rownames(x$informationCriteria$deviance) =
c('AIC (Akaike information criterion):',
'BIC (Bayesian information criterion):',
'CAIC (Consistent AIC):',
'BCC (Browne-Cudeck criterion):'
)
if (!is.na(x$Chi)) {
cat('Information-based fit indices (computed using \u03C7\u00B2)\n')
cat('=================================================\n')
print(round(x$informationCriteria$Chi, digits))
cat("\n")
cat('ECVI (Expected cross-validation index): ', round(x$ECVI, digits),
' [95% CI (', round(x$ECVICI[1], digits), ', ',
round(x$ECVICI[2], digits), ')]', '\n', sep = "")
cat('MECVI (Modified ECVI): ', round(x$ECVI, digits),
' [95% CI (', round(x$MECVICI[1], digits), ', ',
round(x$MECVICI[2], digits), ')]', '\n', sep = "")
cat('\n\n')
}
cat('Information-based fit indices (computed using -2lnL)\n')
cat('====================================================\n')
print(round(x$informationCriteria$deviance, digits))
cat("\n\n")
cat('LISREL and other early fit indices\n')
cat('==================================\n')
cat('** These indices are strongly discouraged due to numerous problems.',
'\n We recommend that you do not report them.\n')
cat('GFI (Goodness of fit index):', round(x$GFI, digits), '\n')
cat('AGFI (Adjusted GFI): ', round(x$AGFI, digits), '\n')
cat('PGFI (Parsimonious GFI): ', round(x$PGFI, digits), '\n')
if (!is.na(x$Chi)) {
cat('NFI (Normed fit index): ', round(x$NFI, digits), '\n')
cat('PNFI (Parsimonious NFI): ', round(x$PNFI, digits), ' [Also called PFI (Parsimonious fit index)]\n')
cat('RFI (Relative fit index): ', round(x$RFI, digits), '\n')
cat('GH (Gamma hat): ', round(x$GH, digits), ' [Estimated population GFI]\n')
}
objectName = regmatches(attr(x, "call"),
regexec("umxFitIndices\\((?:model *= *)?([^,)]+)",
attr(x, "call"))
)[[1]][[2]]
if (is.na(x$Chi)) {
cat('\n\nFor additional fit indices, specify reference models using:',
'\n umxFitIndices(', objectName,
', refModels = mxRefModels(', objectName, ', run = TRUE))', sep="")
}
invisible(x)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.