Nothing
R/MxSummary.R
In OpenMx: Extended Structural Equation Modelling
Defines functions
omxRMSEA
fitStatistics
catFitStatistics
computeFitStatistics
computeFValue
numberObservations
fitfunctionNumberObservations
observedStatistics
observedStatisticsHelper
calculateConstraintsHelper
calculateConstraints
Documented in
omxRMSEA
#
# Copyright 2007-2021 by the individuals mentioned in the source code history
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
calculateConstraints <- function(model, flatModel) {
constraints <- flatModel@constraints
retval <- c()
if (length(constraints) > 0) {
retval <- sapply(constraints, calculateConstraintsHelper, model)
}
return(retval)
}
calculateConstraintsHelper <- function(constraint, model) {
if (constraint@relation == "==") {
leftHandSide <- constraint@formula[[2]]
value <- eval(substitute(mxEval(x, model, compute=TRUE),
list(x = leftHandSide)))
value <- as.matrix(value)
return(nrow(value) * ncol(value))
} else {
return(0)
}
}
observedStatisticsHelper <- function(model, expectation, datalist, historySet) {
if ('numStats' %in% slotNames(expectation)) {
if (length(expectation@numStats) > 0 && !is.na(expectation@numStats)) {
return(list(expectation@numStats, historySet))
}
}
if (length(expectation@data)==0 || is.na(expectation@data) || !.hasSlot(expectation, 'dims')) {
return(list(0, historySet))
}
if (is.numeric(expectation@data)) {
data <- datalist[[expectation@data + 1]]
} else {
data <- model[[expectation@data]]
}
if (is(data, "MxDataDynamic")) {
# need to revisit TODO
return(list(0, historySet))
}
obsStats <- data@observedStats
if (data@type == 'cov') {
if (data@name %in% historySet) {
return (list(0, historySet))
}
n <- nrow(data@observed)
dof <- n * (n + 1) / 2
if (!single.na(data@means)) {
dof <- dof + length(data@means)
}
historySet <- append(data, historySet)
} else if (data@type == 'cor') {
if (data@name %in% historySet) {
return (list(0, historySet))
}
n <- nrow(data@observed)
dof <- n * (n - 1) / 2
if (!single.na(data@means)) {
dof <- dof + length(data@means)
}
historySet <- append(data, historySet)
} else if (is(expectation, "MxExpectationBA81")) { # refactor TODO
if (!is.na(expectation@weightColumn) || !is.na(data@weight)) {
dof <- nrow(data@observed) - 1
} else {
dof <- nrow(rpf::compressDataFrame(data@observed)) - 1
}
historySet <- append(data, historySet)
} else if (!is.null(obsStats[['cov']])) {
if (data@name %in% historySet) {
return (list(0, historySet))
}
numThresh <- sum(!is.na(obsStats[['thresholds']]))
numMeans <- sum(!is.na(obsStats[['means']]))
numSlope <- sum(!is.na(obsStats[['slope']]))
n <- nrow(obsStats[['cov']])
dof <- numThresh + numSlope + numMeans + n*(n+1)/2 - 2*sum(ncol(obsStats[['thresholds']]))
historySet <- append(data, historySet)
} else {
# Incorporate row frequency and weight information? TODO
dof <- 0
observed <- data@observed
for (i in 1:ncol(observed)) {
colname <- colnames(observed)[[i]]
fullname <- paste(data@name, colname, sep='.')
if ((colname %in% expectation@dims) && !(fullname %in% historySet)) {
dof <- dof + sum(!is.na(observed[,i]))
historySet <- append(fullname, historySet)
}
}
}
return(list(dof, historySet))
}
observedStatistics <- function(model, flatModel, constraintOffset) {
datalist <- flatModel@datasets
labelsData <- imxGenerateLabels(model)
expectations <- convertExpectationFunctions(flatModel, model, labelsData, new("MxDirectedGraph"))
retval <- constraintOffset
if (length(expectations) > 0) {
historySet <- character()
for(i in 1:length(expectations)) {
result <- observedStatisticsHelper(model, expectations[[i]], datalist, historySet)
retval <- retval + result[[1]]
historySet <- result[[2]]
}
}
return(retval)
}
fitfunctionNumberObservations <- function(fitfunction) {
if (("numObsAdjust" %in% slotNames(fitfunction)) && !single.na(fitfunction@numObsAdjust)) {
return(fitfunction@numObsAdjust)
} else {
return(0)
}
}
numberObservations <- function(datalist, fitfunctions) {
datalist <- Filter(function(x) !is(x,"MxDataDynamic"), datalist)
dataObservations <- sapply(datalist, slot, name = "numObs")
fitfunctionObservations <- sapply(fitfunctions, fitfunctionNumberObservations)
return(sum(as.numeric(dataObservations), as.numeric(fitfunctionObservations)))
}
computeFValue <- function(datalist, likelihood, chi) {
if(length(datalist) == 0) return(NA)
datalist <- Filter(function(x) !is(x,"MxDataDynamic"), datalist)
if(all(sapply(datalist, function(x)
{length(x@observedStats) > 0 || is(x,"MxDataLegacyWLS") }))) return(chi)
if(all(sapply(datalist, function(x)
{x@type == 'raw'}))) return(likelihood)
if(all(sapply(datalist, function(x)
{x@type == 'cov'}))) return(chi)
if(all(sapply(datalist, function(x)
{x@type == 'cor'}))) return(chi)
return(NA)
}
computeFitStatistics <- function(likelihood, DoF, chi, chiDoF, numObs,
independence, indDoF, saturated=0, satDoF=0) {
if (!is.na(independence) && !is.na(likelihood) && independence < likelihood) {
warning(paste("Your model may be mis-specified (and fit worse than an independence model),",
"or you may be using the wrong independence model, see ?mxRefModels"))
}
CFI <- (independence - indDoF - likelihood + DoF)/(independence - indDoF - saturated + satDoF)
TLI <- NA
rmseaSquared <- 0
RMSEA <- NA
RMSEACI <- c(lower=NA, upper=NA)
RMSEANull <- 0.05
RMSEAClose <- NA
if (!is.na(chiDoF) && chiDoF > 0) {
TLI <- ((independence-saturated)/(indDoF-satDoF) - (chi)/(DoF-satDoF))/((independence-saturated)/(indDoF-satDoF) - 1)
# Here we use N in the denominator as given in the original
# RMSEA paper. The difference between N and N-1 is negligible
# for sample sizes over 30. RMSEA should not be taken seriously
# such small samples anyway.
rmseaSquared <- (chi / (chiDoF) - 1) / numObs
if (length(rmseaSquared) == 0 || is.na(rmseaSquared) ||
is.nan(rmseaSquared)) {
# || (rmseaSquared < 0)) { # changed so 'rmseaSquared < 0' yields zero with comment
RMSEA <- NA
} else {
RMSEA <- ifelse(rmseaSquared < 0, 0.0, sqrt(rmseaSquared))
ci <- try(rmseaConfidenceIntervalHelper(chi, chiDoF, numObs, .025, .975))
if (!inherits(ci, "try-error")) RMSEACI <- ci
RMSEAClose <- rmseaPCloseHelper(chi, chiDoF, numObs, null=RMSEANull)
}
}
list(CFI=CFI, TLI=TLI, RMSEA=RMSEA, RMSEASquared=rmseaSquared, RMSEACI=RMSEACI, RMSEAClose=RMSEAClose, RMSEANull=RMSEANull)
}
catFitStatistics <- function(x) {
cat("CFI:", x$CFI, '\n')
cat("TLI:", x$TLI, ' (also known as NNFI)', '\n')
if (length(x$RMSEASquared) == 1 && !is.na(x$RMSEASquared) && x$RMSEASquared < 0.0) {
cat("RMSEA: ", x$RMSEA, " *(Non-centrality parameter is negative)", " [95% CI (", x$RMSEACI[1], ", ", x$RMSEACI[2], ")]", '\n', sep="")
} else {
cat("RMSEA: ", x$RMSEA, " [95% CI (", x$RMSEACI[1], ", ", x$RMSEACI[2], ")]", '\n', sep="")
}
cat("Prob(RMSEA <= ", x$RMSEANull, "): ", x$RMSEAClose, '\n', sep='')
}
fitStatistics <- function(model, retval) {
datalist <- generateDataList(model)
likelihood <- retval[['Minus2LogLikelihood']]
saturated <- retval[['SaturatedLikelihood']]
independence <- retval[['IndependenceLikelihood']]
if (is.null(model@output$fitUnits)) return(retval)
if (model@output$fitUnits=="-2lnL" && is.null(model@output$chi)) {
chi <- likelihood - saturated
} else {chi <- model@output$chi}
if (is.null(chi)) return(retval)
DoF <- retval$degreesOfFreedom
satDoF <- retval$saturatedDoF
indDoF <- retval$independenceDoF
nParam <- dim(retval$parameters)[1]
Fvalue <- computeFValue(datalist, likelihood, chi)
if (model@output$fitUnits=="-2lnL" && is.null(model@output$chiDoF)) {
chiDoF <- DoF - satDoF # DoF = obsStat-model.ep; satDoF = obsStat-sat.ep; So sat.ep-model.ep == DoF-satDoF
} else {chiDoF <- model@output$chiDoF}
retval[['ChiDoF']] <- chiDoF
retval[['Chi']] <- chi
retval[['p']] <- suppressWarnings(ifelse(chiDoF==0,1.0,pchisq(chi, chiDoF, lower.tail = FALSE)))
retval[['AIC.Mx']] <- Fvalue - 2 * DoF
retval[['BIC.Mx']] <- (Fvalue - DoF * log(retval[['numObs']]))
AIC.p <- Fvalue + 2 * nParam
BIC.p <- (Fvalue + nParam * log(retval[['numObs']]))
sBIC <- (Fvalue + nParam * log((retval[['numObs']]+2)/24))
AICc <- Fvalue + 2*nParam + (2*nParam*(nParam+1))/(retval[['numObs']]-nParam-1)
retval[['satDoF']] <- satDoF
retval[['indDoF']] <- indDoF
IC <- matrix(NA, nrow=2, ncol=3, dimnames=list(c("AIC:", "BIC:"), c('df', 'par', 'sample')))
IC[,'df'] <- c(retval$AIC.Mx, retval$BIC.Mx)
IC[,'par'] <- c(AIC.p, BIC.p)
IC['BIC:','sample'] <- sBIC
IC['AIC:','sample'] <- AICc
if(length(model@output) && (is.null(model@output$fitUnits) || model@output$fitUnits=="-2lnL")){
retval[['informationCriteria']] <- IC
}
retval$fitUnits <- model@output$fitUnits
retval$fit <- if(retval$fitUnits == '-2lnL'){ retval[['Minus2LogLikelihood']] } else if(retval$fitUnits %in% c("r'Wr", "r'wr")) {retval[['Chi']]}
fi <- computeFitStatistics(likelihood, DoF, chi, chiDoF,
retval[['numObs']], independence, indDoF, saturated, satDoF)
for (k in names(fi)) retval[[k]] <- fi[[k]]
return(retval)
}
omxRMSEA <- function(model, lower=.025, upper=.975, null=.05, ...){
smod <- summary(model, ...)
x2 <- smod$Chi
df <- smod$ChiDoF
N <- smod$numObs
rmsea <- smod$RMSEA
ci <- rmseaConfidenceIntervalHelper(chi.squared=x2, df=df, N=N, lower=lower, upper=upper)
pn <- rmseaPCloseHelper(x2, df, N, null)
return(c(ci[1], est.rmsea=rmsea, ci[2], null=null, `Prob(x <= null)`=pn))
}
rmseaConfidenceIntervalHelper <- function(chi.squared, df, N, lower, upper){
# Lower confidence interval
if( pchisq(chi.squared, df=df, ncp=0) >= upper){ #sic
lower.lam <- uniroot(f=pChiSqFun, interval=c(1e-10, 1e4), val=chi.squared,
degf=df, goal=upper, extendInt="upX", maxiter=100L)$root
# solve pchisq(ch, df=df, ncp=x) == upper for x
} else{
lower.lam <- 0
}
# Upper confidence interval
if( pchisq(chi.squared, df=df, ncp=0) >= lower){ #sic
upper.lam <- uniroot(f=pChiSqFun, interval=c(1e-10, 1e4), val=chi.squared,
degf=df, goal=lower, extendInt="upX", maxiter=100L)$root
# solve pchisq(ch, df=df, ncp=x) == lower for x
} else{
upper.lam <- 0
}
lower.rmsea <- sqrt(lower.lam/(N*df))
upper.rmsea <- sqrt(upper.lam/(N*df))
return(c(lower=lower.rmsea, upper=upper.rmsea))
}
pChiSqFun <- function(x, val, degf, goal){
goal - pchisq(val, degf, ncp=x)
}
rmseaPCloseHelper <- function(x2, df, N, null){
1-pchisq(x2, df=df, ncp=N*df*(null)^2)
}
parameterList <- function(model, flatModel) {
parameterListHelper(model, flatModel, FALSE)
}
parameterListHelper <- function(model, flatModel, withModelName) {
ptable <- data.frame()
if(length(model@output) == 0) { return(ptable) }
estimates <- model@output$estimate
errorEstimates <- rep.int(as.numeric(NA), length(estimates))
if (!is.null(model@output$standardErrors)) {
se <- model@output$standardErrors
errorEstimates[match(rownames(se), names(estimates))] <- se
}
matrices <- generateMatrixList(flatModel)
parameters <- flatModel@parameters
if (length(estimates) > 0) {
matrixNames <- names(matrices)
for(i in 1:length(estimates)) {
mLocation <- parameters[[i]][[5]][[1]] + 1
mRow <- parameters[[i]][[5]][[2]] + 1
mCol <- parameters[[i]][[5]][[3]] + 1
mat <- matrices[[mLocation]][[1]]
dn <- dimnames(mat)
if (!is.null(dn[[1]])) mRow <- dn[[1]][mRow]
if (!is.null(dn[[2]])) mCol <- dn[[2]][mCol]
lbound <- parameters[[i]][[1]]
ubound <- parameters[[i]][[2]]
if (withModelName) {
ptable[i, 'model'] <- model@name
}
if (!is.null(names(estimates))) {
ptable[i, 'name'] <- names(estimates)[[i]]
} else {
ptable[i, 'name'] <- paste("p", i, sep="")
}
ptable[i, 'matrix'] <- simplifyName(matrixNames[[mLocation]], model@name)
ptable[i, 'row'] <- mRow
ptable[i, 'col'] <- mCol
ptable[i, 'Estimate'] <- estimates[[i]]
ptable[i, 'Std.Error'] <- errorEstimates[[i]]
ptable[i, 'lbound'] <- lbound
ptable[i, 'ubound'] <- ubound
}
}
return(ptable)
}
computeOptimizationStatistics <- function(model, flatModel, numStats, saturatedDoF, independenceDoF, retval) {
datalist <- flatModel@datasets
labelsData <- imxGenerateLabels(model)
expectations <- convertExpectationFunctions(flatModel, model, labelsData, new("MxDirectedGraph"))
# get estimated parameters
estimates <- model@output$estimate
# should saturated/independence models include means?
if(length(datalist)==1){
type <- datalist[[1]]@type
means <- datalist[[1]]@means
# if there's raw data, then use means in saturated/independence models
if(type=="raw"){
useMeans <- TRUE
} else {
# if there's not raw data, only use means if they're present
if((dim(means)[2]==1)&is.na(means[1,1])){
useMeans <- FALSE
} else{
useMeans <- TRUE
}
}
# number of variables
if(datalist[[1]]@type != 'raw'){
if(datalist[[1]]@type == 'none'){
nvar <- nrow(datalist[[1]]@observedStats$cov)
}
else{
nvar <- dim(datalist[[1]]@observed)[2]
}
} else if( length(expectations) == 1 ) {
nvar <- length(expectations[[1]]@dims)
} else {
nvar <- 0
}
# if there are multiple or zero datalists, then do nothing
} else {
useMeans <- NA
nvar <- 0
}
# how many thresholds does each variable have (needed for saturated and independence DoF calculation)
obj <- expectations
# grab the thresholdLevels object and expected means; punt if there is more than one expectation
if (length(obj)==1){
if ("thresholdLevels" %in% slotNames(obj[[1]])){
thresholdLevels <- obj[[1]]@thresholdLevels
if (length(thresholdLevels)==0){thresholdLevels <- rep(NA, nvar)}
} else {
thresholdLevels <- rep(NA, nvar)
}
} else {
thresholdLevels <- NULL
}
# number of continuous variables, provided there is just one expectation
if (!is.null(thresholdLevels)){
continuous <- sum(is.na(thresholdLevels))
} else{
continuous <- NA
}
# number of thresholds in the model
if (!is.null(thresholdLevels)){
thresh <- sum(thresholdLevels, na.rm=TRUE)
} else{
thresh <- NA
}
# constraints, parameters, model degrees of freedom
retval[['constraints']] <- calculateConstraints(model, flatModel)
retval[['estimatedParameters']] <- nrow(retval$parameters)
if(any(sapply(obj,function(x){is(x,"MxExpectationGREML")}))){
retval[['estimatedParameters']] <- retval[['estimatedParameters']] +
# When REML=FALSE and length(yhat)>0, numFixEff will be zero:
sum(sapply(obj,imxExtractSlot,name="numFixEff"))
}
if (is.null(numStats)) {
retval[['observedStatistics']] <- observedStatistics(model, flatModel, sum(retval$constraints))
} else {
retval[['observedStatistics']] <- numStats
}
retval[['degreesOfFreedom']] <- retval$observedStatistics - retval$estimatedParameters
# calculate or populate saturated degrees of freedom
if(is.null(saturatedDoF)) {
retval[['saturatedDoF']] <- retval$observedStatistics - (nvar * (nvar-1) / 2 + continuous*(1+useMeans) + thresh)
} else {
retval[['saturatedDoF']] <- saturatedDoF
}
#The "saturated model" has no sensible definiton with GREML expectation:
if(any(sapply(obj,function(x){is(x,"MxExpectationGREML")}))){
retval[['saturatedDoF']] <- NA
}
# calculate or populate independence degrees of freedom
if(is.null(independenceDoF)) {
if(!any(sapply(obj,function(x){is(x,"MxExpectationGREML")}))){
# indDoF = 1 df per continuous variable variance + 1 df per continuous mean + 1 df per threshold
retval[['independenceDoF']] <- retval$observedStatistics - (continuous*(1+useMeans) + thresh)
} else{
#TODO: the GREML expectation doesn't currently have a way to know how many phenotypes there are in every case.
#For now, leave the GREML independence model undefined
# #With GREML expectation, the independence model has a variance for each phenotype, and the same fixed effects as the fitted model:
# retval[['independenceDoF']] <-
# retval$observedStatistics - sum(sapply(obj,function(x){length(x@yvars)})) - sum(sapply(obj,imxExtractSlot,name="numFixEff"))
retval[['independenceDoF']] <- NA
}
} else {
retval[['independenceDoF']] <- independenceDoF
}
# set NULLs to NAs
if (is.null(retval$saturatedDoF)) {
retval$SaturatedDoF <- NA
}
if (is.null(retval$independenceDoF)) {
retval$IndependenceDoF <- NA
}
retval[['saturatedParameters']] <- retval[['observedStatistics']] - retval[['saturatedDoF']]
retval[['independenceParameters']] <- retval[['observedStatistics']] - retval[['independenceDoF']]
# calculate fit statistics
retval <- fitStatistics(model, retval)
return(retval)
}
print.summary.mxmodel <- function(x,...) {
cat("Summary of", x$modelName, '\n', '\n')
if(x$verbose==TRUE){
if (length(x$compute)) {
cat("compute plan:\n")
print(x$compute)
cat("\n")
}
if (length(x$dataSummary) > 0) {
cat("data:\n")
print(x$dataSummary)
}
}
if (!is.null(x$npsolMessage)) {
cat(x$npsolMessage,'\n','\n')
if ((x$statusCode == "not convex/red" || x$statusCode == "nonzero gradient/red") &&
(is.na(x$maxRelativeOrdinalError) || x$maxRelativeOrdinalError > 0)) {
cat(paste("Your ordinal model may converge if you reduce mvnRelEps\n",
" try: model <- mxOption(model, 'mvnRelEps', mxOption(model, 'mvnRelEps')/5)\n\n"))
}
}
params <- x$parameters
if (!is.null(params) && nrow(params)) {
cat("free parameters:\n")
params$lbound <- mapply(highlightProblem, params$lbound, params$lboundMet)
params$ubound <- mapply(highlightProblem, params$ubound, params$uboundMet)
params$lbound[is.na(params$lbound)] <- ""
params$ubound[is.na(params$ubound)] <- ""
params$lboundMet <- NULL
params$uboundMet <- NULL
if (!is.null(x$bootstrapSE) && length(x$bootstrapSE) == nrow(params)) {
params[['Std.Error']] <- x$bootstrapSE
}
if (!x$verbose) {
if (all(is.na(params[['Std.Error']]))) {
params[['Std.Error']] <- NULL
}
if (all(params$lbound == "" & params$ubound == "")) {
params[['lbound']] <- NULL
params[['ubound']] <- NULL
}
}
if (!is.null(params[['Std.Error']]) && !is.null(x[['seSuspect']])) {
seCol <- match('Std.Error', colnames(params))
before <- params[1:seCol]
stars <- mapply(highlightProblem, rep('',length(x[['seSuspect']])), x[['seSuspect']])
fullStars <- rep('', nrow(params))
fullStars[match(names(x[['seSuspect']]), params$name)] <- stars
if (length(params) > seCol) {
params <- cbind(before, 'A'=fullStars, params[(seCol+1):length(params)])
} else {
params <- cbind(before, 'A'=fullStars)
}
}
if (!is.null(x$bootstrapQuantile) && nrow(x$bootstrapQuantile) == nrow(params)) {
bq <- x$bootstrapQuantile
params <- cbind(params, bq)
}
cmap <- 1:ncol(params)
isBound <- colnames(params) %in% paste0(c('l','u'),'bound')
params <- params[,c(cmap[!isBound], cmap[isBound]),drop=FALSE]
print(params)
cat('\n')
}
if (!is.null(x$CI) && length(x$CI) > 0) {
cat("confidence intervals:\n")
print(x$CI)
if (any(is.na(x$CI[,c('lbound','ubound')])) && !(x$verbose)) {
cat(" To investigate missing CIs, run summary() again, with verbose=T, to see CI details.", '\n')
}
cat('\n')
} else if (x$CI.Requested) {
cat("To obtain confidence intervals re-run with intervals=TRUE\n\n")
}
if(x$verbose && length(x$CIdetail)){
cat("CI details:\n")
print(x$CIdetail)
cat("\n")
}
if(length(x$GREMLfixeff)>0 && any(sapply(x$GREMLfixeff,length)>0)){
cat("regression coefficients:\n")
print(x$GREMLfixeff)
cat("\n")
}
cat('Model Statistics:', '\n')
EP <- matrix(
c(x$estimatedParameters, x$saturatedParameters, x$independenceParameters,
x$degreesOfFreedom, x$saturatedDoF, x$independenceDoF,
x$Minus2LogLikelihood, x$SaturatedLikelihood, x$IndependenceLikelihood),
nrow=3, ncol=3,
dimnames=list(
c(' Model:', ' Saturated:', 'Independence:'),
c(' | Parameters', ' | Degrees of Freedom', paste0(' | Fit (', x$fitUnits, ' units)'))
)
)
print(EP)
cat('Number of observations/statistics: ', x$numObs, "/", x$observedStatistics, '\n\n', sep="")
constraints <- x$constraints
if(length(constraints) > 0) {
for(i in 1:length(constraints)) {
name <- names(constraints)[[i]]
if (constraints[[i]] == 1) plural <- ''
else plural <- 's'
cat("Constraint", omxQuotes(simplifyName(name, x$modelName)), "contributes",
constraints[[i]], paste("observed statistic", plural, '.', sep=''), "\n")
if(i==length(constraints)){cat("\n")}
}
}
if (!is.null(x$infoDefinite) && !is.na(x$infoDefinite)) {
if (!x$infoDefinite) {
cat("\n** Information matrix is not positive definite (not at a candidate optimum).\n Be suspicious of these results. At minimum, do not trust the standard errors.\n\n")
} else if(x$verbose==TRUE) {
cat("condition number of the information matrix: ", x$conditionNumber, "\n")
}
}
if (x$verbose==TRUE && !is.null(x$maxAbsGradient)) {
cat("maximum absolute gradient: ", x$maxAbsGradient, " (",names(x$maxAbsGradient),")\n")
}
#
# Chi-square goodness of fit test
if(x$verbose==TRUE || (!is.null(x$Chi) && !is.na(x$Chi))) {
chival <- x$Chi
if(is.na(x$SaturatedLikelihood)){
chidof <- NA
} else {
chidof <- x$ChiDoF
}
chipee <- x$p
cat("chi-square: ", "\U03C7\U00B2 ( df=", chidof, " ) = ", chival, ", p = ", chipee, '\n', sep="")
}
#
if(length(x$informationCriteria)){
cat("Information Criteria: \n")
IC <- x$informationCriteria
colnames(IC) <- c(" | df Penalty", " | Parameters Penalty", " | Sample-Size Adjusted")
print(IC)
}
#
# Absolute fit indices
if(x$verbose==TRUE || any(!is.na(c(x$CFI, x$TLI, x$RMSEA)))){
catFitStatistics(x)
}
if(any(is.na(c(x$CFI, x$TLI, x$RMSEA)))){
cat("To get additional fit indices, see help(mxRefModels)\n")
}
#
# Timing information
cat("timestamp:", format(x$timestamp), '\n')
if(x$verbose==TRUE){
cat("frontend time:", format(x$frontendTime), '\n')
cat("backend time:", format(x$backendTime), '\n')
cat("independent submodels time:", format(x$independentTime), '\n')
cat("cpu time:", format(x$cpuTime), '\n')
}
cat("Wall clock time:", format(x$wallTime), "\n")
if (x$verbose==FALSE && !is.null(x$optimizerEngine)) {
cat("optimizer: ", x$optimizerEngine, '\n')
}
cat("OpenMx version number:", format(x$mxVersion), '\n')
cat("Need help? See help(mxSummary)", '\n')
cat('\n')
if (!x$wasRun) {
message("WARNING: This model has not been run yet. Tip: Use\n model = mxRun(model)\nto estimate a model.")
} else if (x$stale) {
message("WARNING: This model was modified since it was run. Summary information may be out-of-date.")
}
}
setLikelihoods <- function(model, saturatedLikelihood, independenceLikelihood, retval) {
# populate saturated -2 log likelihood
if(is.null(saturatedLikelihood)) {
retval$SaturatedLikelihood <- attr(model@fitfunction$result, "SaturatedLikelihood")
} else {
retval$SaturatedLikelihood <- saturatedLikelihood
}
# populate independence -2 log likelihood
if(is.null(independenceLikelihood)) {
retval$IndependenceLikelihood <- attr(model@fitfunction$result, "IndependenceLikelihood")
} else {
retval$IndependenceLikelihood <- independenceLikelihood
}
# populate model -2 log likelihood
retval$Minus2LogLikelihood <- model@output$Minus2LogLikelihood
# set NULLs to NAs
if (is.null(retval$SaturatedLikelihood)) {
retval$SaturatedLikelihood <- NA
}
if (is.null(retval$Minus2LogLikelihood)) {
retval$Minus2LogLikelihood <- NA
}
if (is.null(retval$IndependenceLikelihood)) {
retval$IndependenceLikelihood <- NA
}
return(retval)
}
setNumberObservations <- function(numObs, datalist, fitfunctions, retval) {
if(is.null(numObs)) {
retval$numObs <- numberObservations(datalist, fitfunctions)
} else {
retval$numObs <- numObs
}
return(retval)
}
generateDataSummary <- function(model) {
datalist <- generateDataList(model)
retval <- lapply(datalist, summarize)
return(retval)
}
##' imxEvalByName
##'
##' This is an internal function exported for those people who know
##' what they are doing.
##'
##' @param name name
##' @param model model
##' @param compute compute
##' @param show show
imxEvalByName <- function(name, model, compute=FALSE, show=FALSE) {
if ((length(name) != 1) || typeof(name) != "character") {
stop("'name' argument must be a character argument")
}
if (!is(model, "MxModel")) {
stop("'model' argument must be a MxModel object")
}
if (hasSquareBrackets(name)) {
components <- splitSubstitution(name)
eval(substitute(mxEval(x[y,z], model, compute, show),
list(x = as.name(components[[1]]),
y = parse(text=components[[2]])[[1]],
z = parse(text=components[[3]])[[1]])))
} else {
eval(substitute(mxEval(x, model, compute, show),
list(x = as.name(name))))
}
}
boundsMet <- function(model, retval){
params <- retval$parameters
lbound <- params$lbound
ubound <- params$ubound
estimate <- params$Estimate
threshold <- model@options[["Feasibility tolerance"]]
if (is.null(threshold)){
threshold <- getOption("mxOptions")[["Feasibility tolerance"]]
}
threshold <- as.numeric(threshold)
lboundMet <- mapply(compareBounds, estimate, lbound, MoreArgs=list(threshold))
uboundMet <- mapply(compareBounds, estimate, ubound, MoreArgs=list(threshold))
params$lboundMet <- lboundMet
params$uboundMet <- uboundMet
retval$parameters <- params
return(retval)
}
compareBounds <- function(estimate, bound, threshold){
if (is.na(bound)){
return(FALSE)
}
absDelta <- abs(estimate - bound)
return (absDelta < threshold)
}
highlightProblem <- function(bound, boundMet){
if (boundMet){
if (is.numeric(bound)) bound <- round(bound,4)
return(paste(bound, "!", sep=""))
}
else {
return(bound)
}
}
parseLikelihoodArg <- function(input, arg) {
input <- input[[arg]]
if (is.null(input)) {
return(input)
} else if (is.numeric(input)) {
return(input)
} else if (is(input, "MxModel")) {
name <- input@name
if (is.null(input@fitfunction)) {
stop(paste(omxQuotes(name), "model passed",
"to summary function does not",
"have top-level fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
if (length(input@fitfunction@result) != 1) {
stop(paste(omxQuotes(name), "model passed to summary",
"function does not have a 1x1 matrix",
"result in fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
return(input@fitfunction@result[1,1])
} else if(is.list(input) && length(input)==2) {
stop(paste("List of length two (illegal argument) passed to", omxQuotes(arg),
"argument of summary function. You probably meant to use",
"the refModels argument instead."), call. = FALSE)
} else {
stop(paste("Illegal argument passed to", omxQuotes(arg),
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
parseDfArg <- function(input, arg) {
input <- input[[arg]]
if (is.null(input)) {
return(input)
} else if (is.numeric(input)) {
return(input)
} else if (is(input, "MxModel")) {
name <- input@name
if (is.null(input@fitfunction)) {
stop(paste(omxQuotes(name), "model passed",
"to summary function does not",
"have top-level fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
if (length(input@fitfunction@result) != 1) {
stop(paste(omxQuotes(name), "model passed to summary",
"function does not have a 1x1 matrix",
"result in fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
return(summary(input)$degreesOfFreedom)
} else {
stop(paste("Illegal argument passed to", omxQuotes(arg),
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
refToLikelihood <- function(model) {
if (is(model, "MxModel")) {
if (!model@.wasRun) stop("Reference model must be run to obtain fit indices")
model$output$Minus2LogLikelihood
} else if (is.list(model)) {
model[[1]]
} else {
stop(paste("Illegal argument passed to refModels",
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
refToDof <- function(model) {
if (is(model, "MxModel")) {
if (!model@.wasRun) stop("Reference model must be run to obtain fit indices")
return(summary(model)$degreesOfFreedom)
} else if (is.list(model)) {
model[[2]]
} else {
stop(paste("Illegal argument passed to refModels",
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
summarizeBootstrap <- function(mle, bootData, bq, summaryType) {
if (summaryType == 'quantile') {
t(apply(bootData, 2, quantile, probs=bq, type=1))
} else if (summaryType == 'bcbci') {
zcrit <- qnorm(bq)
out <- matrix(NA, nrow=length(mle), ncol=length(bq),
dimnames=list(names(mle),
sapply(bq, function(x) sprintf("%.1f%%", round(100*min(x), 1)))))
for(i in 1:length(mle)) {
z0 <- qnorm(mean(bootData[,i] <= mle[i]))
for (qx in 1:length(bq)) {
phi <- pnorm(2*z0 + zcrit[qx])
out[i,qx] <- quantile(bootData[,i], probs=phi, type=1)
}
}
out
} else {
warning(paste("boot.SummaryType =", omxQuotes(summaryType),
"is not recognized"))
}
}
summary.MxModel <- function(object, ..., verbose=FALSE) {
out <- try(summary.MxModel.Impl(object, ..., verbose=verbose), silent=TRUE)
if (is(out, "try-error")) {
if (!object@.wasRun) {
# Maybe should never happen?
stop("This model has not been run yet. Tip: Use\n model = mxRun(model)\nto estimate a model.")
} else if (object@.modifiedSinceRun) {
stop(paste("This model is in an inconsistent state.",
"Tip: Use\n model = mxRun(model)\nto re-estimate the model."))
} else {
stop(out)
}
}
out
}
summary.MxModel.Impl <- function(object, ..., verbose) {
model <- object
dotArguments <- list(...)
if (!is.null(dotArguments[["refModels"]])) {
refModels <- dotArguments[["refModels"]]
satModel <- refModels[['Saturated']]
indModel <- refModels[['Independence']]
saturatedLikelihood <- refToLikelihood(satModel)
saturatedDoF <- refToDof(satModel)
independenceLikelihood <- refToLikelihood(indModel)
independenceDoF <- refToDof(indModel)
} else {
saturatedLikelihood <- parseLikelihoodArg(dotArguments, "SaturatedLikelihood")
saturatedDoF <- parseDfArg(dotArguments, "SaturatedDoF")
independenceLikelihood <- parseLikelihoodArg(dotArguments, "IndependenceLikelihood")
independenceDoF <- parseDfArg(dotArguments, "IndependenceDoF")
}
numObs <- dotArguments$numObs
numStats <- dotArguments$numStats
namespace <- imxGenerateNamespace(model)
flatModel <- imxFlattenModel(model, namespace, TRUE)
flatModel <- constraintsToAlgebras(flatModel)
flatModel <- eliminateObjectiveFunctions(flatModel)
convertArguments <- imxCheckVariables(flatModel, namespace)
flatModel <- convertAlgebras(flatModel, convertArguments)
labelsData <- imxGenerateLabels(model)
flatModel <- expectationFunctionConvertEntities(flatModel, namespace, labelsData)
flatModel <- populateDefInitialValues(flatModel)
flatModel <- checkEvaluation(model, flatModel)
dependencies <- cycleDetection(flatModel)
dependencies <- transitiveClosure(flatModel, dependencies)
freeVarGroups <- buildFreeVarGroupList(flatModel)
flatModel <- generateParameterList(flatModel, dependencies, freeVarGroups)
retval <- list(wasRun=model@.wasRun, stale=model@.modifiedSinceRun)
retval$parameters <- parameterList(model, flatModel)
if (!is.null(model@compute$steps[['ND']]) && model@compute$steps[['ND']]$checkGradient &&
!is.null(model@compute$steps[['ND']]$output$gradient)) {
gdetail <- model@compute$steps[['ND']]$output$gradient
retval$seSuspect <- !gdetail[,'symmetric']
names(retval$seSuspect) <- rownames(gdetail)
}
if (is(model@compute, "MxComputeBootstrap")) {
bq <- c(.25,.75)
if (!is.null(dotArguments[["boot.quantile"]])) {
bq <- sort(as.numeric(dotArguments[["boot.quantile"]]))
}
summaryType <- 'bcbci'
if (!is.null(dotArguments[["boot.SummaryType"]])) {
summaryType <- dotArguments[["boot.SummaryType"]]
}
cb <- model@compute
if (!is.null(cb@output$raw) && is.na(cb@only) && cb@output$numParam == nrow(retval$parameters)) {
raw <- cb@output$raw
mask <- raw[,'statusCode'] %in% cb@OK
bootData <- raw[mask, 3:(nrow(retval$parameters)+2), drop=FALSE]
if (sum(mask) < .95*nrow(raw)) {
pct <- round(100*sum(mask) / nrow(raw))
warning(paste0("Only ",pct,"% of the bootstrap replications ",
"converged. Accuracy is much less than the ", nrow(raw),
" replications requested"), call.=FALSE)
}
if (sum(mask) >= 3) {
retval$bootstrapSE <- apply(bootData, 2, sd)
retval$bootstrapQuantile <-
summarizeBootstrap(retval$parameters[, 'Estimate'], bootData, bq, summaryType)
}
}
} else if (any(grep('^boot\\.', names(dotArguments)))) {
warning("No bootstrap data found. See ?mxBootstrap")
}
retval$GREMLfixeff <- GREMLFixEffList(model)
retval$infoDefinite <- model@output$infoDefinite
retval$conditionNumber <- model@output$conditionNumber
if (length(model@output$gradient)) {
agrad <- abs(model@output$gradient)
retval$maxAbsGradient <- agrad[ order(-agrad)[1] ]
}
retval <- boundsMet(model, retval)
retval <- setLikelihoods(model, saturatedLikelihood, independenceLikelihood, retval)
labelsData <- imxGenerateLabels(model)
fitfunctions <- convertFitFunctions(flatModel, model, labelsData, dependencies)
retval <- setNumberObservations(numObs, flatModel@datasets, fitfunctions, retval)
retval <- computeOptimizationStatistics(model, flatModel, numStats, saturatedDoF, independenceDoF, retval)
retval$dataSummary <- generateDataSummary(model)
retval$CI.Requested <- length(model@intervals) > 0
retval$CI <- as.data.frame(model@output$confidenceIntervals)
if (length(retval$CI) && nrow(retval$CI)) {
retval$CI <- cbind(retval$CI, note=apply(retval$CI, 1, function(ci) {
# This should probably take into account whether both bounds
# were requested and consider the optimizer codes also. TODO
if (any(is.na(ci)) || ci[1] == ci[3] || ci[1] >= ci[2] || ci[2] >= ci[3]) {
"!!!"
} else {
""
}
}))
}
retval$CIcodes <- model@output$confidenceIntervalCodes
statusCode <- model@output$status$code
if (!is.null(statusCode)) {
message <- optimizerMessages[[as.character(statusCode)]]
retval[['npsolMessage']] <- message
retval[['statusCode']] <- as.statusCode(statusCode)
retval[['maxRelativeOrdinalError']] <- model@output[['maxRelativeOrdinalError']]
}
if( .hasSlot(model,"compute") && length(model$compute$steps$CI) ){
retval$CIdetail <- model$compute$steps$CI$output$detail
}
retval$timestamp <- model@output$timestamp
retval$frontendTime <- model@output$frontendTime
retval$backendTime <- model@output$backendTime
retval$independentTime <- model@output$independentTime
retval$wallTime <- model@output$wallTime
retval$cpuTime <- model@output$cpuTime
retval$mxVersion <- model@output$mxVersion
retval$modelName <- model@name
plan <- model@compute
if (is(plan, "MxComputeSequence")) {
gd <- plan$steps[['GD']]
if (is(gd, "MxComputeGradientDescent")) {
retval$optimizerEngine <- gd$engine
}
}
retval$verbose <- verbose
class(retval) <- "summary.mxmodel"
return(retval)
}
assertModelRunAndFresh <- function(model) {
warnModelCreatedByOldVersion(model)
if (.hasSlot(model,".wasRun") && !model@.wasRun) stop("This model has not been run yet. Tip: Use\n model = mxRun(model)\nto estimate a model.")
if (.hasSlot(model,".modifiedSinceRun") && model@.modifiedSinceRun) {
msg <- paste("MxModel", omxQuotes(model@name), "was modified",
"since it was run.")
warning(msg)
}
}
assertModelFreshlyRun <- function(model) {
warnModelCreatedByOldVersion(model)
if (model@.wasRun && model@.modifiedSinceRun) {
msg <- paste("MxModel", omxQuotes(model@name), "was modified",
"since it was run.")
warning(msg)
}
}
logLik.MxModel <- function(object, ...) {
model <- object
moreModels <- list(...)
assertModelFreshlyRun(model)
ll <- NA
if (length(model@output) && !is.null(model@output$fit) &&
!is.null(model@output$fitUnits) ) {
if(model@output$fitUnits=="-2lnL"){
ll <- -0.5*model@output$fit
} else if(model@output$fitUnits %in% c("r'Wr", "r'wr")) {
ll <- -0.5*model@output$chi
}
#TODO: this doesn't count "implicit" free parameters that are "profiled out":
attr(ll, "df") <- length(model@output$estimate)
} else {
attr(ll,"df") <- NA
}
namespace <- imxGenerateNamespace(model)
flatModel <- imxFlattenModel(model, namespace, TRUE)
labelsData <- imxGenerateLabels(model)
fitfunctions <- convertFitFunctions(flatModel, model, labelsData, new("MxDirectedGraph"))
nobs <- numberObservations(flatModel@datasets, fitfunctions)
if (!is.na(nobs)) {
attr(ll,"nobs") <- nobs
}
class(ll) <- "logLik"
if (length(moreModels)) {
c(list(ll), lapply(moreModels, logLik.MxModel))
} else {
ll
}
}
.standardizeParams <- function(x=NULL, model, Apos, Spos, Mpos, give.matrices=FALSE){
if(is.null(x)){x <- omxGetParameters(model)}
param <- omxGetParameters(model)
paramNames <- names(param)
model <- omxSetParameters(model, values=x, labels=paramNames, free=TRUE)
model_A <- model[[model$expectation$A]] #<--the A matrix might not be named "A".
A <- list( model_A$values, model_A$result )
A <- A[[which.max(c( length(A[[1]]), length(A[[2]]) ))]]
model_S <- model[[model$expectation$S]] #<--Likewise for S
S <- list( model_S$values, model_S$result )
S <- S[[which.max(c( length(S[[1]]), length(S[[2]]) ))]]
M <- NULL
if(!is.null(Mpos)){
model_M <- model[[model$expectation$M]]
M <- list( model_M$values, model_M$result )
M <- M[[which.max(c( length(M[[1]]), length(M[[2]]) ))]]
}
else{model_M <- NULL}
I <- diag(1, nrow(A))
ImAInv <- solve(I-A)
SD <- sqrt(diag(ImAInv %*% S %*% t(ImAInv)))
SD <- diag(SD,nrow=length(SD)) #<--Needed if line immediately above is a scalar.
InvSD <- 1/diag(SD)
InvSD <- diag(InvSD,nrow=length(InvSD))
Az <- InvSD %*% A %*% SD
Sz <- InvSD %*% S %*% InvSD
if(!is.null(M)){
Mz <- M %*% InvSD
sparam <- c(Az[!is.na(Apos)],Sz[!is.na(Spos)],Mz[!is.na(Mpos)])
names(sparam) <- c(Apos[!is.na(Apos)],Spos[!is.na(Spos)],Mpos[!is.na(Mpos)])
if(!give.matrices){return(sparam)}
else{return(list(sparam=sparam,Az=Az,Sz=Sz,Mz=Mz))}
}
else{
sparam <- c(Az[!is.na(Apos)],Sz[!is.na(Spos)])
names(sparam) <- c(Apos[!is.na(Apos)],Spos[!is.na(Spos)])
if(!give.matrices){return(sparam)}
else{return(list(sparam=sparam,Az=Az,Sz=Sz))}
}
}
.mxStandardizeRAMhelper <- function(model,SE=FALSE,ParamsCov,inde.subs.flag=FALSE,ignoreSubmodels=FALSE){
#Recur the function for the appropriate submodels, if any:
if(length(model@submodels) && !ignoreSubmodels){
return(lapply(
model@submodels[which(
sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0})
)],
.mxStandardizeRAMhelper,SE=SE,ParamsCov=ParamsCov,inde.subs.flag=inde.subs.flag,ignoreSubmodels=FALSE))
}
#Get A and S:
model_A <- model[[model$expectation$A]] #<--Necessary because the A matrix might not be named "A".
A <- list( model_A$values, model_A$result )
A <- A[[which.max(c( length(A[[1]]), length(A[[2]]) ))]]
model_S <- model[[model$expectation$S]] #<--Likewise for S
S <- list( model_S$values, model_S$result )
S <- S[[which.max(c( length(S[[1]]), length(S[[2]]) ))]]
M <- NULL
if(!is.na(model$expectation$M)){
model_M <- model[[model$expectation$M]]
M <- list( model_M$values, model_M$result )
M <- M[[which.max(c( length(M[[1]]), length(M[[2]]) ))]]
}
#Find positions of nonzero paths:
Apos <- matrix(NA,nrow=nrow(A),ncol=ncol(A),dimnames=dimnames(A))
Spos <- matrix(NA,nrow=nrow(S),ncol=ncol(S),dimnames=dimnames(S))
A_need_pos <- which(A!=0,arr.ind=T)
S_need_pos <- which(S!=0,arr.ind=T)
S_need_pos <- subset(S_need_pos, S_need_pos[,1]>=S_need_pos[,2]) #<--Lower tri only
if(!is.null(M)){
Mpos <- matrix(NA,nrow=1,ncol=ncol(M),dimnames=dimnames(M))
M_need_pos <- which(M!=0)
}
else{
Mpos <- NULL
M_need_pos <- NULL
}
numelem <- nrow(A_need_pos)+nrow(S_need_pos)+length(M_need_pos)
#Create output object:
out <- data.frame(name=vector(mode="character",length=numelem),label=vector(mode="character",length=numelem),
matrix=vector(mode="character",length=numelem),
row=vector(mode="character",length=numelem),col=vector(mode="character",length=numelem),
Raw.Value=vector(mode="numeric",length=numelem),
Raw.SE=vector(mode="numeric",length=numelem),
Std.Value=vector(mode="numeric",length=numelem),
Std.SE=vector(mode="numeric",length=numelem),stringsAsFactors=FALSE)
out$label <- NA
j <- 1
#Create position strings where needed and begin to populate output:
if(nrow(A_need_pos)>0){
for(i in 1:nrow(A_need_pos)){
Apos[A_need_pos[i,1],A_need_pos[i,2]] <- out$name[j] <- paste(
model@name,".A[",A_need_pos[i,1],",",A_need_pos[i,2],"]",sep="")
if(!isAllNa(model_A$labels)){
out$label[j] <- model_A$labels[A_need_pos[i,1],A_need_pos[i,2]]
}
out$matrix[j] <- "A"
out$Raw.Value[j] <- A[A_need_pos[i,1],A_need_pos[i,2]]
out$row[j] <- ifelse(length(rownames(A))>0,rownames(A)[A_need_pos[i,1]],A_need_pos[i,1])
out$col[j] <- ifelse(length(colnames(A))>0,colnames(A)[A_need_pos[i,2]],A_need_pos[i,2])
j <- j+1
}
}
if(nrow(S_need_pos)>0){
for(i in 1:nrow(S_need_pos)){
Spos[S_need_pos[i,1],S_need_pos[i,2]] <- out$name[j] <- paste(
model@name,".S[",S_need_pos[i,1],",",S_need_pos[i,2],"]",sep="")
if(!isAllNa(model_S$labels)){
out$label[j] <- model_S$labels[S_need_pos[i,1],S_need_pos[i,2]]
}
out$matrix[j] <- "S"
out$Raw.Value[j] <- S[S_need_pos[i,1],S_need_pos[i,2]]
out$row[j] <- ifelse(length(rownames(S))>0,rownames(S)[S_need_pos[i,1]],S_need_pos[i,1])
out$col[j] <- ifelse(length(colnames(S))>0,colnames(S)[S_need_pos[i,2]],S_need_pos[i,2])
j <- j+1
}
}
if(length(M_need_pos)){
for(i in 1:length(M_need_pos)){
Mpos[1,M_need_pos[i]] <- out$name[j] <- paste(
model@name,".M[1,",M_need_pos[i],"]",sep="")
if(!is.null(model_M$labels) && !isAllNa(model_M$labels)){
out$label[j] <- model_M$labels[1,M_need_pos[i]]
}
out$matrix[j] <- "M"
out$Raw.Value[j] <- M[1,M_need_pos[i]]
out$row[j] <- 1
out$col[j] <- ifelse(length(colnames(M))>0,colnames(M)[M_need_pos[i]],M_need_pos[i])
j <- j+1
}
}
#Get standardized values:
freeparams <- omxGetParameters(model)
paramnames <- names(freeparams)
zout <- .standardizeParams(x=freeparams,model=model,Apos=Apos,Spos=Spos,Mpos=Mpos)
#Compute SEs, or assign them 'not requested' values, as the case may be:
if(SE){
if(!all(paramnames %in% rownames(ParamsCov))){
stop(paste(
"some free parameter labels do not appear in the dimnames of the parameter estimates' covariance matrix;\n",
"are you running mxStandardizeRAMpaths() on a dependent submodel instead of on its multigroup container model?\n",
"the missing parameter labels are:\n",
paste(paramnames[!(paramnames %in% rownames(ParamsCov))],collapse=", "),sep=""))
}
#From Mike Hunter's delta method example:
covParam <- ParamsCov[paramnames,paramnames,drop=FALSE]#<--submodel will usually not contain all free param.s
jacStand <- numDeriv::jacobian(func=.standardizeParams, x=freeparams, model=model, Apos=Apos, Spos=Spos, Mpos=Mpos)
covSparam <- jacStand %*% covParam %*% t(jacStand)
dimnames(covSparam) <- list(names(zout),names(zout))
SEs <- sqrt(diag(covSparam))
#SEs[diag(covSparam)<.Machine$double.eps] <- 0
}
else{SEs <- rep("not_requested",length(zout)); names(SEs) <- names(zout)}
#Add standardized values and SEs to output:
out$Std.Value <- zout
out$Std.SE <- SEs
#Pull in raw SEs if requested:
if(SE){
for(i in 1:numelem){
if( (out$name[i] %in% paramnames) |
(out$label[i] %in% paramnames) ){
tdiags <- covParam[ifelse(is.na(out$label[i]),out$name[i],out$label[i]),
ifelse(is.na(out$label[i]),out$name[i],out$label[i])]
if (length(tdiags) == 1) {
# For diag, R will return a square identity matrix of size given by the scalar
if(tdiags < 0 || is.na(tdiags)) {
warning("Some diagonal elements of the repeated-sampling covariance matrix of the point estimates are less than zero or NA.\nThat's weird. Raise an eyebrow at these standard errors.")
}
} else {
if(any(diag(tdiags) < 0) || any(is.na(tdiags))){
warning("Some diagonal elements of the repeated-sampling covariance matrix of the point estimates are less than zero or NA.\nThat's weird. Raise an eyebrow at these standard errors.")
}
}
out$Raw.SE[i] <- suppressWarnings(sqrt(tdiags))
}
else{
# # The elements of `out$name` should all reference scalars:
# tmp <- try(mxSE(x=out$name[i],model=model,forceName=TRUE,silent=TRUE))
# if(inherits(tmp,"try-error")){
# # Paths with the same label will have the same raw value and raw-value SE:
# tmp <- try(mxSE(x=out$label[i],model=model,forceName=TRUE,silent=TRUE))
# }
# if( !inherits(tmp,"try-error") && is.numeric(tmp) ){
# out$Raw.SE[i] <- tmp
# }
}
}
}
else{out$Raw.SE <- "not_requested"}
return(out)
}
mxStandardizeRAMpaths <- function(model, SE=FALSE, cov=NULL){
assertModelFreshlyRun(model)
if(imxHasDefinitionVariable(model)){
warning("'model' (or one of its submodels) contains definition variables; interpret results of mxStandardizeRAMpaths() cautiously")
}
#If SE=T,need to check for independent submodels because they will have their own Hessians;
#recur main function as appropriate:
inde.subs.flag <- FALSE
if(SE & length(model$submodels)>0){
RAM.subs <- (sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
inde.subs <- sapply(model@submodels,function(x){x@independent})==TRUE &
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
if(sum(inde.subs)>0){
out2 <- NULL
#if ALL submodels are either independent RAM models or non-RAM models:
if(all(RAM.subs==inde.subs)){
out <- lapply(model@submodels[which(inde.subs)],mxStandardizeRAMpaths,SE=T)
if(length(out)==0){stop(paste("model '",model@name,"' contains no submodels that use RAM expectation",sep=""))}
return(out)
}
else{out2 <- lapply(model@submodels[which(inde.subs)],mxStandardizeRAMpaths,SE=T)}
inde.subs.flag <- TRUE
}}
covParam <- NULL
if(SE){
#If user requests SEs and provided no covariance matrix, check to be sure SEs can and should be computed:
if(!length(cov)){
# if(length(model@constraints)>0){
# msg <- paste("standard errors will not be computed because model '",model@name,"' contains at least one mxConstraint",sep="")
# warning(msg)
# SE <- FALSE
# }
if(SE & length(model@output$vcov)==0){
if(!model@.wasRun){
msg <- paste("standard errors will not be computed because model '",model@name,"' has not yet been run, and no matrix was provided for argument 'cov'",sep="")
warning(msg)
SE <- FALSE
}
else{
warning("argument 'SE=TRUE' requires model to have a nonempty 'vcov' output slot, or a non-NULL value for argument 'cov'; continuing with 'SE' coerced to 'FALSE'")
SE <- FALSE
}}
libraries <- rownames(installed.packages())
pkgcheck <- ("numDeriv" %in% libraries)
if(SE & !pkgcheck){
warning("argument 'SE=TRUE' requires package 'numDeriv' to be installed; continuing with 'SE' coerced to 'FALSE'")
SE <- FALSE
}
if(SE){
covParam <- vcov(model)
}
}
#If user requests SEs and provided a covariance matrix:
else{
#Conceivably, the user could provide a sampling covariance matrix that IS valid in the presence of MxConstraints...
if(length(model@constraints)>0){
#msg <- paste("standard errors may be invalid because model '",model@name,"' contains at least one mxConstraint",sep="")
#warning(msg)
}
#Sanity checks on the value of argument 'cov':
if(!is.matrix(cov)){ #<--Is it a matrix?
cov <- try(as.matrix(cov),silent=T)
if("try-error" %in% class(cov) || !is.matrix(cov)){ #<--If its not a matrix, can it be coerced to one?
stop("non-NULL value to argument 'cov' must be (or be coercible to) a matrix")
}
}
if(nrow(cov)!=ncol(cov)){ #<--Is it square?
msg <- paste("non-NULL value to argument 'cov' must be a square matrix; it has ",nrow(cov)," rows and ",ncol(cov)," columns",sep="")
stop(msg)
}
#Do its row and column names match?:
if(!length(rownames(cov)) || !length(colnames(cov)) || any(rownames(cov) != colnames(cov))){
stop("non-NULL value to argument 'cov' must have matching and complete rownames and colnames")
}
paramnames <- names(omxGetParameters(model))
if(nrow(cov) != length(paramnames)){ #<--Do its dimensions match the number of free parameters?
msg <- paste("value of argument 'cov' has dimension ",nrow(cov),", but '",model@name,"' has ",length(paramnames)," free parameters",sep="")
stop(msg)
}
covnames <- colnames(cov)
#Do its row and column names match the free-parameter labels (ignoring permutations)?:
if(any(sort(covnames) != sort(paramnames))){
msg <- paste("the dimnames of the matrix provided for argument 'cov' do not match the free-parameter labels of '",model@name,"'",sep="")
stop(msg)
}
covParam <- cov[paramnames,paramnames]
}
}
#Check if single-group model uses RAM expectation, and proceed if so:
if(length(model@submodels)==0){
if (!inherits(model$expectation, "MxExpectationRAM")) {stop(paste("model '",model@name,"' does not use RAM expectation",sep=""))}
return(.mxStandardizeRAMhelper(model=model,SE=SE,ParamsCov=covParam))
}
#Handle multi-group model:
if(length(model@submodels)>0){
out <- NULL
if (inherits(model$expectation, "MxExpectationRAM")) {
out <- list(.mxStandardizeRAMhelper(model=model,SE=SE,ParamsCov=covParam,ignoreSubmodels=TRUE))
names(out)[1] <- model@name
}
if(!inde.subs.flag){
out <- c(out,lapply(
model@submodels[which(
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
)],
.mxStandardizeRAMhelper,SE=SE,ParamsCov=covParam))
if(length(out)==0){stop(paste("model '",model@name,"' does not use RAM expectation",sep=""))}
return(out)
}
else{
out1 <- lapply(
model@submodels[which(
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0})) &
!sapply(model@submodels,function(x){x@independent})
)],
.mxStandardizeRAMhelper,SE=SE,ParamsCov=covParam)
out <- c(out,out1,out2)
if(length(out)==0){stop(paste("model '",model@name,"' contains no submodels that use RAM expectation",sep=""))}
out <- out[names(model@submodels[which(
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
)])]
return(out)
}}}
#Alias using proper camel-casing:
mxStandardizeRAMPaths <- mxStandardizeRAMpaths
mxBootstrapStdizeRAMpaths <- function(model, bq=c(.25,.75), method=c('bcbci','quantile'), returnRaw=FALSE){
bq <- c(min(bq),max(bq))
if(!is(model, "MxModel")) {
stop("'model' argument must be a MxModel object")
}
if(!length(model@expectation) || !inherits(model@expectation, "MxExpectationRAM")) {
msg <- paste(
"MxModel ",omxQuotes(model@name),
" does not use RAM expectation\n(to use mxBootstrapStdizeRAMpaths() on a RAM submodel, run the function directly on that submodel",sep="")
stop(msg)
}
method <- match.arg(method)
realstdpaths <- .mxStandardizeRAMhelper(model=model,SE=FALSE,ParamsCov=NULL,inde.subs.flag=FALSE,ignoreSubmodels=TRUE)
rawParams <- as.matrix(omxGetBootstrapReplications(model))
#The tricky thing is that the output length of mxStandardizeRAMpaths() is not guaranteed to be the same for every replication...
outputlist <- vector("list",nrow(rawParams))
conformableFlag <- TRUE
for(i in 1:nrow(rawParams)){
modelcurr <- omxSetParameters(model,labels=colnames(rawParams),values=rawParams[i,])
stdpaths <- .mxStandardizeRAMhelper(model=modelcurr,SE=FALSE,ParamsCov=NULL,inde.subs.flag=FALSE,ignoreSubmodels=TRUE)
outputlist[[i]] <- stdpaths$Std.Value
names(outputlist[[i]]) <- stdpaths$name
if(conformableFlag && (nrow(stdpaths)!=nrow(realstdpaths) || !all(stdpaths$name==realstdpaths$name)) ){
conformableFlag <- FALSE
}
}
if( !conformableFlag ){
if(returnRaw){
warning("names of nonzero paths varied among bootstrap replications; returning raw list of standardized paths")
return(outputlist)
}
else{stop("names of nonzero paths varied among bootstrap replications, and argument 'returnRaw' is FALSE")}
}
else{
outmtx <- matrix(NA_real_,nrow=nrow(rawParams),ncol=nrow(realstdpaths))
colnames(outmtx) <- realstdpaths$name
for(i in 1:nrow(rawParams)){
outmtx[i,] <- as.vector(outputlist[[i]])
}
if(returnRaw){return(outmtx)}
}
out <- data.frame(realstdpaths$name,realstdpaths$label,realstdpaths$matrix,realstdpaths$row,realstdpaths$col,
realstdpaths$Std.Value,apply(outmtx,2,sd),numeric(length(realstdpaths$name)),numeric(length(realstdpaths$name)))
colnames(out) <- c("name","label","matrix","row","col","Std.Value","Boot.SE",
sprintf("%.1f%%", round(100*min(bq), 1)),sprintf("%.1f%%", round(100*max(bq), 1)))
if(method=="quantile"){
out[,8] <- as.vector(apply(outmtx,2,quantile,probs=min(bq),type=1))
out[,9] <- as.vector(apply(outmtx,2,quantile,probs=max(bq),type=1))
}
else if(method=="bcbci"){
zcrit <- qnorm(bq)
for(i in 1:nrow(realstdpaths)){
z0 <- qnorm(mean(outmtx[,i] <= realstdpaths$Std.Value[i]))
for (qx in 1:2){
phi <- pnorm(2*z0 + zcrit[qx])
out[i,7+qx] <- quantile(outmtx[,i], probs=phi, type=1)
}
}
}
else{warning("unrecognized value provided for argument 'method'")}
rownames(out) <- NULL
return(out)
}
coef.MxModel <- function(object, ...) omxGetParameters(object)
Try the OpenMx package in your browser
Any scripts or data that you put into this service are public.
OpenMx documentation built on June 8, 2025, 9:33 p.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 omxRMSEA fitStatistics catFitStatistics computeFitStatistics computeFValue numberObservations fitfunctionNumberObservations observedStatistics observedStatisticsHelper calculateConstraintsHelper calculateConstraints
Documented in omxRMSEA
#
# Copyright 2007-2021 by the individuals mentioned in the source code history
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
calculateConstraints <- function(model, flatModel) {
constraints <- flatModel@constraints
retval <- c()
if (length(constraints) > 0) {
retval <- sapply(constraints, calculateConstraintsHelper, model)
}
return(retval)
}
calculateConstraintsHelper <- function(constraint, model) {
if (constraint@relation == "==") {
leftHandSide <- constraint@formula[[2]]
value <- eval(substitute(mxEval(x, model, compute=TRUE),
list(x = leftHandSide)))
value <- as.matrix(value)
return(nrow(value) * ncol(value))
} else {
return(0)
}
}
observedStatisticsHelper <- function(model, expectation, datalist, historySet) {
if ('numStats' %in% slotNames(expectation)) {
if (length(expectation@numStats) > 0 && !is.na(expectation@numStats)) {
return(list(expectation@numStats, historySet))
}
}
if (length(expectation@data)==0 || is.na(expectation@data) || !.hasSlot(expectation, 'dims')) {
return(list(0, historySet))
}
if (is.numeric(expectation@data)) {
data <- datalist[[expectation@data + 1]]
} else {
data <- model[[expectation@data]]
}
if (is(data, "MxDataDynamic")) {
# need to revisit TODO
return(list(0, historySet))
}
obsStats <- data@observedStats
if (data@type == 'cov') {
if (data@name %in% historySet) {
return (list(0, historySet))
}
n <- nrow(data@observed)
dof <- n * (n + 1) / 2
if (!single.na(data@means)) {
dof <- dof + length(data@means)
}
historySet <- append(data, historySet)
} else if (data@type == 'cor') {
if (data@name %in% historySet) {
return (list(0, historySet))
}
n <- nrow(data@observed)
dof <- n * (n - 1) / 2
if (!single.na(data@means)) {
dof <- dof + length(data@means)
}
historySet <- append(data, historySet)
} else if (is(expectation, "MxExpectationBA81")) { # refactor TODO
if (!is.na(expectation@weightColumn) || !is.na(data@weight)) {
dof <- nrow(data@observed) - 1
} else {
dof <- nrow(rpf::compressDataFrame(data@observed)) - 1
}
historySet <- append(data, historySet)
} else if (!is.null(obsStats[['cov']])) {
if (data@name %in% historySet) {
return (list(0, historySet))
}
numThresh <- sum(!is.na(obsStats[['thresholds']]))
numMeans <- sum(!is.na(obsStats[['means']]))
numSlope <- sum(!is.na(obsStats[['slope']]))
n <- nrow(obsStats[['cov']])
dof <- numThresh + numSlope + numMeans + n*(n+1)/2 - 2*sum(ncol(obsStats[['thresholds']]))
historySet <- append(data, historySet)
} else {
# Incorporate row frequency and weight information? TODO
dof <- 0
observed <- data@observed
for (i in 1:ncol(observed)) {
colname <- colnames(observed)[[i]]
fullname <- paste(data@name, colname, sep='.')
if ((colname %in% expectation@dims) && !(fullname %in% historySet)) {
dof <- dof + sum(!is.na(observed[,i]))
historySet <- append(fullname, historySet)
}
}
}
return(list(dof, historySet))
}
observedStatistics <- function(model, flatModel, constraintOffset) {
datalist <- flatModel@datasets
labelsData <- imxGenerateLabels(model)
expectations <- convertExpectationFunctions(flatModel, model, labelsData, new("MxDirectedGraph"))
retval <- constraintOffset
if (length(expectations) > 0) {
historySet <- character()
for(i in 1:length(expectations)) {
result <- observedStatisticsHelper(model, expectations[[i]], datalist, historySet)
retval <- retval + result[[1]]
historySet <- result[[2]]
}
}
return(retval)
}
fitfunctionNumberObservations <- function(fitfunction) {
if (("numObsAdjust" %in% slotNames(fitfunction)) && !single.na(fitfunction@numObsAdjust)) {
return(fitfunction@numObsAdjust)
} else {
return(0)
}
}
numberObservations <- function(datalist, fitfunctions) {
datalist <- Filter(function(x) !is(x,"MxDataDynamic"), datalist)
dataObservations <- sapply(datalist, slot, name = "numObs")
fitfunctionObservations <- sapply(fitfunctions, fitfunctionNumberObservations)
return(sum(as.numeric(dataObservations), as.numeric(fitfunctionObservations)))
}
computeFValue <- function(datalist, likelihood, chi) {
if(length(datalist) == 0) return(NA)
datalist <- Filter(function(x) !is(x,"MxDataDynamic"), datalist)
if(all(sapply(datalist, function(x)
{length(x@observedStats) > 0 || is(x,"MxDataLegacyWLS") }))) return(chi)
if(all(sapply(datalist, function(x)
{x@type == 'raw'}))) return(likelihood)
if(all(sapply(datalist, function(x)
{x@type == 'cov'}))) return(chi)
if(all(sapply(datalist, function(x)
{x@type == 'cor'}))) return(chi)
return(NA)
}
computeFitStatistics <- function(likelihood, DoF, chi, chiDoF, numObs,
independence, indDoF, saturated=0, satDoF=0) {
if (!is.na(independence) && !is.na(likelihood) && independence < likelihood) {
warning(paste("Your model may be mis-specified (and fit worse than an independence model),",
"or you may be using the wrong independence model, see ?mxRefModels"))
}
CFI <- (independence - indDoF - likelihood + DoF)/(independence - indDoF - saturated + satDoF)
TLI <- NA
rmseaSquared <- 0
RMSEA <- NA
RMSEACI <- c(lower=NA, upper=NA)
RMSEANull <- 0.05
RMSEAClose <- NA
if (!is.na(chiDoF) && chiDoF > 0) {
TLI <- ((independence-saturated)/(indDoF-satDoF) - (chi)/(DoF-satDoF))/((independence-saturated)/(indDoF-satDoF) - 1)
# Here we use N in the denominator as given in the original
# RMSEA paper. The difference between N and N-1 is negligible
# for sample sizes over 30. RMSEA should not be taken seriously
# such small samples anyway.
rmseaSquared <- (chi / (chiDoF) - 1) / numObs
if (length(rmseaSquared) == 0 || is.na(rmseaSquared) ||
is.nan(rmseaSquared)) {
# || (rmseaSquared < 0)) { # changed so 'rmseaSquared < 0' yields zero with comment
RMSEA <- NA
} else {
RMSEA <- ifelse(rmseaSquared < 0, 0.0, sqrt(rmseaSquared))
ci <- try(rmseaConfidenceIntervalHelper(chi, chiDoF, numObs, .025, .975))
if (!inherits(ci, "try-error")) RMSEACI <- ci
RMSEAClose <- rmseaPCloseHelper(chi, chiDoF, numObs, null=RMSEANull)
}
}
list(CFI=CFI, TLI=TLI, RMSEA=RMSEA, RMSEASquared=rmseaSquared, RMSEACI=RMSEACI, RMSEAClose=RMSEAClose, RMSEANull=RMSEANull)
}
catFitStatistics <- function(x) {
cat("CFI:", x$CFI, '\n')
cat("TLI:", x$TLI, ' (also known as NNFI)', '\n')
if (length(x$RMSEASquared) == 1 && !is.na(x$RMSEASquared) && x$RMSEASquared < 0.0) {
cat("RMSEA: ", x$RMSEA, " *(Non-centrality parameter is negative)", " [95% CI (", x$RMSEACI[1], ", ", x$RMSEACI[2], ")]", '\n', sep="")
} else {
cat("RMSEA: ", x$RMSEA, " [95% CI (", x$RMSEACI[1], ", ", x$RMSEACI[2], ")]", '\n', sep="")
}
cat("Prob(RMSEA <= ", x$RMSEANull, "): ", x$RMSEAClose, '\n', sep='')
}
fitStatistics <- function(model, retval) {
datalist <- generateDataList(model)
likelihood <- retval[['Minus2LogLikelihood']]
saturated <- retval[['SaturatedLikelihood']]
independence <- retval[['IndependenceLikelihood']]
if (is.null(model@output$fitUnits)) return(retval)
if (model@output$fitUnits=="-2lnL" && is.null(model@output$chi)) {
chi <- likelihood - saturated
} else {chi <- model@output$chi}
if (is.null(chi)) return(retval)
DoF <- retval$degreesOfFreedom
satDoF <- retval$saturatedDoF
indDoF <- retval$independenceDoF
nParam <- dim(retval$parameters)[1]
Fvalue <- computeFValue(datalist, likelihood, chi)
if (model@output$fitUnits=="-2lnL" && is.null(model@output$chiDoF)) {
chiDoF <- DoF - satDoF # DoF = obsStat-model.ep; satDoF = obsStat-sat.ep; So sat.ep-model.ep == DoF-satDoF
} else {chiDoF <- model@output$chiDoF}
retval[['ChiDoF']] <- chiDoF
retval[['Chi']] <- chi
retval[['p']] <- suppressWarnings(ifelse(chiDoF==0,1.0,pchisq(chi, chiDoF, lower.tail = FALSE)))
retval[['AIC.Mx']] <- Fvalue - 2 * DoF
retval[['BIC.Mx']] <- (Fvalue - DoF * log(retval[['numObs']]))
AIC.p <- Fvalue + 2 * nParam
BIC.p <- (Fvalue + nParam * log(retval[['numObs']]))
sBIC <- (Fvalue + nParam * log((retval[['numObs']]+2)/24))
AICc <- Fvalue + 2*nParam + (2*nParam*(nParam+1))/(retval[['numObs']]-nParam-1)
retval[['satDoF']] <- satDoF
retval[['indDoF']] <- indDoF
IC <- matrix(NA, nrow=2, ncol=3, dimnames=list(c("AIC:", "BIC:"), c('df', 'par', 'sample')))
IC[,'df'] <- c(retval$AIC.Mx, retval$BIC.Mx)
IC[,'par'] <- c(AIC.p, BIC.p)
IC['BIC:','sample'] <- sBIC
IC['AIC:','sample'] <- AICc
if(length(model@output) && (is.null(model@output$fitUnits) || model@output$fitUnits=="-2lnL")){
retval[['informationCriteria']] <- IC
}
retval$fitUnits <- model@output$fitUnits
retval$fit <- if(retval$fitUnits == '-2lnL'){ retval[['Minus2LogLikelihood']] } else if(retval$fitUnits %in% c("r'Wr", "r'wr")) {retval[['Chi']]}
fi <- computeFitStatistics(likelihood, DoF, chi, chiDoF,
retval[['numObs']], independence, indDoF, saturated, satDoF)
for (k in names(fi)) retval[[k]] <- fi[[k]]
return(retval)
}
omxRMSEA <- function(model, lower=.025, upper=.975, null=.05, ...){
smod <- summary(model, ...)
x2 <- smod$Chi
df <- smod$ChiDoF
N <- smod$numObs
rmsea <- smod$RMSEA
ci <- rmseaConfidenceIntervalHelper(chi.squared=x2, df=df, N=N, lower=lower, upper=upper)
pn <- rmseaPCloseHelper(x2, df, N, null)
return(c(ci[1], est.rmsea=rmsea, ci[2], null=null, `Prob(x <= null)`=pn))
}
rmseaConfidenceIntervalHelper <- function(chi.squared, df, N, lower, upper){
# Lower confidence interval
if( pchisq(chi.squared, df=df, ncp=0) >= upper){ #sic
lower.lam <- uniroot(f=pChiSqFun, interval=c(1e-10, 1e4), val=chi.squared,
degf=df, goal=upper, extendInt="upX", maxiter=100L)$root
# solve pchisq(ch, df=df, ncp=x) == upper for x
} else{
lower.lam <- 0
}
# Upper confidence interval
if( pchisq(chi.squared, df=df, ncp=0) >= lower){ #sic
upper.lam <- uniroot(f=pChiSqFun, interval=c(1e-10, 1e4), val=chi.squared,
degf=df, goal=lower, extendInt="upX", maxiter=100L)$root
# solve pchisq(ch, df=df, ncp=x) == lower for x
} else{
upper.lam <- 0
}
lower.rmsea <- sqrt(lower.lam/(N*df))
upper.rmsea <- sqrt(upper.lam/(N*df))
return(c(lower=lower.rmsea, upper=upper.rmsea))
}
pChiSqFun <- function(x, val, degf, goal){
goal - pchisq(val, degf, ncp=x)
}
rmseaPCloseHelper <- function(x2, df, N, null){
1-pchisq(x2, df=df, ncp=N*df*(null)^2)
}
parameterList <- function(model, flatModel) {
parameterListHelper(model, flatModel, FALSE)
}
parameterListHelper <- function(model, flatModel, withModelName) {
ptable <- data.frame()
if(length(model@output) == 0) { return(ptable) }
estimates <- model@output$estimate
errorEstimates <- rep.int(as.numeric(NA), length(estimates))
if (!is.null(model@output$standardErrors)) {
se <- model@output$standardErrors
errorEstimates[match(rownames(se), names(estimates))] <- se
}
matrices <- generateMatrixList(flatModel)
parameters <- flatModel@parameters
if (length(estimates) > 0) {
matrixNames <- names(matrices)
for(i in 1:length(estimates)) {
mLocation <- parameters[[i]][[5]][[1]] + 1
mRow <- parameters[[i]][[5]][[2]] + 1
mCol <- parameters[[i]][[5]][[3]] + 1
mat <- matrices[[mLocation]][[1]]
dn <- dimnames(mat)
if (!is.null(dn[[1]])) mRow <- dn[[1]][mRow]
if (!is.null(dn[[2]])) mCol <- dn[[2]][mCol]
lbound <- parameters[[i]][[1]]
ubound <- parameters[[i]][[2]]
if (withModelName) {
ptable[i, 'model'] <- model@name
}
if (!is.null(names(estimates))) {
ptable[i, 'name'] <- names(estimates)[[i]]
} else {
ptable[i, 'name'] <- paste("p", i, sep="")
}
ptable[i, 'matrix'] <- simplifyName(matrixNames[[mLocation]], model@name)
ptable[i, 'row'] <- mRow
ptable[i, 'col'] <- mCol
ptable[i, 'Estimate'] <- estimates[[i]]
ptable[i, 'Std.Error'] <- errorEstimates[[i]]
ptable[i, 'lbound'] <- lbound
ptable[i, 'ubound'] <- ubound
}
}
return(ptable)
}
computeOptimizationStatistics <- function(model, flatModel, numStats, saturatedDoF, independenceDoF, retval) {
datalist <- flatModel@datasets
labelsData <- imxGenerateLabels(model)
expectations <- convertExpectationFunctions(flatModel, model, labelsData, new("MxDirectedGraph"))
# get estimated parameters
estimates <- model@output$estimate
# should saturated/independence models include means?
if(length(datalist)==1){
type <- datalist[[1]]@type
means <- datalist[[1]]@means
# if there's raw data, then use means in saturated/independence models
if(type=="raw"){
useMeans <- TRUE
} else {
# if there's not raw data, only use means if they're present
if((dim(means)[2]==1)&is.na(means[1,1])){
useMeans <- FALSE
} else{
useMeans <- TRUE
}
}
# number of variables
if(datalist[[1]]@type != 'raw'){
if(datalist[[1]]@type == 'none'){
nvar <- nrow(datalist[[1]]@observedStats$cov)
}
else{
nvar <- dim(datalist[[1]]@observed)[2]
}
} else if( length(expectations) == 1 ) {
nvar <- length(expectations[[1]]@dims)
} else {
nvar <- 0
}
# if there are multiple or zero datalists, then do nothing
} else {
useMeans <- NA
nvar <- 0
}
# how many thresholds does each variable have (needed for saturated and independence DoF calculation)
obj <- expectations
# grab the thresholdLevels object and expected means; punt if there is more than one expectation
if (length(obj)==1){
if ("thresholdLevels" %in% slotNames(obj[[1]])){
thresholdLevels <- obj[[1]]@thresholdLevels
if (length(thresholdLevels)==0){thresholdLevels <- rep(NA, nvar)}
} else {
thresholdLevels <- rep(NA, nvar)
}
} else {
thresholdLevels <- NULL
}
# number of continuous variables, provided there is just one expectation
if (!is.null(thresholdLevels)){
continuous <- sum(is.na(thresholdLevels))
} else{
continuous <- NA
}
# number of thresholds in the model
if (!is.null(thresholdLevels)){
thresh <- sum(thresholdLevels, na.rm=TRUE)
} else{
thresh <- NA
}
# constraints, parameters, model degrees of freedom
retval[['constraints']] <- calculateConstraints(model, flatModel)
retval[['estimatedParameters']] <- nrow(retval$parameters)
if(any(sapply(obj,function(x){is(x,"MxExpectationGREML")}))){
retval[['estimatedParameters']] <- retval[['estimatedParameters']] +
# When REML=FALSE and length(yhat)>0, numFixEff will be zero:
sum(sapply(obj,imxExtractSlot,name="numFixEff"))
}
if (is.null(numStats)) {
retval[['observedStatistics']] <- observedStatistics(model, flatModel, sum(retval$constraints))
} else {
retval[['observedStatistics']] <- numStats
}
retval[['degreesOfFreedom']] <- retval$observedStatistics - retval$estimatedParameters
# calculate or populate saturated degrees of freedom
if(is.null(saturatedDoF)) {
retval[['saturatedDoF']] <- retval$observedStatistics - (nvar * (nvar-1) / 2 + continuous*(1+useMeans) + thresh)
} else {
retval[['saturatedDoF']] <- saturatedDoF
}
#The "saturated model" has no sensible definiton with GREML expectation:
if(any(sapply(obj,function(x){is(x,"MxExpectationGREML")}))){
retval[['saturatedDoF']] <- NA
}
# calculate or populate independence degrees of freedom
if(is.null(independenceDoF)) {
if(!any(sapply(obj,function(x){is(x,"MxExpectationGREML")}))){
# indDoF = 1 df per continuous variable variance + 1 df per continuous mean + 1 df per threshold
retval[['independenceDoF']] <- retval$observedStatistics - (continuous*(1+useMeans) + thresh)
} else{
#TODO: the GREML expectation doesn't currently have a way to know how many phenotypes there are in every case.
#For now, leave the GREML independence model undefined
# #With GREML expectation, the independence model has a variance for each phenotype, and the same fixed effects as the fitted model:
# retval[['independenceDoF']] <-
# retval$observedStatistics - sum(sapply(obj,function(x){length(x@yvars)})) - sum(sapply(obj,imxExtractSlot,name="numFixEff"))
retval[['independenceDoF']] <- NA
}
} else {
retval[['independenceDoF']] <- independenceDoF
}
# set NULLs to NAs
if (is.null(retval$saturatedDoF)) {
retval$SaturatedDoF <- NA
}
if (is.null(retval$independenceDoF)) {
retval$IndependenceDoF <- NA
}
retval[['saturatedParameters']] <- retval[['observedStatistics']] - retval[['saturatedDoF']]
retval[['independenceParameters']] <- retval[['observedStatistics']] - retval[['independenceDoF']]
# calculate fit statistics
retval <- fitStatistics(model, retval)
return(retval)
}
print.summary.mxmodel <- function(x,...) {
cat("Summary of", x$modelName, '\n', '\n')
if(x$verbose==TRUE){
if (length(x$compute)) {
cat("compute plan:\n")
print(x$compute)
cat("\n")
}
if (length(x$dataSummary) > 0) {
cat("data:\n")
print(x$dataSummary)
}
}
if (!is.null(x$npsolMessage)) {
cat(x$npsolMessage,'\n','\n')
if ((x$statusCode == "not convex/red" || x$statusCode == "nonzero gradient/red") &&
(is.na(x$maxRelativeOrdinalError) || x$maxRelativeOrdinalError > 0)) {
cat(paste("Your ordinal model may converge if you reduce mvnRelEps\n",
" try: model <- mxOption(model, 'mvnRelEps', mxOption(model, 'mvnRelEps')/5)\n\n"))
}
}
params <- x$parameters
if (!is.null(params) && nrow(params)) {
cat("free parameters:\n")
params$lbound <- mapply(highlightProblem, params$lbound, params$lboundMet)
params$ubound <- mapply(highlightProblem, params$ubound, params$uboundMet)
params$lbound[is.na(params$lbound)] <- ""
params$ubound[is.na(params$ubound)] <- ""
params$lboundMet <- NULL
params$uboundMet <- NULL
if (!is.null(x$bootstrapSE) && length(x$bootstrapSE) == nrow(params)) {
params[['Std.Error']] <- x$bootstrapSE
}
if (!x$verbose) {
if (all(is.na(params[['Std.Error']]))) {
params[['Std.Error']] <- NULL
}
if (all(params$lbound == "" & params$ubound == "")) {
params[['lbound']] <- NULL
params[['ubound']] <- NULL
}
}
if (!is.null(params[['Std.Error']]) && !is.null(x[['seSuspect']])) {
seCol <- match('Std.Error', colnames(params))
before <- params[1:seCol]
stars <- mapply(highlightProblem, rep('',length(x[['seSuspect']])), x[['seSuspect']])
fullStars <- rep('', nrow(params))
fullStars[match(names(x[['seSuspect']]), params$name)] <- stars
if (length(params) > seCol) {
params <- cbind(before, 'A'=fullStars, params[(seCol+1):length(params)])
} else {
params <- cbind(before, 'A'=fullStars)
}
}
if (!is.null(x$bootstrapQuantile) && nrow(x$bootstrapQuantile) == nrow(params)) {
bq <- x$bootstrapQuantile
params <- cbind(params, bq)
}
cmap <- 1:ncol(params)
isBound <- colnames(params) %in% paste0(c('l','u'),'bound')
params <- params[,c(cmap[!isBound], cmap[isBound]),drop=FALSE]
print(params)
cat('\n')
}
if (!is.null(x$CI) && length(x$CI) > 0) {
cat("confidence intervals:\n")
print(x$CI)
if (any(is.na(x$CI[,c('lbound','ubound')])) && !(x$verbose)) {
cat(" To investigate missing CIs, run summary() again, with verbose=T, to see CI details.", '\n')
}
cat('\n')
} else if (x$CI.Requested) {
cat("To obtain confidence intervals re-run with intervals=TRUE\n\n")
}
if(x$verbose && length(x$CIdetail)){
cat("CI details:\n")
print(x$CIdetail)
cat("\n")
}
if(length(x$GREMLfixeff)>0 && any(sapply(x$GREMLfixeff,length)>0)){
cat("regression coefficients:\n")
print(x$GREMLfixeff)
cat("\n")
}
cat('Model Statistics:', '\n')
EP <- matrix(
c(x$estimatedParameters, x$saturatedParameters, x$independenceParameters,
x$degreesOfFreedom, x$saturatedDoF, x$independenceDoF,
x$Minus2LogLikelihood, x$SaturatedLikelihood, x$IndependenceLikelihood),
nrow=3, ncol=3,
dimnames=list(
c(' Model:', ' Saturated:', 'Independence:'),
c(' | Parameters', ' | Degrees of Freedom', paste0(' | Fit (', x$fitUnits, ' units)'))
)
)
print(EP)
cat('Number of observations/statistics: ', x$numObs, "/", x$observedStatistics, '\n\n', sep="")
constraints <- x$constraints
if(length(constraints) > 0) {
for(i in 1:length(constraints)) {
name <- names(constraints)[[i]]
if (constraints[[i]] == 1) plural <- ''
else plural <- 's'
cat("Constraint", omxQuotes(simplifyName(name, x$modelName)), "contributes",
constraints[[i]], paste("observed statistic", plural, '.', sep=''), "\n")
if(i==length(constraints)){cat("\n")}
}
}
if (!is.null(x$infoDefinite) && !is.na(x$infoDefinite)) {
if (!x$infoDefinite) {
cat("\n** Information matrix is not positive definite (not at a candidate optimum).\n Be suspicious of these results. At minimum, do not trust the standard errors.\n\n")
} else if(x$verbose==TRUE) {
cat("condition number of the information matrix: ", x$conditionNumber, "\n")
}
}
if (x$verbose==TRUE && !is.null(x$maxAbsGradient)) {
cat("maximum absolute gradient: ", x$maxAbsGradient, " (",names(x$maxAbsGradient),")\n")
}
#
# Chi-square goodness of fit test
if(x$verbose==TRUE || (!is.null(x$Chi) && !is.na(x$Chi))) {
chival <- x$Chi
if(is.na(x$SaturatedLikelihood)){
chidof <- NA
} else {
chidof <- x$ChiDoF
}
chipee <- x$p
cat("chi-square: ", "\U03C7\U00B2 ( df=", chidof, " ) = ", chival, ", p = ", chipee, '\n', sep="")
}
#
if(length(x$informationCriteria)){
cat("Information Criteria: \n")
IC <- x$informationCriteria
colnames(IC) <- c(" | df Penalty", " | Parameters Penalty", " | Sample-Size Adjusted")
print(IC)
}
#
# Absolute fit indices
if(x$verbose==TRUE || any(!is.na(c(x$CFI, x$TLI, x$RMSEA)))){
catFitStatistics(x)
}
if(any(is.na(c(x$CFI, x$TLI, x$RMSEA)))){
cat("To get additional fit indices, see help(mxRefModels)\n")
}
#
# Timing information
cat("timestamp:", format(x$timestamp), '\n')
if(x$verbose==TRUE){
cat("frontend time:", format(x$frontendTime), '\n')
cat("backend time:", format(x$backendTime), '\n')
cat("independent submodels time:", format(x$independentTime), '\n')
cat("cpu time:", format(x$cpuTime), '\n')
}
cat("Wall clock time:", format(x$wallTime), "\n")
if (x$verbose==FALSE && !is.null(x$optimizerEngine)) {
cat("optimizer: ", x$optimizerEngine, '\n')
}
cat("OpenMx version number:", format(x$mxVersion), '\n')
cat("Need help? See help(mxSummary)", '\n')
cat('\n')
if (!x$wasRun) {
message("WARNING: This model has not been run yet. Tip: Use\n model = mxRun(model)\nto estimate a model.")
} else if (x$stale) {
message("WARNING: This model was modified since it was run. Summary information may be out-of-date.")
}
}
setLikelihoods <- function(model, saturatedLikelihood, independenceLikelihood, retval) {
# populate saturated -2 log likelihood
if(is.null(saturatedLikelihood)) {
retval$SaturatedLikelihood <- attr(model@fitfunction$result, "SaturatedLikelihood")
} else {
retval$SaturatedLikelihood <- saturatedLikelihood
}
# populate independence -2 log likelihood
if(is.null(independenceLikelihood)) {
retval$IndependenceLikelihood <- attr(model@fitfunction$result, "IndependenceLikelihood")
} else {
retval$IndependenceLikelihood <- independenceLikelihood
}
# populate model -2 log likelihood
retval$Minus2LogLikelihood <- model@output$Minus2LogLikelihood
# set NULLs to NAs
if (is.null(retval$SaturatedLikelihood)) {
retval$SaturatedLikelihood <- NA
}
if (is.null(retval$Minus2LogLikelihood)) {
retval$Minus2LogLikelihood <- NA
}
if (is.null(retval$IndependenceLikelihood)) {
retval$IndependenceLikelihood <- NA
}
return(retval)
}
setNumberObservations <- function(numObs, datalist, fitfunctions, retval) {
if(is.null(numObs)) {
retval$numObs <- numberObservations(datalist, fitfunctions)
} else {
retval$numObs <- numObs
}
return(retval)
}
generateDataSummary <- function(model) {
datalist <- generateDataList(model)
retval <- lapply(datalist, summarize)
return(retval)
}
##' imxEvalByName
##'
##' This is an internal function exported for those people who know
##' what they are doing.
##'
##' @param name name
##' @param model model
##' @param compute compute
##' @param show show
imxEvalByName <- function(name, model, compute=FALSE, show=FALSE) {
if ((length(name) != 1) || typeof(name) != "character") {
stop("'name' argument must be a character argument")
}
if (!is(model, "MxModel")) {
stop("'model' argument must be a MxModel object")
}
if (hasSquareBrackets(name)) {
components <- splitSubstitution(name)
eval(substitute(mxEval(x[y,z], model, compute, show),
list(x = as.name(components[[1]]),
y = parse(text=components[[2]])[[1]],
z = parse(text=components[[3]])[[1]])))
} else {
eval(substitute(mxEval(x, model, compute, show),
list(x = as.name(name))))
}
}
boundsMet <- function(model, retval){
params <- retval$parameters
lbound <- params$lbound
ubound <- params$ubound
estimate <- params$Estimate
threshold <- model@options[["Feasibility tolerance"]]
if (is.null(threshold)){
threshold <- getOption("mxOptions")[["Feasibility tolerance"]]
}
threshold <- as.numeric(threshold)
lboundMet <- mapply(compareBounds, estimate, lbound, MoreArgs=list(threshold))
uboundMet <- mapply(compareBounds, estimate, ubound, MoreArgs=list(threshold))
params$lboundMet <- lboundMet
params$uboundMet <- uboundMet
retval$parameters <- params
return(retval)
}
compareBounds <- function(estimate, bound, threshold){
if (is.na(bound)){
return(FALSE)
}
absDelta <- abs(estimate - bound)
return (absDelta < threshold)
}
highlightProblem <- function(bound, boundMet){
if (boundMet){
if (is.numeric(bound)) bound <- round(bound,4)
return(paste(bound, "!", sep=""))
}
else {
return(bound)
}
}
parseLikelihoodArg <- function(input, arg) {
input <- input[[arg]]
if (is.null(input)) {
return(input)
} else if (is.numeric(input)) {
return(input)
} else if (is(input, "MxModel")) {
name <- input@name
if (is.null(input@fitfunction)) {
stop(paste(omxQuotes(name), "model passed",
"to summary function does not",
"have top-level fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
if (length(input@fitfunction@result) != 1) {
stop(paste(omxQuotes(name), "model passed to summary",
"function does not have a 1x1 matrix",
"result in fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
return(input@fitfunction@result[1,1])
} else if(is.list(input) && length(input)==2) {
stop(paste("List of length two (illegal argument) passed to", omxQuotes(arg),
"argument of summary function. You probably meant to use",
"the refModels argument instead."), call. = FALSE)
} else {
stop(paste("Illegal argument passed to", omxQuotes(arg),
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
parseDfArg <- function(input, arg) {
input <- input[[arg]]
if (is.null(input)) {
return(input)
} else if (is.numeric(input)) {
return(input)
} else if (is(input, "MxModel")) {
name <- input@name
if (is.null(input@fitfunction)) {
stop(paste(omxQuotes(name), "model passed",
"to summary function does not",
"have top-level fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
if (length(input@fitfunction@result) != 1) {
stop(paste(omxQuotes(name), "model passed to summary",
"function does not have a 1x1 matrix",
"result in fitfunction in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
return(summary(input)$degreesOfFreedom)
} else {
stop(paste("Illegal argument passed to", omxQuotes(arg),
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
refToLikelihood <- function(model) {
if (is(model, "MxModel")) {
if (!model@.wasRun) stop("Reference model must be run to obtain fit indices")
model$output$Minus2LogLikelihood
} else if (is.list(model)) {
model[[1]]
} else {
stop(paste("Illegal argument passed to refModels",
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
refToDof <- function(model) {
if (is(model, "MxModel")) {
if (!model@.wasRun) stop("Reference model must be run to obtain fit indices")
return(summary(model)$degreesOfFreedom)
} else if (is.list(model)) {
model[[2]]
} else {
stop(paste("Illegal argument passed to refModels",
"argument of summary function in",
deparse(width.cutoff = 400L, sys.call(-1))), call. = FALSE)
}
}
summarizeBootstrap <- function(mle, bootData, bq, summaryType) {
if (summaryType == 'quantile') {
t(apply(bootData, 2, quantile, probs=bq, type=1))
} else if (summaryType == 'bcbci') {
zcrit <- qnorm(bq)
out <- matrix(NA, nrow=length(mle), ncol=length(bq),
dimnames=list(names(mle),
sapply(bq, function(x) sprintf("%.1f%%", round(100*min(x), 1)))))
for(i in 1:length(mle)) {
z0 <- qnorm(mean(bootData[,i] <= mle[i]))
for (qx in 1:length(bq)) {
phi <- pnorm(2*z0 + zcrit[qx])
out[i,qx] <- quantile(bootData[,i], probs=phi, type=1)
}
}
out
} else {
warning(paste("boot.SummaryType =", omxQuotes(summaryType),
"is not recognized"))
}
}
summary.MxModel <- function(object, ..., verbose=FALSE) {
out <- try(summary.MxModel.Impl(object, ..., verbose=verbose), silent=TRUE)
if (is(out, "try-error")) {
if (!object@.wasRun) {
# Maybe should never happen?
stop("This model has not been run yet. Tip: Use\n model = mxRun(model)\nto estimate a model.")
} else if (object@.modifiedSinceRun) {
stop(paste("This model is in an inconsistent state.",
"Tip: Use\n model = mxRun(model)\nto re-estimate the model."))
} else {
stop(out)
}
}
out
}
summary.MxModel.Impl <- function(object, ..., verbose) {
model <- object
dotArguments <- list(...)
if (!is.null(dotArguments[["refModels"]])) {
refModels <- dotArguments[["refModels"]]
satModel <- refModels[['Saturated']]
indModel <- refModels[['Independence']]
saturatedLikelihood <- refToLikelihood(satModel)
saturatedDoF <- refToDof(satModel)
independenceLikelihood <- refToLikelihood(indModel)
independenceDoF <- refToDof(indModel)
} else {
saturatedLikelihood <- parseLikelihoodArg(dotArguments, "SaturatedLikelihood")
saturatedDoF <- parseDfArg(dotArguments, "SaturatedDoF")
independenceLikelihood <- parseLikelihoodArg(dotArguments, "IndependenceLikelihood")
independenceDoF <- parseDfArg(dotArguments, "IndependenceDoF")
}
numObs <- dotArguments$numObs
numStats <- dotArguments$numStats
namespace <- imxGenerateNamespace(model)
flatModel <- imxFlattenModel(model, namespace, TRUE)
flatModel <- constraintsToAlgebras(flatModel)
flatModel <- eliminateObjectiveFunctions(flatModel)
convertArguments <- imxCheckVariables(flatModel, namespace)
flatModel <- convertAlgebras(flatModel, convertArguments)
labelsData <- imxGenerateLabels(model)
flatModel <- expectationFunctionConvertEntities(flatModel, namespace, labelsData)
flatModel <- populateDefInitialValues(flatModel)
flatModel <- checkEvaluation(model, flatModel)
dependencies <- cycleDetection(flatModel)
dependencies <- transitiveClosure(flatModel, dependencies)
freeVarGroups <- buildFreeVarGroupList(flatModel)
flatModel <- generateParameterList(flatModel, dependencies, freeVarGroups)
retval <- list(wasRun=model@.wasRun, stale=model@.modifiedSinceRun)
retval$parameters <- parameterList(model, flatModel)
if (!is.null(model@compute$steps[['ND']]) && model@compute$steps[['ND']]$checkGradient &&
!is.null(model@compute$steps[['ND']]$output$gradient)) {
gdetail <- model@compute$steps[['ND']]$output$gradient
retval$seSuspect <- !gdetail[,'symmetric']
names(retval$seSuspect) <- rownames(gdetail)
}
if (is(model@compute, "MxComputeBootstrap")) {
bq <- c(.25,.75)
if (!is.null(dotArguments[["boot.quantile"]])) {
bq <- sort(as.numeric(dotArguments[["boot.quantile"]]))
}
summaryType <- 'bcbci'
if (!is.null(dotArguments[["boot.SummaryType"]])) {
summaryType <- dotArguments[["boot.SummaryType"]]
}
cb <- model@compute
if (!is.null(cb@output$raw) && is.na(cb@only) && cb@output$numParam == nrow(retval$parameters)) {
raw <- cb@output$raw
mask <- raw[,'statusCode'] %in% cb@OK
bootData <- raw[mask, 3:(nrow(retval$parameters)+2), drop=FALSE]
if (sum(mask) < .95*nrow(raw)) {
pct <- round(100*sum(mask) / nrow(raw))
warning(paste0("Only ",pct,"% of the bootstrap replications ",
"converged. Accuracy is much less than the ", nrow(raw),
" replications requested"), call.=FALSE)
}
if (sum(mask) >= 3) {
retval$bootstrapSE <- apply(bootData, 2, sd)
retval$bootstrapQuantile <-
summarizeBootstrap(retval$parameters[, 'Estimate'], bootData, bq, summaryType)
}
}
} else if (any(grep('^boot\\.', names(dotArguments)))) {
warning("No bootstrap data found. See ?mxBootstrap")
}
retval$GREMLfixeff <- GREMLFixEffList(model)
retval$infoDefinite <- model@output$infoDefinite
retval$conditionNumber <- model@output$conditionNumber
if (length(model@output$gradient)) {
agrad <- abs(model@output$gradient)
retval$maxAbsGradient <- agrad[ order(-agrad)[1] ]
}
retval <- boundsMet(model, retval)
retval <- setLikelihoods(model, saturatedLikelihood, independenceLikelihood, retval)
labelsData <- imxGenerateLabels(model)
fitfunctions <- convertFitFunctions(flatModel, model, labelsData, dependencies)
retval <- setNumberObservations(numObs, flatModel@datasets, fitfunctions, retval)
retval <- computeOptimizationStatistics(model, flatModel, numStats, saturatedDoF, independenceDoF, retval)
retval$dataSummary <- generateDataSummary(model)
retval$CI.Requested <- length(model@intervals) > 0
retval$CI <- as.data.frame(model@output$confidenceIntervals)
if (length(retval$CI) && nrow(retval$CI)) {
retval$CI <- cbind(retval$CI, note=apply(retval$CI, 1, function(ci) {
# This should probably take into account whether both bounds
# were requested and consider the optimizer codes also. TODO
if (any(is.na(ci)) || ci[1] == ci[3] || ci[1] >= ci[2] || ci[2] >= ci[3]) {
"!!!"
} else {
""
}
}))
}
retval$CIcodes <- model@output$confidenceIntervalCodes
statusCode <- model@output$status$code
if (!is.null(statusCode)) {
message <- optimizerMessages[[as.character(statusCode)]]
retval[['npsolMessage']] <- message
retval[['statusCode']] <- as.statusCode(statusCode)
retval[['maxRelativeOrdinalError']] <- model@output[['maxRelativeOrdinalError']]
}
if( .hasSlot(model,"compute") && length(model$compute$steps$CI) ){
retval$CIdetail <- model$compute$steps$CI$output$detail
}
retval$timestamp <- model@output$timestamp
retval$frontendTime <- model@output$frontendTime
retval$backendTime <- model@output$backendTime
retval$independentTime <- model@output$independentTime
retval$wallTime <- model@output$wallTime
retval$cpuTime <- model@output$cpuTime
retval$mxVersion <- model@output$mxVersion
retval$modelName <- model@name
plan <- model@compute
if (is(plan, "MxComputeSequence")) {
gd <- plan$steps[['GD']]
if (is(gd, "MxComputeGradientDescent")) {
retval$optimizerEngine <- gd$engine
}
}
retval$verbose <- verbose
class(retval) <- "summary.mxmodel"
return(retval)
}
assertModelRunAndFresh <- function(model) {
warnModelCreatedByOldVersion(model)
if (.hasSlot(model,".wasRun") && !model@.wasRun) stop("This model has not been run yet. Tip: Use\n model = mxRun(model)\nto estimate a model.")
if (.hasSlot(model,".modifiedSinceRun") && model@.modifiedSinceRun) {
msg <- paste("MxModel", omxQuotes(model@name), "was modified",
"since it was run.")
warning(msg)
}
}
assertModelFreshlyRun <- function(model) {
warnModelCreatedByOldVersion(model)
if (model@.wasRun && model@.modifiedSinceRun) {
msg <- paste("MxModel", omxQuotes(model@name), "was modified",
"since it was run.")
warning(msg)
}
}
logLik.MxModel <- function(object, ...) {
model <- object
moreModels <- list(...)
assertModelFreshlyRun(model)
ll <- NA
if (length(model@output) && !is.null(model@output$fit) &&
!is.null(model@output$fitUnits) ) {
if(model@output$fitUnits=="-2lnL"){
ll <- -0.5*model@output$fit
} else if(model@output$fitUnits %in% c("r'Wr", "r'wr")) {
ll <- -0.5*model@output$chi
}
#TODO: this doesn't count "implicit" free parameters that are "profiled out":
attr(ll, "df") <- length(model@output$estimate)
} else {
attr(ll,"df") <- NA
}
namespace <- imxGenerateNamespace(model)
flatModel <- imxFlattenModel(model, namespace, TRUE)
labelsData <- imxGenerateLabels(model)
fitfunctions <- convertFitFunctions(flatModel, model, labelsData, new("MxDirectedGraph"))
nobs <- numberObservations(flatModel@datasets, fitfunctions)
if (!is.na(nobs)) {
attr(ll,"nobs") <- nobs
}
class(ll) <- "logLik"
if (length(moreModels)) {
c(list(ll), lapply(moreModels, logLik.MxModel))
} else {
ll
}
}
.standardizeParams <- function(x=NULL, model, Apos, Spos, Mpos, give.matrices=FALSE){
if(is.null(x)){x <- omxGetParameters(model)}
param <- omxGetParameters(model)
paramNames <- names(param)
model <- omxSetParameters(model, values=x, labels=paramNames, free=TRUE)
model_A <- model[[model$expectation$A]] #<--the A matrix might not be named "A".
A <- list( model_A$values, model_A$result )
A <- A[[which.max(c( length(A[[1]]), length(A[[2]]) ))]]
model_S <- model[[model$expectation$S]] #<--Likewise for S
S <- list( model_S$values, model_S$result )
S <- S[[which.max(c( length(S[[1]]), length(S[[2]]) ))]]
M <- NULL
if(!is.null(Mpos)){
model_M <- model[[model$expectation$M]]
M <- list( model_M$values, model_M$result )
M <- M[[which.max(c( length(M[[1]]), length(M[[2]]) ))]]
}
else{model_M <- NULL}
I <- diag(1, nrow(A))
ImAInv <- solve(I-A)
SD <- sqrt(diag(ImAInv %*% S %*% t(ImAInv)))
SD <- diag(SD,nrow=length(SD)) #<--Needed if line immediately above is a scalar.
InvSD <- 1/diag(SD)
InvSD <- diag(InvSD,nrow=length(InvSD))
Az <- InvSD %*% A %*% SD
Sz <- InvSD %*% S %*% InvSD
if(!is.null(M)){
Mz <- M %*% InvSD
sparam <- c(Az[!is.na(Apos)],Sz[!is.na(Spos)],Mz[!is.na(Mpos)])
names(sparam) <- c(Apos[!is.na(Apos)],Spos[!is.na(Spos)],Mpos[!is.na(Mpos)])
if(!give.matrices){return(sparam)}
else{return(list(sparam=sparam,Az=Az,Sz=Sz,Mz=Mz))}
}
else{
sparam <- c(Az[!is.na(Apos)],Sz[!is.na(Spos)])
names(sparam) <- c(Apos[!is.na(Apos)],Spos[!is.na(Spos)])
if(!give.matrices){return(sparam)}
else{return(list(sparam=sparam,Az=Az,Sz=Sz))}
}
}
.mxStandardizeRAMhelper <- function(model,SE=FALSE,ParamsCov,inde.subs.flag=FALSE,ignoreSubmodels=FALSE){
#Recur the function for the appropriate submodels, if any:
if(length(model@submodels) && !ignoreSubmodels){
return(lapply(
model@submodels[which(
sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0})
)],
.mxStandardizeRAMhelper,SE=SE,ParamsCov=ParamsCov,inde.subs.flag=inde.subs.flag,ignoreSubmodels=FALSE))
}
#Get A and S:
model_A <- model[[model$expectation$A]] #<--Necessary because the A matrix might not be named "A".
A <- list( model_A$values, model_A$result )
A <- A[[which.max(c( length(A[[1]]), length(A[[2]]) ))]]
model_S <- model[[model$expectation$S]] #<--Likewise for S
S <- list( model_S$values, model_S$result )
S <- S[[which.max(c( length(S[[1]]), length(S[[2]]) ))]]
M <- NULL
if(!is.na(model$expectation$M)){
model_M <- model[[model$expectation$M]]
M <- list( model_M$values, model_M$result )
M <- M[[which.max(c( length(M[[1]]), length(M[[2]]) ))]]
}
#Find positions of nonzero paths:
Apos <- matrix(NA,nrow=nrow(A),ncol=ncol(A),dimnames=dimnames(A))
Spos <- matrix(NA,nrow=nrow(S),ncol=ncol(S),dimnames=dimnames(S))
A_need_pos <- which(A!=0,arr.ind=T)
S_need_pos <- which(S!=0,arr.ind=T)
S_need_pos <- subset(S_need_pos, S_need_pos[,1]>=S_need_pos[,2]) #<--Lower tri only
if(!is.null(M)){
Mpos <- matrix(NA,nrow=1,ncol=ncol(M),dimnames=dimnames(M))
M_need_pos <- which(M!=0)
}
else{
Mpos <- NULL
M_need_pos <- NULL
}
numelem <- nrow(A_need_pos)+nrow(S_need_pos)+length(M_need_pos)
#Create output object:
out <- data.frame(name=vector(mode="character",length=numelem),label=vector(mode="character",length=numelem),
matrix=vector(mode="character",length=numelem),
row=vector(mode="character",length=numelem),col=vector(mode="character",length=numelem),
Raw.Value=vector(mode="numeric",length=numelem),
Raw.SE=vector(mode="numeric",length=numelem),
Std.Value=vector(mode="numeric",length=numelem),
Std.SE=vector(mode="numeric",length=numelem),stringsAsFactors=FALSE)
out$label <- NA
j <- 1
#Create position strings where needed and begin to populate output:
if(nrow(A_need_pos)>0){
for(i in 1:nrow(A_need_pos)){
Apos[A_need_pos[i,1],A_need_pos[i,2]] <- out$name[j] <- paste(
model@name,".A[",A_need_pos[i,1],",",A_need_pos[i,2],"]",sep="")
if(!isAllNa(model_A$labels)){
out$label[j] <- model_A$labels[A_need_pos[i,1],A_need_pos[i,2]]
}
out$matrix[j] <- "A"
out$Raw.Value[j] <- A[A_need_pos[i,1],A_need_pos[i,2]]
out$row[j] <- ifelse(length(rownames(A))>0,rownames(A)[A_need_pos[i,1]],A_need_pos[i,1])
out$col[j] <- ifelse(length(colnames(A))>0,colnames(A)[A_need_pos[i,2]],A_need_pos[i,2])
j <- j+1
}
}
if(nrow(S_need_pos)>0){
for(i in 1:nrow(S_need_pos)){
Spos[S_need_pos[i,1],S_need_pos[i,2]] <- out$name[j] <- paste(
model@name,".S[",S_need_pos[i,1],",",S_need_pos[i,2],"]",sep="")
if(!isAllNa(model_S$labels)){
out$label[j] <- model_S$labels[S_need_pos[i,1],S_need_pos[i,2]]
}
out$matrix[j] <- "S"
out$Raw.Value[j] <- S[S_need_pos[i,1],S_need_pos[i,2]]
out$row[j] <- ifelse(length(rownames(S))>0,rownames(S)[S_need_pos[i,1]],S_need_pos[i,1])
out$col[j] <- ifelse(length(colnames(S))>0,colnames(S)[S_need_pos[i,2]],S_need_pos[i,2])
j <- j+1
}
}
if(length(M_need_pos)){
for(i in 1:length(M_need_pos)){
Mpos[1,M_need_pos[i]] <- out$name[j] <- paste(
model@name,".M[1,",M_need_pos[i],"]",sep="")
if(!is.null(model_M$labels) && !isAllNa(model_M$labels)){
out$label[j] <- model_M$labels[1,M_need_pos[i]]
}
out$matrix[j] <- "M"
out$Raw.Value[j] <- M[1,M_need_pos[i]]
out$row[j] <- 1
out$col[j] <- ifelse(length(colnames(M))>0,colnames(M)[M_need_pos[i]],M_need_pos[i])
j <- j+1
}
}
#Get standardized values:
freeparams <- omxGetParameters(model)
paramnames <- names(freeparams)
zout <- .standardizeParams(x=freeparams,model=model,Apos=Apos,Spos=Spos,Mpos=Mpos)
#Compute SEs, or assign them 'not requested' values, as the case may be:
if(SE){
if(!all(paramnames %in% rownames(ParamsCov))){
stop(paste(
"some free parameter labels do not appear in the dimnames of the parameter estimates' covariance matrix;\n",
"are you running mxStandardizeRAMpaths() on a dependent submodel instead of on its multigroup container model?\n",
"the missing parameter labels are:\n",
paste(paramnames[!(paramnames %in% rownames(ParamsCov))],collapse=", "),sep=""))
}
#From Mike Hunter's delta method example:
covParam <- ParamsCov[paramnames,paramnames,drop=FALSE]#<--submodel will usually not contain all free param.s
jacStand <- numDeriv::jacobian(func=.standardizeParams, x=freeparams, model=model, Apos=Apos, Spos=Spos, Mpos=Mpos)
covSparam <- jacStand %*% covParam %*% t(jacStand)
dimnames(covSparam) <- list(names(zout),names(zout))
SEs <- sqrt(diag(covSparam))
#SEs[diag(covSparam)<.Machine$double.eps] <- 0
}
else{SEs <- rep("not_requested",length(zout)); names(SEs) <- names(zout)}
#Add standardized values and SEs to output:
out$Std.Value <- zout
out$Std.SE <- SEs
#Pull in raw SEs if requested:
if(SE){
for(i in 1:numelem){
if( (out$name[i] %in% paramnames) |
(out$label[i] %in% paramnames) ){
tdiags <- covParam[ifelse(is.na(out$label[i]),out$name[i],out$label[i]),
ifelse(is.na(out$label[i]),out$name[i],out$label[i])]
if (length(tdiags) == 1) {
# For diag, R will return a square identity matrix of size given by the scalar
if(tdiags < 0 || is.na(tdiags)) {
warning("Some diagonal elements of the repeated-sampling covariance matrix of the point estimates are less than zero or NA.\nThat's weird. Raise an eyebrow at these standard errors.")
}
} else {
if(any(diag(tdiags) < 0) || any(is.na(tdiags))){
warning("Some diagonal elements of the repeated-sampling covariance matrix of the point estimates are less than zero or NA.\nThat's weird. Raise an eyebrow at these standard errors.")
}
}
out$Raw.SE[i] <- suppressWarnings(sqrt(tdiags))
}
else{
# # The elements of `out$name` should all reference scalars:
# tmp <- try(mxSE(x=out$name[i],model=model,forceName=TRUE,silent=TRUE))
# if(inherits(tmp,"try-error")){
# # Paths with the same label will have the same raw value and raw-value SE:
# tmp <- try(mxSE(x=out$label[i],model=model,forceName=TRUE,silent=TRUE))
# }
# if( !inherits(tmp,"try-error") && is.numeric(tmp) ){
# out$Raw.SE[i] <- tmp
# }
}
}
}
else{out$Raw.SE <- "not_requested"}
return(out)
}
mxStandardizeRAMpaths <- function(model, SE=FALSE, cov=NULL){
assertModelFreshlyRun(model)
if(imxHasDefinitionVariable(model)){
warning("'model' (or one of its submodels) contains definition variables; interpret results of mxStandardizeRAMpaths() cautiously")
}
#If SE=T,need to check for independent submodels because they will have their own Hessians;
#recur main function as appropriate:
inde.subs.flag <- FALSE
if(SE & length(model$submodels)>0){
RAM.subs <- (sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
inde.subs <- sapply(model@submodels,function(x){x@independent})==TRUE &
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
if(sum(inde.subs)>0){
out2 <- NULL
#if ALL submodels are either independent RAM models or non-RAM models:
if(all(RAM.subs==inde.subs)){
out <- lapply(model@submodels[which(inde.subs)],mxStandardizeRAMpaths,SE=T)
if(length(out)==0){stop(paste("model '",model@name,"' contains no submodels that use RAM expectation",sep=""))}
return(out)
}
else{out2 <- lapply(model@submodels[which(inde.subs)],mxStandardizeRAMpaths,SE=T)}
inde.subs.flag <- TRUE
}}
covParam <- NULL
if(SE){
#If user requests SEs and provided no covariance matrix, check to be sure SEs can and should be computed:
if(!length(cov)){
# if(length(model@constraints)>0){
# msg <- paste("standard errors will not be computed because model '",model@name,"' contains at least one mxConstraint",sep="")
# warning(msg)
# SE <- FALSE
# }
if(SE & length(model@output$vcov)==0){
if(!model@.wasRun){
msg <- paste("standard errors will not be computed because model '",model@name,"' has not yet been run, and no matrix was provided for argument 'cov'",sep="")
warning(msg)
SE <- FALSE
}
else{
warning("argument 'SE=TRUE' requires model to have a nonempty 'vcov' output slot, or a non-NULL value for argument 'cov'; continuing with 'SE' coerced to 'FALSE'")
SE <- FALSE
}}
libraries <- rownames(installed.packages())
pkgcheck <- ("numDeriv" %in% libraries)
if(SE & !pkgcheck){
warning("argument 'SE=TRUE' requires package 'numDeriv' to be installed; continuing with 'SE' coerced to 'FALSE'")
SE <- FALSE
}
if(SE){
covParam <- vcov(model)
}
}
#If user requests SEs and provided a covariance matrix:
else{
#Conceivably, the user could provide a sampling covariance matrix that IS valid in the presence of MxConstraints...
if(length(model@constraints)>0){
#msg <- paste("standard errors may be invalid because model '",model@name,"' contains at least one mxConstraint",sep="")
#warning(msg)
}
#Sanity checks on the value of argument 'cov':
if(!is.matrix(cov)){ #<--Is it a matrix?
cov <- try(as.matrix(cov),silent=T)
if("try-error" %in% class(cov) || !is.matrix(cov)){ #<--If its not a matrix, can it be coerced to one?
stop("non-NULL value to argument 'cov' must be (or be coercible to) a matrix")
}
}
if(nrow(cov)!=ncol(cov)){ #<--Is it square?
msg <- paste("non-NULL value to argument 'cov' must be a square matrix; it has ",nrow(cov)," rows and ",ncol(cov)," columns",sep="")
stop(msg)
}
#Do its row and column names match?:
if(!length(rownames(cov)) || !length(colnames(cov)) || any(rownames(cov) != colnames(cov))){
stop("non-NULL value to argument 'cov' must have matching and complete rownames and colnames")
}
paramnames <- names(omxGetParameters(model))
if(nrow(cov) != length(paramnames)){ #<--Do its dimensions match the number of free parameters?
msg <- paste("value of argument 'cov' has dimension ",nrow(cov),", but '",model@name,"' has ",length(paramnames)," free parameters",sep="")
stop(msg)
}
covnames <- colnames(cov)
#Do its row and column names match the free-parameter labels (ignoring permutations)?:
if(any(sort(covnames) != sort(paramnames))){
msg <- paste("the dimnames of the matrix provided for argument 'cov' do not match the free-parameter labels of '",model@name,"'",sep="")
stop(msg)
}
covParam <- cov[paramnames,paramnames]
}
}
#Check if single-group model uses RAM expectation, and proceed if so:
if(length(model@submodels)==0){
if (!inherits(model$expectation, "MxExpectationRAM")) {stop(paste("model '",model@name,"' does not use RAM expectation",sep=""))}
return(.mxStandardizeRAMhelper(model=model,SE=SE,ParamsCov=covParam))
}
#Handle multi-group model:
if(length(model@submodels)>0){
out <- NULL
if (inherits(model$expectation, "MxExpectationRAM")) {
out <- list(.mxStandardizeRAMhelper(model=model,SE=SE,ParamsCov=covParam,ignoreSubmodels=TRUE))
names(out)[1] <- model@name
}
if(!inde.subs.flag){
out <- c(out,lapply(
model@submodels[which(
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
)],
.mxStandardizeRAMhelper,SE=SE,ParamsCov=covParam))
if(length(out)==0){stop(paste("model '",model@name,"' does not use RAM expectation",sep=""))}
return(out)
}
else{
out1 <- lapply(
model@submodels[which(
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0})) &
!sapply(model@submodels,function(x){x@independent})
)],
.mxStandardizeRAMhelper,SE=SE,ParamsCov=covParam)
out <- c(out,out1,out2)
if(length(out)==0){stop(paste("model '",model@name,"' contains no submodels that use RAM expectation",sep=""))}
out <- out[names(model@submodels[which(
(sapply(model@submodels,function(x){class(x$expectation)})=="MxExpectationRAM" |
sapply(model@submodels,function(x){length(x@submodels)>0}))
)])]
return(out)
}}}
#Alias using proper camel-casing:
mxStandardizeRAMPaths <- mxStandardizeRAMpaths
mxBootstrapStdizeRAMpaths <- function(model, bq=c(.25,.75), method=c('bcbci','quantile'), returnRaw=FALSE){
bq <- c(min(bq),max(bq))
if(!is(model, "MxModel")) {
stop("'model' argument must be a MxModel object")
}
if(!length(model@expectation) || !inherits(model@expectation, "MxExpectationRAM")) {
msg <- paste(
"MxModel ",omxQuotes(model@name),
" does not use RAM expectation\n(to use mxBootstrapStdizeRAMpaths() on a RAM submodel, run the function directly on that submodel",sep="")
stop(msg)
}
method <- match.arg(method)
realstdpaths <- .mxStandardizeRAMhelper(model=model,SE=FALSE,ParamsCov=NULL,inde.subs.flag=FALSE,ignoreSubmodels=TRUE)
rawParams <- as.matrix(omxGetBootstrapReplications(model))
#The tricky thing is that the output length of mxStandardizeRAMpaths() is not guaranteed to be the same for every replication...
outputlist <- vector("list",nrow(rawParams))
conformableFlag <- TRUE
for(i in 1:nrow(rawParams)){
modelcurr <- omxSetParameters(model,labels=colnames(rawParams),values=rawParams[i,])
stdpaths <- .mxStandardizeRAMhelper(model=modelcurr,SE=FALSE,ParamsCov=NULL,inde.subs.flag=FALSE,ignoreSubmodels=TRUE)
outputlist[[i]] <- stdpaths$Std.Value
names(outputlist[[i]]) <- stdpaths$name
if(conformableFlag && (nrow(stdpaths)!=nrow(realstdpaths) || !all(stdpaths$name==realstdpaths$name)) ){
conformableFlag <- FALSE
}
}
if( !conformableFlag ){
if(returnRaw){
warning("names of nonzero paths varied among bootstrap replications; returning raw list of standardized paths")
return(outputlist)
}
else{stop("names of nonzero paths varied among bootstrap replications, and argument 'returnRaw' is FALSE")}
}
else{
outmtx <- matrix(NA_real_,nrow=nrow(rawParams),ncol=nrow(realstdpaths))
colnames(outmtx) <- realstdpaths$name
for(i in 1:nrow(rawParams)){
outmtx[i,] <- as.vector(outputlist[[i]])
}
if(returnRaw){return(outmtx)}
}
out <- data.frame(realstdpaths$name,realstdpaths$label,realstdpaths$matrix,realstdpaths$row,realstdpaths$col,
realstdpaths$Std.Value,apply(outmtx,2,sd),numeric(length(realstdpaths$name)),numeric(length(realstdpaths$name)))
colnames(out) <- c("name","label","matrix","row","col","Std.Value","Boot.SE",
sprintf("%.1f%%", round(100*min(bq), 1)),sprintf("%.1f%%", round(100*max(bq), 1)))
if(method=="quantile"){
out[,8] <- as.vector(apply(outmtx,2,quantile,probs=min(bq),type=1))
out[,9] <- as.vector(apply(outmtx,2,quantile,probs=max(bq),type=1))
}
else if(method=="bcbci"){
zcrit <- qnorm(bq)
for(i in 1:nrow(realstdpaths)){
z0 <- qnorm(mean(outmtx[,i] <= realstdpaths$Std.Value[i]))
for (qx in 1:2){
phi <- pnorm(2*z0 + zcrit[qx])
out[i,7+qx] <- quantile(outmtx[,i], probs=phi, type=1)
}
}
}
else{warning("unrecognized value provided for argument 'method'")}
rownames(out) <- NULL
return(out)
}
coef.MxModel <- function(object, ...) omxGetParameters(object)
Try the OpenMx package in your browser
Any scripts or data that you put into this service are public.
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.