#
# Copyright 2007-2019 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.
#' @import OpenMx
# #' MxRegularizedModel
# ##'
# ##' This is an internal class and should not be used directly.
# ##'
# ##' @aliases
# ##' $<-,MxRegularizedModel-method
# ##' [[<-,MxRegularizedModel-method
setClass(Class = "MxRegularizedModel",
slots = representation(
regularizations = "list"
),
contains = "MxModel")
publicMxRegularizedModelSlots <- c("name", "submodels", "output", "compute", "options", "intervals", "penalty") # accessible via model$
visibleMxRegularizedModelSlots <- c("name", "options", "compute", "output", "intervals", "penalty")
settableMxRegularizedModelSlots <- c("submodels", "compute", "options", "penalty")
setMethod("$", "MxRegularizedModel", function(x, name) {
result <- OpenMx:::imxExtractMethod(x, name)
if(name %in% publicMxRegularizedModelSlots) {
result <- OpenMx:::imxExtractSlot(x, name)
} else {
result <- eval(substitute(x$submodels[[1]]$nrm,list(nrm=name)))
}
return(result)
})
setReplaceMethod("$", "MxRegularizedModel",
function(x, name, value) {
if(name %in% settableMxRegularizedModelSlots) {
return(imxReplaceMethod(x, name, value))
} else {
eval(substitute(x@submodels[[1]]$nrm <- value, list(nrm=name)))
return(x)
}
}
)
imxInitRegularizedModel <- function(submodel1) {
model <- new("MxRegularizedModel")
model@submodels[[submodel1$name]] <- submodel1
if(is.null(model@algebras[['penaltyFunction']])) {
model@algebras[['penaltyFunction']] <- mxAlgebraFromString(paste0(submodel1$name,".fitfunction"), name="penaltyFunction")
}
if(is.null(model@fitfunction)) {
model@fitfunction <- mxFitFunctionAlgebra("penaltyFunction")
}
model@name <- paste0(submodel1$name, "_Reg")
return(model)
}
setMethod("imxVerifyModel", "MxRegularizedModel",
function(model) {
if (length(model@submodels) == 0) {
msg <- paste("The regularized model", omxQuotes(model@name),
"does not regularize anything.")
stop(msg, call. = FALSE)
}
fitfunction <- model@fitfunction
if (is.null(fitfunction) || !is(fitfunction, "MxFitFunctionAlgebra")) {
msg <- paste("The regularized model", omxQuotes(model@name),
"is somehow malformed. Sorry!")
stop(msg, call. = FALSE)
}
if (length(model@submodels) > 0) {
return(all(sapply(model@submodels, imxVerifyModel)))
}
return(TRUE)
})
setMethod("print", "MxRegularizedModel", function(x,...) {
displayRegModel(x, TRUE)
})
#' Create a (regularizable) mxModel
#'
#' Identical to OpenMx::mxModel, except the result is an MxRegularizedModel.
#'
#' @importFrom stats sd
#' @export
#'
mxModel <- function(model = NA, ..., manifestVars = NA, latentVars = NA,
remove = FALSE, independent = NA, type = NA, name = NA) {
retval <- regFirstArgument(model, name) # Reimplemented `firstArgument()` function from MxModel
first <- retval[[1]]
regModel <- retval[[2]]
name <- retval[[3]]
regModel@submodels[[1]] <- OpenMx:::typeArgument(regModel@submodels[[1]], type)
if(remove==TRUE) {
# Removal punts upstream. Use with care.
return(OpenMx::mxModel(regModel, first, ..., name=name, manifestVars=manifestVars,
latentVars=latentVars, remove=remove, independent=independent, type=type))
}
# Filter products of variables.
nam <- names(list(...))
filter <- nam %in% c('product', 'productVars')
if(any(filter)) {
stop("Error NYI: Products of variables are not yet supported by mxRegSEM.")
}
# Handle the magic first argument:
lst <- c(first, list(...))
lst <- unlist(lst)
# Handle each thing as it comes through:
filter <- sapply(lst, is, "MxModel")
submodels <- lst[filter]
lst <- lst[!filter]
filter <- sapply(lst, is, "MxPenalty")
if(any(filter)) {
regModel <- mxRegularizedAddPenalty(regModel, lst[filter])
}
lst <- lst[!filter]
filter <- sapply(lst, is, "MxCompute")
if(any(filter)) {
regModel <- OpenMx::mxModel(regModel, lst[filter])
}
lst <- lst[!filter]
filter <- sapply(lst, is, "MxData")
regModel@submodels[[1]] <- imxModelBuilder(model=regModel@submodels[[1]],
lst=lst, name=name,
manifestVars=manifestVars,
latentVars=latentVars,
productVars=character(),
submodels=submodels,
remove=remove,
independent=independent)
# Handle renaming
regModel <- regUpdateNames(regModel)
return(regModel)
}
regUpdateNames <- function(regModel) {
if(!is(regModel, "MxRegularizedModel")) {
return(regModel)
}
# Handle renaming
regModel@name <- paste0(regModel@submodels[[1]]@name, "_Reg")
# Namespace maintenance. Note: use mxRename to really rename.
names(regModel@submodels)[[1]] <- regModel@submodels[[1]]@name
return(regModel)
}
# Helpers to help build regularized model objects
regFirstArgument <- function(model, name) {
first <- NULL
defaultType <- imxModelTypes[[getOption("mxDefaultType")]]
# Check if it's a regularized model already:
if(is(model, "MxRegularizedModel")) {
regModel <- model
} else if (is(model, "MxModel")) {
regModel <- imxInitRegularizedModel(model)
} else {
# Creating a new model: need to make both model and regmodel here.
if (OpenMx:::single.na(model)) {
# Nothing doing here.
} else if (typeof(model) == "character") {
name <- model
} else if (isS4(model)) {
first <- model
} else {
first <- list(model)
}
if (length(name) > 0 && is.na(name)) {
name <- imxUntitledName()
}
imxVerifyName(name, -1)
model <- new(defaultType, name)
regModel <- imxInitRegularizedModel(model)
}
return(list(first, regModel, name))
}
mxRegularizedAddPenalty <- function(regModel, penalties) {
pNames <- sapply(penalties, getElement, "name")
regNames <- sapply(regModel@regularizations, getElement, "name")
conflicts <- pNames %in% regNames
if(any(conflicts)) {
for(aMatch in which(conflicts)) {
# Remove original regularizer; replace below.
regModel@regularizations[[aMatch]] <- NULL
if(!is.null(regModel@matrices[[paste0(aMatch, "_Params")]])) {
regModel@matrices[[paste0(aMatch, "_Params")]] <- NULL
}
if(!is.null(regModel@matrices[[paste0(aMatch, "_Penalty")]])) {
regModel@matrices[[paste0(aMatch, "_Penalty")]] <- NULL
}
regModel@algebras[[aMatch]] <- NULL
}
}
paramVals <- omxGetParameters(regModel, fetch="values")
# Add algebras for each new penalty:
for(aPenalty in penalties) {
# Register penalty
regModel@regularizations[[aPenalty$name]] <- aPenalty
# Generate parameter vectors and create the penalty matrix
reg_params <- aPenalty$params
penalty_Mat <- NA
if(!is.null(reg_params) && length(reg_params) > 0) {
pVals <- paramVals[reg_params]
paramsFree <- rep(TRUE, length(pVals))
paramsFree[grep("[\\[]", reg_params)] <- FALSE
# Doesn't handle def vars yet.
if(length(grep("data[\\.]", reg_params))>0) {
stop("Error NYI in regularization penalty ", aPenalty$name,
": cannot regularize definition vars.")
}
# Add penalty to the regModel
paramsMat <- mxMatrix("Full", nrow=length(paramsFree),
ncol=1, values=pVals,
free=paramsFree, labels=reg_params,
name=paste0(aPenalty$name, "_Params"))
regModel <- OpenMx::mxModel(regModel, paramsMat)
# Add hyperparameters
if(length(aPenalty@hyperparameters) > 0) {
penaltyMat <- mxMatrix("Full", nrow=length(aPenalty@hyperparameters),
ncol=1, free=FALSE,
values=unlist(aPenalty@hyperparameters),
name=paste0(aPenalty$name, "_Penalty"),
labels=paste(aPenalty@name,
names(aPenalty@hyperparameters),
sep="_"),
)
regModel <- OpenMx::mxModel(regModel, penaltyMat)
}
}
regModel <- OpenMx::mxModel(regModel,
imxRegularizationTypes[[aPenalty@type]](aPenalty, paramsMat, penaltyMat))
}
# Rebuild global penalty algebra:
penaltyString <- paste0(regModel@submodels[[1]]$name, ".fitfunction")
if(!is.null(regModel@regularizations) && length(regModel@regularizations) != 0) {
penaltyString <- paste(penaltyString, names(regModel@regularizations), sep="+")
}
regModel@algebras[["penaltyFunction"]] <- mxAlgebraFromString(penaltyString, name="penaltyFunction")
return(regModel)
}
#' Display a regularized model:
displayRegModel <- function(regmodel, expand = FALSE) {
omxQuotes <- OpenMx::omxQuotes # Better done with namespace package.
model <- regmodel@submodels[[1]]
cat("MxModel", omxQuotes(model@name), '(Regularized)\n')
cat("type :", OpenMx::imxTypeName(model), '\n')
cat("$matrices :", omxQuotes(names(model@matrices)), '\n')
cat("$algebras :", omxQuotes(names(model@algebras)), '\n')
cat("$constraints :", omxQuotes(names(model@constraints)), '\n')
cat("$intervals :", omxQuotes(names(model@intervals)), '\n')
cat("$regularizers :", omxQuotes(names(regmodel@regularizers)), '\n')
# latentVars and manifestVars should really be considered
# an implementation detail of RAM and LISREL type models with paths.
# We currently do not return anything when an attempt
# is made to access these slots using $ notation.
# The proper thing to do is probably to ignore the @ slots
# and extract the variable information from the RAM F matrix dimnames
# or equivalent LISREL matrix because the slots are only
# used by mxPath. Matrix constructed models do not use the @ slots.
if (length(model@latentVars) == 0 || OpenMx::imxTypeName(model) %in% "default") {
cat("$latentVars : none\n")
} else if (is.character(model@latentVars)) {
cat("$latentVars :", omxQuotes(model@latentVars), '\n')
} else {
cat("$latentVars :\n")
print(format(model@latentVars))
}
if (length(model@manifestVars) == 0 || imxTypeName(model) %in% "default") {
cat("$manifestVars : none\n")
} else if (is.character(model@manifestVars)) {
cat("$manifestVars :", omxQuotes(model@manifestVars), '\n')
} else {
cat("$manifestVars :\n")
print(format(model@manifestVars))
}
data <- model@data
if (is.null(data)) {
cat("$data : NULL\n")
} else {
if (is(data, "MxDataDynamic")) {
cat("$data type:", omxQuotes(data@type), '\n')
cat("$data$expectation :", omxQuotes(data@expectation), "\n")
} else {
cat("$data :", nrow(data@observed),
"x", ncol(data@observed), "\n")
if(length(data@means) == 1 && is.na(data@means)) {
cat("$data means : NA\n")
} else {
cat("$data means : 1 x", length(data@means), "\n")
}
cat("$data type:", omxQuotes(data@type), '\n')
}
}
cat("$submodels :", omxQuotes(names(model@submodels)), '\n')
expectation <- model@expectation
fitfunction <- model@fitfunction
regularization <- regmodel@regularization
compute <- model@compute
if (is.null(expectation)) {
expectationType <- "NULL"
} else {
expectationType <- class(expectation)[[1]]
}
cat("$expectation :", expectationType, '\n')
if (is.null(fitfunction)) {
fitfunctionType <- "NULL"
} else {
fitfunctionType <- class(fitfunction)[[1]]
}
cat("$fitfunction :", fitfunctionType, '\n')
if (is.null(compute)) {
computeType <- "NULL"
} else {
computeType <- class(compute)[[1]]
}
cat("$compute :", computeType, '\n')
cat("$independent :", model@independent, '\n')
cat("$options :", OpenMx::printOptions(model@options), '\n')
cat("$output :", length(model@output) > 0, '\n')
if(is(regmodel, "MxRegularizedModel")) {
if (is.null(regularization) || length(regularization) == 0) {
regularizationType <- "NULL"
} else {
regularizationType <- paste(sapply(regularization, getElement, "type"), sep=" & ")
}
cat("$regularization :", regularizationType)
if(length(regmodel@matrices) > 0) {
cat("\n--------REGULARIZATION MATRICES--------\n")
lapply(regmodel@matrices, print)
}
if(length(regmodel@algebras) > 0) {
cat("\n--------REGULARIZATION ALGEBRAS--------\n")
lapply(regmodel@algebras, print)
}
if(!is.null(regmodel@submodels[[1]])) {
if (is.null(regmodel@submodels[[1]]@expectation)) {
regularizationType <- "NULL"
} else {
expectationType <- class(regmodel@submodels[[1]]@expectation)[[1]]
}
cat("\t-regularized expectation :", expectationType, '\n')
if (is.null(regmodel@submodels[[1]]@fitfunction)) {
fitfunctionType <- "NULL"
} else {
fitfunctionType <- class(regmodel@submodels[[1]]@fitfunction)[[1]]
}
cat("\t-regularized fitfunction :", fitfunctionType, '\n')
}
}
if(expand) {
if(length(model@matrices) > 0) {
cat("\n--------MATRICES--------\n")
lapply(model@matrices, print)
}
if(length(model@algebras) > 0) {
cat("\n--------ALGEBRAS--------\n")
lapply(model@algebras, print)
}
if(length(model@constraints) > 0) {
cat("\n--------CONSTRAINTS--------\n")
lapply(model@constraints, print)
}
if(!is.null(model@data) > 0) {
cat("\n--------DATA--------\n")
print(model@data)
}
if(!is.null(model@expectation) > 0) {
cat("\n--------EXPECTATION FUNCTION--------\n")
print(model@expectation)
}
if(!is.null(model@fitfunction) > 0) {
cat("\n--------FIT FUNCTION--------\n")
print(model@fitfunction)
}
if(!is.null(model@compute) > 0) {
cat("\n--------COMPUTE--------\n")
print(model@compute)
}
if(length(model@output) > 0) {
cat("\n--------OUTPUT--------\n")
print(model@output)
}
if(length(model@submodels) > 0) {
cat("\n--------SUBMODELS--------\n")
lapply(model@submodels, print)
}
if(length(model@options) > 0) {
cat("\n--------OPTIONS--------\n")
print(model@options)
}
}
invisible(model)
}
# Overwrite renaming:
##' Not documented. See OpenMx::mxRename.
##' @export
mxRename <- function(model, newname, oldname = NA) {
model <- OpenMx::mxRename(model, newname, oldname)
# Apply namespace corrections.
model <- regUpdateNames(model)
return(model)
}
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)
}
}
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)
}
}
#' computeOptimizationStatistics
#'
#' This is overridden from OpenMx. This function is exported for people
#' who already know what they are doing.
#'
computeOptimizationStatistics <- function(model, numStats, useSubmodels,
saturatedDoF, independenceDoF,
retval) {
# get estimated parameters
estimates <- model@output$estimate
# should saturated/independence models include means?
if(length(model@runstate$datalist)==1){
type <- model@runstate$datalist[[1]]@type
means <- model@runstate$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(model@runstate$datalist[[1]]@type != 'raw'){
nvar <- dim(model@runstate$datalist[[1]]@observed)[2]
} else if( length(model@runstate$expectations) == 1 ) {
nvar <- length(model@runstate$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)
# grab the expectation
obj <- model@runstate$expectation
# 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']] <- OpenMx:::calculateConstraints(model, useSubmodels)
retval[['estimatedParameters']] <- nrow(retval$parameters)
if(any(sapply(obj,function(x){"MxExpectationGREML" %in% class(x)}))){
retval[['estimatedParameters']] <- retval[['estimatedParameters']] +
sum(sapply(obj,OpenMx:::imxExtractSlot,name="numFixEff"))
}
if(!is.null(retval[['regularizedParameters']])) {
retval$estimatedParameters <- retval$estimatedParameters-retval$regularizedParameters
}
# browser()
# TODO: Intervene here!
if (is.null(numStats)) {
retval[['observedStatistics']] <- OpenMx:::observedStatistics(model, useSubmodels, 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){"MxExpectationGREML" %in% class(x)}))){
retval[['saturatedDoF']] <- NA
}
# calculate or populate independence degrees of freedom
if(is.null(independenceDoF)) {
if(!any(sapply(obj,function(x){"MxExpectationGREML" %in% class(x)}))){
# 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 <- OpenMx:::fitStatistics(model, useSubmodels, retval)
return(retval)
}
#' Summarize an MxRegularizedModel
#'
#' Create a summary of an MxRegularizedModel.
#'
#' @param regFit the regularized model to summarize
#' @param ... values passed to mxSummarize
#' @param verbose gives more summary feedback
#' @param epsilon how close to zero is close enough? Regularization limits for DF correction
#' @return Summary object summarizing the MxRegularizedModel
#' @import OpenMx
#' @export
setMethod("summary", "MxRegularizedModel", function(object, ..., verbose=FALSE, computeRefs=FALSE, epsilon = 1e-6) {
# summary.MxRegularizedModel <- function(object, ..., verbose=FALSE, epsilon = 1e-6) {
regFit <- object
# default penalty_function ("guess") looks for a "lasso" object in the container model.
dotArguments <- list(...)
if(length(regFit$submodels) == 1) {
model <- regFit$submodels[[1]]
} else {
model <- regFit
}
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")
if(computeRefs && (is.null(saturatedDoF) || is.null(independenceDoF))) {
refModels <- mxRefModels(regFit$submodels[[1]], run=FALSE)
satModel <- OpenMx::mxRun(refModels[['Saturated']], silent=TRUE)
indModel <- OpenMx::mxRun(refModels[['Independence']], silent=TRUE)
saturatedLikelihood <- refToLikelihood(satModel)
saturatedDoF <- refToDof(satModel)
independenceLikelihood <- refToLikelihood(indModel)
independenceDoF <- refToDof(indModel)
}
}
numObs <- dotArguments$numObs
numStats <- dotArguments$numStats
# Handle regularization elements
regList <- c()
for(penalty in regFit@regularizations) {
# TODO: Manage uniqueness in a better way
regged <- omxGetParameters(regFit, free=TRUE, labels=penalty@params)
zeroedOut <- names(regged)[abs(regged) < epsilon]
regList <- union(regList, zeroedOut)
}
# Win back free parameters for any that's less than zero.
regOffset <- length(regList)
useSubmodels <- dotArguments$indep
if (is.null(useSubmodels)) { useSubmodels <- TRUE }
retval <- list(wasRun=model@.wasRun, stale=model@.modifiedSinceRun)
retval$parameters <- OpenMx:::parameterList(regFit, useSubmodels)
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 <- regFit@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 <-
OpenMx:::summarizeBootstrap(retval$parameters[, 'Estimate'], bootData, bq, summaryType)
}
}
} else if (any(grep('^boot\\.', names(dotArguments)))) {
warning("No bootstrap data found. See ?mxBootstrap")
}
retval$GREMLfixeff <- OpenMx:::GREMLFixEffList(model)
retval$infoDefinite <- regFit@output$infoDefinite
retval$conditionNumber <- regFit@output$conditionNumber
if (length(regFit@output$gradient)) {
agrad <- abs(regFit@output$gradient)
retval$maxAbsGradient <- agrad[ order(-agrad)[1] ]
}
retval <- OpenMx:::boundsMet(model, retval)
retval <- OpenMx:::setLikelihoods(regFit, saturatedLikelihood, independenceLikelihood, retval)
retval[["Minus2LogLikelihood"]] <- mxEvalByName(paste0(model$name, ".fitfunction"), regFit, compute=TRUE)[1,1]
retval <- OpenMx:::setNumberObservations(numObs, regFit@runstate$datalist, regFit@runstate$fitfunctions, retval)
retval[["regularizedParameters"]] <- regOffset
retval[["zeroedParameters"]] <- regList
retval <- computeOptimizationStatistics(regFit, numStats, useSubmodels, saturatedDoF, independenceDoF, retval)
retval$dataSummary <- OpenMx:::generateDataSummary(model, useSubmodels)
retval$CI <- as.data.frame(regFit@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 <- regFit@output$confidenceIntervalCodes
statusCode <- regFit@output$status$code
if (!is.null(statusCode)) {
message <- OpenMx:::optimizerMessages[[as.character(statusCode)]]
retval[['npsolMessage']] <- message
retval[['statusCode']] <- OpenMx:::as.statusCode(statusCode)
retval[['maxRelativeOrdinalError']] <- regFit@output[['maxRelativeOrdinalError']]
}
if( .hasSlot(regFit,"compute") && length(regFit$compute$steps$CI) ){
retval$CIdetail <- regFit$compute$steps$CI$output$detail
}
retval$timestamp <- regFit@output$timestamp
retval$frontendTime <- regFit@output$frontendTime
retval$backendTime <- regFit@output$backendTime
retval$independentTime <- regFit@output$independentTime
retval$wallTime <- regFit@output$wallTime
retval$cpuTime <- regFit@output$cpuTime
retval$mxVersion <- regFit@output$mxVersion
retval$modelName <- model@name
plan <- regFit@runstate$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)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.