Nothing
#
# 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.
##' MxInterval
##'
##' @description
##' This is an internal class and should not be used directly.
##'
##' @aliases
##' $,MxInterval-method
##' $<-,MxInterval-method
##' print,MxInterval-method
##' show,MxInterval-method
##' @seealso
##' \link{mxCI}
##' @rdname MxInterval-class
setClass(Class = "MxInterval",
representation = representation(
reference = "character",
lowerdelta = "numeric",
upperdelta = "numeric",
boundAdj = "logical"
))
setMethod("initialize", "MxInterval",
function(.Object, reference, lowerdelta, upperdelta, boundAdj) {
.Object@reference <- reference
.Object@lowerdelta <- lowerdelta
.Object@upperdelta <- upperdelta
.Object@boundAdj <- boundAdj
return(.Object)
}
)
createNewInterval <- function(reference, lowerdelta, upperdelta, boundAdj) {
return(new("MxInterval", reference, lowerdelta, upperdelta, boundAdj))
}
setMethod("$", "MxInterval", imxExtractSlot)
setReplaceMethod("$", "MxInterval",
function(x, name, value) {
return(imxReplaceSlot(x, name, value, check=TRUE))
}
)
setMethod("names", "MxInterval", slotNames)
expandSingleInterval <- function(interval) {
references <- interval@reference
if (length(references) == 1) {
return(interval)
} else {
return(lapply(references, createNewInterval,
interval@lowerdelta, interval@upperdelta, interval@boundAdj))
}
}
expandIntervals <- function(intervals) {
if (length(intervals) == 0) {
return(intervals)
}
retval <- lapply(intervals, expandSingleInterval)
retval <- unlist(retval)
return(retval)
}
mxCI <- function(reference, interval = 0.95, type = c('both', 'lower', 'upper'), ..., boundAdj=TRUE) {
prohibitDotdotdot(list(...))
if (!is.numeric(interval) || interval < 0 || interval > 1) {
stop("'interval' must be a numeric value between 0 and 1")
}
if (identical(type, c('both', 'lower', 'upper'))) {
type <- 'both'
}
if (!is.character(type) || length(type) != 1 || !(type %in% c('both', 'lower', 'upper'))) {
stop("'type' must be either 'both' or 'lower' or 'upper'")
}
if (type == 'both') {
lowerValue <- qchisq(interval, 1)
upperValue <- qchisq(interval, 1)
} else if (type == 'lower') {
lowerValue <- qchisq(interval, 1)
upperValue <- as.numeric(NA)
} else if (type == 'upper') {
lowerValue <- as.numeric(NA)
upperValue <- qchisq(interval, 1)
}
return(confidenceIntervalHelper(reference, lowerValue, upperValue, boundAdj))
}
confidenceIntervalHelper <- function(reference, lowerdelta, upperdelta, boundAdj) {
if (single.na(lowerdelta)) { lowerdelta <- as.numeric(NA) }
if (single.na(upperdelta)) { upperdelta <- as.numeric(NA) }
if (!is.character(reference) || length(reference) < 1 || any(is.na(reference))) {
stop("'reference' argument must be a character vector")
}
if (!is.na(lowerdelta) && (!is.numeric(lowerdelta) || length(lowerdelta) != 1 || lowerdelta < 0)) {
stop("'lowerdelta' argument must be a non-negative numeric value")
}
if (!is.na(upperdelta) && (!is.numeric(upperdelta) || length(upperdelta) != 1 || upperdelta < 0)) {
stop("'upperdelta' argument must be a non-negative numeric value")
}
retval <- createNewInterval(reference, lowerdelta, upperdelta, as.logical(boundAdj))
return(retval)
}
modelAddIntervals <- function(model, intervals) {
if (length(intervals) == 0) {
return(model)
}
iNames <- names(intervals)
for(i in 1:length(intervals)) {
model@intervals[[iNames[[i]]]] <- intervals[[i]]
}
return(model)
}
modelRemoveIntervals <- function(model, intervals) {
if (length(intervals) == 0) {
return(model)
}
iNames <- names(intervals)
for(i in 1:length(intervals)) {
model@intervals[[iNames[[i]]]] <- NULL
}
return(model)
}
generateIntervalList <- function(flatModel, modelname, parameters, labelsData) {
retval <- lapply(flatModel@intervals, generateIntervalListHelper,
flatModel, modelname, parameters, labelsData)
names(retval) <- NULL
retval <- unlist(retval, recursive = FALSE)
# Might have requested a CI of the same parameter in different submodels
retval <- retval[!duplicated(names(retval))]
if (length(retval) > 0) {
approveWLSIntervals(flatModel, modelname)
}
return(retval)
}
makeIntervalReference <- function(entityNumber, row, col, iobj) {
return(c(entityNumber, row - 1, col - 1, iobj@lowerdelta, iobj@upperdelta, iobj@boundAdj))
}
generateIntervalListHelper <- function(interval, flatModel, modelname,
parameters, labelsData) {
reference <- interval@reference
# length(reference)==1 because of expandSingleInterval
pindex <- match(reference, names(parameters))
if (!is.na(pindex)) {
location <- parameters[[pindex]][[5]]
location[[1]] <- - location[[1]] - 1
retval <- list()
retval[[reference]] <-
c(location, interval@lowerdelta, interval@upperdelta, interval@boundAdj)
return(retval)
}
entity <- flatModel[[reference]]
retval <- list()
if (!is.null(entity)) {
entityNumber <- imxLocateIndex(flatModel, reference,
paste("confidence interval", reference))
if (is(entity, "MxAlgebra")) {
newName <- reference
retval[[newName]] <- makeIntervalReference(entityNumber, 0L, 0L, interval)
} else {
if(!any(entity@free)){return(retval)}
rows <- nrow(entity)
cols <- ncol(entity)
free <- entity@free
for(i in 1:rows) {
for(j in 1:cols) {
if (imxSymmetricMatrix(entity) && i > j) next
if (!free[i, j]) next
if(!is.na(entity@labels[i,j]) ){
newName <- entity@labels[i,j]
} else {
newName <- paste(reference, '[', i, ',', j, ']', sep = '')
}
retval[[newName]] <- makeIntervalReference(entityNumber, i, j, interval)
}
}
}
return(retval)
} else if (hasSquareBrackets(reference)) {
components <- splitSubstitution(reference)
entityName <- components[[1]]
rows <- eval(parse(text = components[[2]]), envir = globalenv())
cols <- eval(parse(text = components[[3]]), envir = globalenv())
entityValue <- eval(substitute(mxEval(x, flatModel, compute=TRUE, .extraBack=7L),
list(x = as.symbol(entityName))))
if (is.null(rows)) {
rows <- 1:nrow(entityValue)
}
if (is.null(cols)) {
cols <- 1:ncol(entityValue)
}
entity <- flatModel[[entityName]]
entityNumber <- imxLocateIndex(flatModel, entityName,
paste("confidence interval", reference))
retval <- list()
for(i in rows) {
for(j in cols) {
newName <- paste(entityName, '[', i, ',', j, ']', sep = '')
retval[[newName]] <- makeIntervalReference(entityNumber, i, j, interval)
}
}
return(retval)
} else {
for (entityName in names(flatModel@matrices)) {
entity <- flatModel[[entityName]]
free <- entity@free
rows <- nrow(free)
cols <- ncol(free)
for(i in 1:rows) {
for(j in 1:cols) {
if (free[i, j]) next
label <- entity@labels[i,j]
if (is.na(label) || label != reference) next
entityNumber <- imxLocateIndex(flatModel, entityName,
paste("confidence interval", reference))
retval[[label]] <- makeIntervalReference(entityNumber, i, j, interval)
return(retval)
}
}
}
stop(paste("Unknown reference to", omxQuotes(reference),
"detected in a confidence interval",
"specification in model", omxQuotes(modelname),
"\nYou should check spelling (case-sensitive), and also addressing the right model: to refer to an algebra",
"\nSee help(mxCI) to see how to refer to an algebra in a submodel.\n",
"FYI, I got as far as: ", deparse(width.cutoff = 400L, sys.call(-3))), call. = FALSE)
}
}
displayInterval <- function(object) {
cat("MxInterval", '\n')
cat("$reference: ", object@reference, '\n')
cat("$lowerdelta: ", object@lowerdelta, '\n')
cat("$upperdelta: ", object@upperdelta, '\n')
cat("$boundAdj: ", object@boundAdj, '\n')
}
setMethod("print", "MxInterval", function(x,...) { displayInterval(x) })
setMethod("show", "MxInterval", function(object) { displayInterval(object) })
expandConfidenceIntervals <- function(model, intervals) {
expansion <- lapply(intervals, expandConfidenceIntervalsHelper, model)
return(unlist(expansion))
}
createSimilarInterval <- function(name, template) {
new("MxInterval", name, template@lowerdelta, template@upperdelta, template@boundAdj)
}
expandConfidenceIntervalsHelper <- function(interval, model) {
reference <- interval@reference
entity <- model[[reference]]
retval <- list()
if(!is.null(entity)) {
entityValue <- eval(substitute(mxEval(x, model, compute=TRUE),
list(x = as.symbol(reference))))
rows <- nrow(entityValue)
cols <- ncol(entityValue)
if (is(entity, "MxMatrix")) {
free <- entity@free
if (imxSymmetricMatrix(entity)) {
free[lower.tri(free, diag=FALSE)] <- FALSE
}
} else {
free <- matrix(TRUE, rows, cols)
}
for(i in 1:rows) {
for(j in 1:cols) {
if (free[i, j]) {
newName <- paste(reference, '[', i, ',', j, ']', sep = '')
retval <- c(retval, createSimilarInterval(newName, interval))
}
}
}
} else if (hasSquareBrackets(reference)) {
components <- splitSubstitution(reference)
entityName <- components[[1]]
rows <- eval(parse(text = components[[2]]), envir = globalenv())
cols <- eval(parse(text = components[[3]]), envir = globalenv())
entityValue <- eval(substitute(mxEval(x, model, compute=TRUE),
list(x = as.symbol(entityName))))
if (is.null(rows)) {
rows <- 1:nrow(entityValue)
}
if (is.null(cols)) {
cols <- 1:ncol(entityValue)
}
entity <- model[[entityName]]
if (is.null(rows)) {
rows <- 1:nrow(entityValue)
}
if (is.null(cols)) {
cols <- 1:ncol(entityValue)
}
entity <- model[[entityName]]
if (is(entity, "MxMatrix")) {
free <- entity@free
if (imxSymmetricMatrix(entity)) {
free[lower.tri(free, diag=FALSE)] <- FALSE
}
} else {
free <- matrix(TRUE, max(rows), max(cols))
}
for(i in rows) {
for(j in cols) {
if (free[i, j]) {
newName <- paste(entityName, '[', i, ',', j, ']', sep = '')
retval <- c(retval, createSimilarInterval(newName, interval))
}
}
}
} else {
retval <- interval
}
return(retval)
}
removeAllIntervals <- function(model) {
model@intervals <- list()
model@submodels <- lapply(model@submodels, removeAllIntervals)
return(model)
}
omxParallelCI <- function(model, run = TRUE, verbose=0L, independentSubmodels=TRUE,
optimizer=mxOption(NULL,"Default optimizer")) {
if(missing(model) || !is(model, "MxModel")) {
stop("first argument must be a MxModel object")
}
assertModelRunAndFresh(model)
if(length(model@output) == 0) {
stop("'model' argument to omxParallelCI must be a fitted model")
}
if( !(optimizer %in%c("NPSOL","CSOLNP","SLSQP")) ){
stop("argument 'optimizer' must be one of 'NPSOL', 'CSOLNP', or 'SLSQP'")
}
if(!independentSubmodels){
optionList <- generateOptionsList(model, TRUE)
ctype <- ifelse(optimizer=="SLSQP","ineq","none")
ciOpt <- mxComputeGradientDescent(
verbose=verbose, engine=optimizer,
nudgeZeroStarts=FALSE,
gradientAlgo=optionList[['Gradient algorithm']],
gradientIterations=imxAutoOptionValue('Gradient iterations',optionList),
gradientStepSize=imxAutoOptionValue('Gradient step size',optionList))
if (ctype == 'ineq') {
ciOpt <- mxComputeTryHard(plan=ciOpt, scale=0.05, verbose=verbose)
}
pciplan <- mxComputeSequence(
steps=list(
CI=mxComputeConfidenceInterval(
constraintType=ctype,
verbose=verbose,
plan=ciOpt
)))
pciplan
model <- mxModel(model, pciplan)
if(!run){
return(model)
}
model <- mxRun(model, suppressWarnings=TRUE)
model@compute@.persist <- FALSE
return(model)
}
else{
namespace <- imxGenerateNamespace(model)
flatModel <- imxFlattenModel(model, namespace, FALSE)
intervals <- flatModel@intervals
if (length(intervals) == 0) return(model)
intervals <- expandConfidenceIntervals(model, intervals)
template <- model
template <- removeAllIntervals(template)
optionList <- generateOptionsList(model, TRUE)
ctype <- ifelse(optimizer=="SLSQP","ineq","none")
ciOpt <- mxComputeGradientDescent(
verbose=verbose, engine=optimizer,
nudgeZeroStarts=FALSE,
gradientAlgo=optionList[['Gradient algorithm']],
gradientIterations=imxAutoOptionValue('Gradient iterations',optionList),
gradientStepSize=imxAutoOptionValue('Gradient step size',optionList))
if (ctype == 'ineq') {
ciOpt <- mxComputeTryHard(plan=ciOpt, scale=0.05, verbose=verbose)
}
pciplan <- mxComputeSequence(
steps=list(
CI=mxComputeConfidenceInterval(
constraintType=ctype,
verbose=verbose,
plan=ciOpt
)))
modelname <- model@name
container <- mxModel(paste(modelname, "container", sep = "_"))
submodels <- list()
for(i in 1:length(intervals)) {
interval <- intervals[[i]]
newmodel <- mxModel(template, interval, pciplan, independent = TRUE)
newmodel <- mxRename(newmodel, paste("interval", i, sep = ""))
newmodel <- mxOption(newmodel, "Number of Threads", 1)
submodels <- c(submodels, newmodel)
}
container <- mxModel(container, submodels)
if (!run) {
return(container)
}
container <- mxRun(container, intervals = TRUE, suppressWarnings = TRUE)
tableCI <- matrix(as.numeric(NA), 0, 3)
tableCodes <- matrix(0, 0, 2)
dimnames(tableCI) <- list(NULL, c('lbound', 'estimate', 'ubound'))
dimnames(tableCodes) <- list(NULL, c('lbound', 'ubound'))
submodels <- container@submodels
for(i in 1:length(submodels)) {
submodel <- submodels[[i]]
submodel <- mxRename(submodel, modelname)
tableCI <- rbind(tableCI, submodel@output$confidenceIntervals)
tableCodes <- rbind(tableCodes, submodel@output$confidenceIntervalCodes)
}
model@output$confidenceIntervals <- tableCI
model@output$confidenceIntervalCodes <- tableCodes
model@output$frontendTime <- container@output$frontendTime
model@output$backendTime <- container@output$backendTime
model@output$independentTime <- container@output$independentTime
model@output$wallTime <- container@output$wallTime
model@output$timestamp <- container@output$timestamp
model@output$cpuTime <- container@output$cpuTime
return(model)
}
}
omxRunCI <- function(model, verbose=0L, optimizer="SLSQP"){
return(omxParallelCI(model=model, run=TRUE, verbose=verbose, independentSubmodels=FALSE, optimizer=optimizer))
}
confint.MxModel <- function(object, parm, level = 0.95, ...) {
allParam <- names(coef(object))
if (missing(parm)) parm <- allParam
assertModelRunAndFresh(object)
se <- c(object$output[['standardErrors']])
if (is.null(se)) stop("No standard errors are available")
names(se) <- allParam
notFound <- is.na(match(parm, allParam))
if (any(notFound)) {
stop(paste("Parameter", omxQuotes(parm[notFound]),
"not recognized"))
}
se <- se[parm]
if (any(is.na(se))) stop(paste("Some standard errors are NA:",
omxQuotes(names(se)[is.na(se)])))
est <- coef(object)[parm]
alpha <- (1-level)/2
thresh <- qnorm(1 - alpha)
df <- data.frame(lower=est - se * thresh,
upper=est + se * thresh)
rownames(df) <- names(est)
colnames(df) <- sapply(c(alpha, 1-alpha),
function(x) sprintf("%.1f%%", round(100*min(x), 1)))
message("Wald type confidence intervals (see ?mxCI for likelihood-based CIs)")
df
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.