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.
setClass(Class = "BaseExpectationNormal",
contains = "MxBaseExpectation",
representation = representation(
expectedCovariance = "MxOptionalCharOrNumber",
expectedMean = "MxOptionalCharOrNumber",
expectedFullCovariance = "MxOptionalCharOrNumber",
expectedFullMean = "MxOptionalCharOrNumber",
thresholds = "MxCharOrNumber",
threshnames = "character",
discrete = "MxCharOrNumber",
discreteSpec = "MxOptionalMatrix",
.discreteCheckCount = "logical",
selectionPlan = "MxOptionalDataFrame",
selectionVector = "MxCharOrNumber"))
setMethod("qualifyNames", signature("BaseExpectationNormal"),
function(.Object, modelname, namespace) {
for (sl in c(paste0('expected', c('Covariance','Mean')),
'thresholds', 'discrete', 'selectionVector')) {
if (!.hasSlot(.Object, sl) || is.null(slot(.Object, sl))) next
slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
}
.Object
})
setMethod("genericExpDependencies", signature("BaseExpectationNormal"),
function(.Object, dependencies) {
sources <- c(.Object$discrete, .Object@thresholds, .Object$selectionVector)
sources <- sources[!is.na(sources)]
dependencies <- imxAddDependency(sources, .Object@name, dependencies)
return(dependencies)
})
setMethod("genericExpRename", signature("BaseExpectationNormal"),
function(.Object, oldname, newname) {
for (sl in c(paste0('expected', c('Covariance','Mean')),
'thresholds', 'discrete', 'selectionVector')) {
if (!.hasSlot(.Object, sl) || is.null(slot(.Object, sl))) next
slot(.Object,sl) <- renameReference(slot(.Object,sl), oldname, newname)
}
.Object
})
setMethod("genericExpFunConvert", signature("BaseExpectationNormal"),
function(.Object, flatModel, model, labelsData, dependencies) {
for (sl in c(paste0('expected', c('Covariance','Mean')))) {
matName <- slot(.Object,sl)
if (length(matName) == 0) next
mat <- flatModel[[matName]]
if (any(mat@free)) {
msg <- paste('Free parameters are not allowed in matrix',
omxQuotes(matName))
stop(msg, call.=FALSE)
}
if (any(!is.na(mat@labels))) {
msg <- paste('Labels are not allowed in matrix', omxQuotes(matName))
stop(msg, call.=FALSE)
}
}
tname <- .Object$thresholds
dname <- .Object$discrete
if (!is.na(dname)) {
dmat <- flatModel[[dname]]
if (length(colnames(dmat)) == 0) {
stop(paste(dname, "must have column names"))
}
ds <- .Object$discreteSpec
if (is.null(ds)) stop("discrete parameters but no discreteSpec")
if (length(colnames(ds)) != length(colnames(dmat))) {
stop(paste("discreteSpec must have the same number of columns",
"as", dname))
}
if (any(colnames(ds) != colnames(dmat))) {
stop(paste("discreteSpec and", dname, "must have the same column names"))
}
}
if (!is.na(tname) && !is.na(dname)) {
dup <- intersect(colnames(flatModel[[tname]]), colnames(flatModel[[dname]]))
if (length(dup)) {
msg <- paste('Manifest', omxQuotes(dup), 'cannot be modelled as',
'a threshold and as a discrete variable at the',
'same time')
stop(msg, call.=FALSE)
}
}
if (.hasSlot(.Object, 'selectionVector')) {
svname <- .Object$selectionVector
if (!is.na(svname)) {
sv <- flatModel[[svname]]
if (ncol(sv) != 1) stop("selectionVector must be a column vector")
svSpec <- .Object@selectionPlan
if (nrow(sv) != nrow(svSpec)) {
msg <- paste('selectionPlan must have the same number of rows',
'as the selectionVector')
stop(msg, call.=FALSE)
}
}
}
.Object
})
setMethod("genericNameToNumber", signature("BaseExpectationNormal"),
function(.Object, flatModel, model) {
name <- .Object@name
for (sl in c(paste0('expected', c('Covariance','Mean')),
'data', 'thresholds', 'discrete', 'selectionVector')) {
if (!.hasSlot(.Object, sl)) next
slot(.Object,sl) <- imxLocateIndex(flatModel, slot(.Object,sl), name)
}
.Object
})
setMethod("genericGetExpected", signature("BaseExpectationNormal"),
function(.Object, model, what, defvar.row=1, subname=model@name) {
ret <- list()
if ('thresholds' %in% what) {
thrname <- .Object@thresholds
thr <- matrix( , 0 , 0)
if(!single.na(thrname)){
thrname <- .modifyDottedName(subname, thrname, sep=".")
thr <- mxEvalByName(thrname, model, compute=TRUE, defvar.row=defvar.row)
tnames <- .Object@threshnames
if(!single.na(tnames)){
colnames(thr) <- tnames
}
}
disname <- .Object@discrete
if (!single.na(disname)) {
npart <- mxGetExpected(model, c('means','covariance'),
defvar.row=defvar.row, subname=subname)
means <- npart[['means']]
if (!is.null(means) && ncol(means)==1 && colnames(means)=='one')
means <- t(means)
cov <- npart[['covariance']]
ds <- .Object@discreteSpec
disname <- .modifyDottedName(subname, disname, sep=".")
dis <- mxEvalByName(disname, model, compute=TRUE, defvar.row=defvar.row)
if (length(colnames(ds)) != length(colnames(dis))) {
stop(paste("discreteSpec must have the same number of columns",
"as", disname))
}
if (any(colnames(ds) != colnames(dis))) {
stop(paste("discreteSpec and", disname, "must have the same column names"))
}
# discreteSpec:
# number of outcomes
# 1 poisson ; 2 negative binomial prob ; 3 neg binomial mu
# discrete 3 rows:
# zero inflation
# poisson: lambda; nbinom: size, prob/mu
thList <- vector("list", ncol(ds))
names(thList) <- colnames(ds)
for (cx in colnames(ds)) {
if (is.na(ds[1,cx])) {
stop(paste("discreteSpec[1,",cx,"], the maximum observed count,",
"cannot be NA when calcuating the expected distribution"))
}
outcome <- 1L:as.integer(ds[1,cx]) - 1L
zif <- min(max(dis[1,cx],0),1)
if (ds[2,cx] == 1) {
pr <- ppois(outcome, dis[2,cx])
} else if (ds[2,cx] == 2) {
pr <- pnbinom(outcome, dis[2,cx], prob=dis[3,cx])
} else if (ds[2,cx] == 3) {
pr <- pnbinom(outcome, dis[2,cx], mu=dis[3,cx])
} else { stop(paste("Unknown discrete distribution code", ds[2,cx],
"in column", cx)) }
m1 <- 0
if (!is.null(means)) m1 <- means[1,cx]
thList[[cx]] <- p2z(zif + (1-zif) * pr) * sqrt(cov[cx,cx]) + m1
}
thBlock <- mxSimplify2Array(thList)
comb <- matrix(NA,
nrow=max(nrow(thr), nrow(thBlock)),
ncol=(ncol(thr) + ncol(thBlock)))
if (length(thr)) {
comb[1:nrow(thr),1:ncol(thr)] <- thr
comb[1:nrow(thBlock),ncol(thr) + 1:ncol(thBlock)] <- thBlock
colnames(comb) <- c(colnames(thr), colnames(thBlock))
thr <- comb
} else {
thr <- thBlock
}
}
ret[['thresholds']] <- thr
}
ret
})
constrainCorData <- function(ex, size, model, flatModel) {
path <- unlist(strsplit(ex@name, imxSeparatorChar, fixed = TRUE))
if (path[1] == model@name) {
submodel <- model
} else {
submodel <- model[[ path[1] ]]
}
newObj <- model@.newobjects
# Don't deal with inherited MxData. That's an indication of
# an advanced user who knows what they are doing.
if (!is.null(submodel@data) && submodel@data@type == 'cor' && is.null(ex@expectedCovariance)) {
for (sl in c('expectedCovariance', 'unit_1_by_nObs', 'constraintForCorData')) {
if (!is.null(submodel[[sl]])) {
msg <- paste("In model", omxQuotes(submodel@name),
"data has type='cor' with", omxQuotes(class(ex)[1]),
"but an object named", omxQuotes(sl),
"already exists preventing the addition of a constraint",
"on the diagonal of the expected covariance")
stop(msg, call.=FALSE)
}
}
ex@expectedCovariance <- 'expectedCovariance'
corConstraint <- list(
ex,
mxMatrix(name="expectedCovariance", nrow=size, ncol=size),
mxMatrix(name= "unit_1_by_nObs", type="Unit", ncol=1, nrow=size),
mxConstraintFromString(name="constraintForCorData", "diag2vec(expectedCovariance) == unit_1_by_nObs"))
submodel <- mxModel(submodel, corConstraint)
newObj <- TRUE
}
if (path[1] == model@name) {
model <- submodel
} else {
model <- mxModel(model, submodel)
}
model@.newobjects <- newObj
model
}
setClass(Class = "MxExpectationNormal",
representation = representation(
covariance = "MxCharOrNumber",
means = "MxCharOrNumber",
dims = "character",
ExpCov = "matrix",
ExpMean = "matrix",
numStats = "numeric"),
contains = "BaseExpectationNormal")
setMethod("initialize", "MxExpectationNormal",
function(.Object, covariance, means, dims, thresholds, threshnames, discrete,
discreteSpec, selectionVector, selectionPlan,
data = as.integer(NA), name = 'expectation') {
.Object@name <- name
.Object@covariance <- covariance
.Object@means <- means
.Object@data <- data
.Object@thresholds <- thresholds
.Object@dims <- dims
.Object@threshnames <- threshnames
.Object@discrete <- discrete
.Object@discreteSpec <- discreteSpec
.Object@.discreteCheckCount <- TRUE
.Object@selectionVector <- selectionVector
.Object@selectionPlan <- selectionPlan
.Object@ExpCov <- matrix()
.Object@ExpMean <- matrix()
return(.Object)
}
)
setMethod("qualifyNames", signature("MxExpectationNormal"),
function(.Object, modelname, namespace) {
.Object <- callNextMethod()
.Object@name <- imxIdentifier(modelname, .Object@name)
for (s in c("covariance", "means", "data")) {
if (is.null(slot(.Object, s))) next
slot(.Object, s) <-
imxConvertIdentifier(slot(.Object, s), modelname, namespace)
}
return(.Object)
})
setMethod("genericExpDependencies", signature("MxExpectationNormal"),
function(.Object, dependencies) {
dependencies <- callNextMethod()
sources <- c(.Object@covariance, .Object@means)
sources <- sources[!is.na(sources)]
dependencies <- imxAddDependency(sources, .Object@name, dependencies)
return(dependencies)
})
setMethod("genericExpConvertEntities", "MxExpectationNormal",
function(.Object, flatModel, namespace, labelsData) {
flatModel <- updateExpectationDimnames(.Object, flatModel, labelsData)
flatModel <- updateThresholdDimnames(.Object, flatModel, labelsData)
return(flatModel)
}
)
setMethod("genericExpRename", signature("MxExpectationNormal"),
function(.Object, oldname, newname) {
.Object <- callNextMethod()
for (s in c("covariance", "means", "data")) {
if (is.null(slot(.Object, s))) next
slot(.Object, s) <- renameReference(slot(.Object, s), oldname, newname)
}
return(.Object)
})
.modifyDottedName <- function(sub, obj, sep='.'){
if(single.na(obj)) {return(obj)}
if(length(obj) > 1) {stop('Please pass me one object name at a time.')}
if(length(strsplit(obj, split=sep, fixed=TRUE)[[1]]) > 1 ){
return(obj)
} else {
return(paste(sub, obj, sep=sep))
}
}
setMethod("genericGetExpected", signature("MxExpectationNormal"),
function(.Object, model, what, defvar.row=1, subname=model@name) {
ret <- callNextMethod()
if (any(c('covariance','covariances') %in% what)) {
covname <- .modifyDottedName(subname, .Object@covariance)
cov <- mxEvalByName(covname, model, compute=TRUE, defvar.row=defvar.row, .extraBack=7L)
dnames <- .Object$dims
if(!single.na(dnames)){
colnames(cov) <- dnames
rownames(cov) <- dnames
}
ret[['covariance']] <- cov
}
if (any(c('means', 'mean') %in% what)) {
meanname <- .Object@means
if(!single.na(meanname)){
meanname <- .modifyDottedName(subname, meanname, sep=".")
mean <- mxEvalByName(meanname, model, compute=TRUE, defvar.row=defvar.row)
dnames <- .Object$dims
if(!single.na(dnames)){
colnames(mean) <- dnames
}
} else {mean <- matrix( , 0, 0)}
ret[['means']] <- mean
}
ret
})
getThresholdMask <- function(model, cols, subname) {
if (subname != model$name) {
d1 <- model[[subname]]@data
} else {
d1 <- model@data
}
if (is.null(d1)) stop(paste("Cannot find observed thresholds, model",
omxQuotes(subname), "has no data"))
if (d1@type == 'raw') {
lev <- sapply(d1$observed[,cols], function(x) length(levels(x)))
mask <- matrix(FALSE, max(lev)-1, length(cols))
colnames(mask) <- cols
for (lx in 1:length(cols)) mask[1:(lev[lx]-1), lx] <- TRUE
mask
} else {
th <- d1@thresholds
if (single.na(th)) stop(paste("Observed data in model", omxQuotes(model$name), "has no thresholds"))
!is.na(th)
}
}
setMethod("genericGetExpectedVector", signature("BaseExpectationNormal"),
function(.Object, model, defvar.row=1, subname=model@name) {
ret <- genericGetExpected(.Object, model, c('covariance', 'means', 'thresholds'), defvar.row, subname)
cov <- ret[['covariance']]
mns <- ret[['means']]
if (is.null(mns)) stop("mns is null")
nv <- nrow(cov)
covNames <- paste0('cov', vech(outer(1:nv, 1:nv, FUN=paste, sep='_')))
mnsNames <- paste0('mean', 1:nv)
thr <- ret[['thresholds']]
if (prod(dim(thr)) == 0) {
v <- c(vech(cov), mns[!is.na(mns)])
names(v) <- c(covNames, mnsNames[!is.na(mns)])
} else {
thrNames <- outer(paste0('thr', 1:nrow(thr)), 1:ncol(thr), paste, sep='_')
dth <- getThresholdMask(model, colnames(thr), subname)
v <- c(vech(cov), mns[!is.na(mns)], thr[dth])
names(v) <- c(covNames, mnsNames[!is.na(mns)], thrNames[dth])
}
return(v)
})
setMethod("genericGetExpectedStandVector", signature("BaseExpectationNormal"),
function(.Object, model, defvar.row=1, subname=model@name) {
ret <- suppressWarnings(genericGetExpected(
.Object, model, c('covariance', 'means', 'thresholds', 'slope'), defvar.row, subname))
cov <- ret[['covariance']]
mns <- ret[['means']]
if (is.null(mns)) stop("mns is null")
nv <- nrow(cov)
thr <- ret[['thresholds']]
if (prod(dim(thr)) != 0) {
thrNames <- outer(paste0('thr', 1:nrow(thr)), 1:ncol(thr), paste, sep='_')
dth <- getThresholdMask(model, colnames(thr), subname)
}
v <- .standardizeCovMeansThresholds(cov, mns, thr, ret[['slope']], dth, vector=TRUE)
return(v)
})
.standardizeCovMeansThresholds <- function(cov, means, thresholds, slope, dth, vector=FALSE){
mnames <- colnames(cov)
if (length(mnames) == 0) {
stop("I give up. Have no idea how to standardize this expectation.\nYour covariance matrix must have dimnames.")
}
if (prod(dim(thresholds)) == 0) {
ordInd <- c()
} else {
ordInd <- match(colnames(thresholds), mnames)
thresholds <- matrix( (c(thresholds) - rep(means[ordInd], each=nrow(thresholds)) ) / rep(sqrt(diag(cov)[ordInd]), each=nrow(thresholds)), nrow=nrow(thresholds), ncol=ncol(thresholds) )
means[ordInd] <- means[ordInd] - means[ordInd]
cov <- .ordinalCov2Cor(cov, ordInd)
}
if(!vector){
return(list(cov=cov, means=means, thresholds=thresholds, slope=slope))
} else {
v <- c()
vn <- c()
if (length(means)) {
for (vx in 1:length(mnames)) {
tcol <- which(vx == ordInd)
if (length(tcol) == 0) {
v <- c(v, means[vx])
vn <- c(vn, mnames[vx])
} else {
tcount <- sum(dth[,tcol])
v <- c(v, thresholds[1:tcount,tcol])
vn <- c(vn, paste0(mnames[vx], 't', 1:tcount))
}
}
}
if (length(slope)) {
v <- c(v, slope)
sname <- apply(expand.grid(row=rownames(slope),
col=colnames(slope))[,c(2,1)], 1, paste, collapse='_')
vn <- c(vn, sname)
}
for (vx in 1:length(mnames)) {
if (any(vx == ordInd)) next
v <- c(v, cov[vx,vx])
vn <- c(vn, paste0('var_', mnames[vx]))
}
v <- c(v, vechs(cov))
nv <- length(mnames)
vn <- c(vn, paste0('poly_', vechs(outer(mnames[1:nv], mnames[1:nv], FUN=paste, sep='_'))))
names(v) <- vn
return(v)
}
}
.ordinalCov2Cor <- function(cov, ordInd){
dim <- ncol(cov)
egOutCov <- matrix(0, nrow=dim, ncol=dim)
stddev <- sqrt(diag(cov))
if(is.logical(ordInd)){notOrdInd <- !ordInd} else {notOrdInd <- -ordInd}
stddev[notOrdInd] <- 1
for(i in 1:dim) {
for(j in 1:i) {
egOutCov[i,j] = cov[i, j] / (stddev[i] * stddev[j]);
egOutCov[j,i] = egOutCov[i,j]
}
}
diag(egOutCov)[ordInd] <- 1
return(egOutCov)
}
imxGetExpectationComponent <- function(model, component, defvar.row=1, subname=model$name)
{
warnModelCreatedByOldVersion(model)
if (!is(model[[subname]], "MxModel")) {
stop(paste("Submodel", subname, "in model", model$name, "is not found"))
}
if(is.null(model[[subname]]$expectation) && (class(model[[subname]]$fitfunction) %in% "MxFitFunctionMultigroup") ){
submNames <- sapply(strsplit(model$fitfunction$groups, ".", fixed=TRUE), "[", 1)
got <- list()
for(amod in submNames){
got[[amod]] <- imxGetExpectationComponent(model, component, defvar.row, subname=amod)
}
if(component=='vector' || tolower(component)=='standvector'){got <- unlist(got)}
got
} else if (length(component) == 1 && component == 'vector') {
genericGetExpectedVector(model[[subname]]$expectation, model, defvar.row, subname)
} else if (length(component) == 1 && tolower(component) == 'standvector') {
genericGetExpectedStandVector(model[[subname]]$expectation, model, defvar.row, subname)
} else {
got <- genericGetExpected(model[[subname]]$expectation, model, component, defvar.row, subname)
if (length(got) == 1) {
got[[1]]
} else {
got
}
}
}
mxGetExpected <- imxGetExpectationComponent
sse <- function(x){sum(x^2)}
#' Estimate the Jacobian of manifest model with respect to parameters
#'
#' The manifest model excludes any latent variables or processes. For
#' RAM and LISREL models, the manifest model contains only the
#' manifest variables with free means, covariance, and thresholds.
#'
#' @details
#' The Jacobian is estimated by the central finite difference.
#'
#' If the \code{standardize} argument is TRUE, then the Jacobian is for the standardized model.
#' For Normal expectations the standardized manifest model has the covariances returned as correlations, the variances returned as ones, the means returned as zeros, and the thresholds are returned as z-scores.
#' For the thresholds the z-scores are computed by using the model-implied means and variances.
#'
#' @param model an mxModel
#' @param defvar.row which row to use for definition variables
#' @param standardize logical, whether or not to standardize the parameters
#' @return a matrix with manifests in the rows and original parameters in the columns
#' @seealso \link{mxGetExpected}
omxManifestModelByParameterJacobian <- function(model, defvar.row=1, standardize=FALSE) {
theParams <- omxGetParameters(model)
if (standardize) {
ex <- 'expectation'
if (is.null(model$expectation) && is(model$fitfunction, "MxFitFunctionMultigroup")) {
submNames <- sapply(strsplit(model$fitfunction$groups, ".", fixed=TRUE), "[", 1)
ex <- paste0(submNames, ".expectation")
}
tmpModel <- mxModel(model, mxComputeJacobian(defvar.row=defvar.row, of=ex))
tmpModel <- mxRun(tmpModel, silent=TRUE)
jac <- tmpModel$compute$output$jacobian
dimnames(jac) <- list(names(mxGetExpected(model, 'standVector')), names(theParams))
} else {
jac <- numDeriv::jacobian(func=.mat2param, x=theParams, method.args=list(r=2), model=model, defvar.row=defvar.row, standardize=standardize)
dimnames(jac) <- list(names(mxGetExpected(model, 'vector')), names(theParams))
}
return(jac)
}
mxCheckIdentification <- function(model, details=TRUE){
warnModelCreatedByOldVersion(model)
notAllowedFits <- c("MxFitFunctionAlgebra", "MxFitFunctionRow", "MxFitFunctionR")
if( class(model$fitfunction) %in% notAllowedFits ){
msg <- paste("Identification check is not possible for models with", omxQuotes(notAllowedFits), 'fit functions.\n', "If you have a multigroup model, use mxFitFunctionMultigroup.")
stop(msg, call.=FALSE)
}
if(imxHasDefinitionVariable(model)){
stop("Beep beep ribby ribby. I found definition variables in your model.\nI might not give you the identification answer you're looking for in this case. See ?mxCheckIdentification.")
}
eps <- 1e-17
theParams <- omxGetParameters(model)
jac <- omxManifestModelByParameterJacobian(model)
if(imxHasConstraint(model)){
# Better if there was a separate compute step to estimate the constraintJacobian
tmpModel <- mxModel(model, mxComputeSequence(list(mxComputeNumericDeriv(hessian=FALSE), mxComputeReportDeriv())))
tmpModel <- mxRun(tmpModel, silent=TRUE)
cjac <- tmpModel$output$constraintJacobian
# drop model name from constraint name
if(length(rownames(cjac))){
rownames(cjac) <- sapply(strsplit(rownames(cjac), fixed=TRUE, split=imxSeparatorChar), '[', 2)
}
} else {
cjac <- matrix(, nrow=0, ncol=length(theParams))
colnames(cjac) <- names(theParams)
}
# Concatenate via rows the model and constraint Jacobians
jac <- rbind(jac, cjac)
# Check that rank of jac == length(theParams)
rank <- qr(jac)$rank
if(rank == length(theParams)){
message("Model is locally identified")
stat <- TRUE
} else {
message("Model is not locally identified")
stat <- FALSE
}
if(details == TRUE){
jacOC <- Null(t(jac)) # Orthogonal complement of t(jac), i.e. the basis for the null space of the column space of jac
nidp <- names(theParams)[apply(jacOC, 1, sse) > eps] # non-identified free params have non-zero rows
if(length(nidp) == 0) {
nidp <- "None"
}
} else {
nidp <- "Not Requested"
}
return(list(status=stat, jacobian=jac, non_identified_parameters=nidp))
}
.mat2param <- function(x, model, defvar.row=1, standardize=FALSE){
paramNames <- names(omxGetParameters(model))
model <- omxSetParameters(model, values=x, labels=paramNames, free=TRUE)
if(!standardize){
got <- mxGetExpected(model, 'vector', defvar.row)
} else {
got <- mxGetExpected(model, 'standVector', defvar.row)
}
return(got)
}
setGeneric("genericGenerateData",
function(.Object, model, nrows, subname, empirical, returnModel, use.miss,
.backend, nrowsProportion, silent) {
return(standardGeneric("genericGenerateData"))
})
setMethod("genericGenerateData", signature("MxExpectationNormal"),
function(.Object, model, nrows, subname, empirical, returnModel, use.miss,
.backend, nrowsProportion, silent) {
return(generateNormalData(model, nrows, subname, empirical, returnModel,
use.miss, .backend, nrowsProportion, silent))
})
.rmvnorm <- function(nrow, mean, sigma, empirical) {
if (!empirical) {
mvtnorm::rmvnorm(nrow, mean, sigma)
} else {
MASS::mvrnorm(nrow, mu=mean, Sigma=sigma, empirical=TRUE)
}
}
calcNumRows <- function(nrows, nrowsProportion, origRows, subname) {
if (is.null(nrows)) {
if (is.null(origRows)) stop("You must specify nrows")
if (!is.null(nrowsProportion)) {
nrows <- round(origRows * nrowsProportion)
if (nrows < 1) {
nrows <- 1
warning(paste("In model", omxQuotes(subname),
"nrowsProportion is too small. Rounding up to 1 row"))
}
} else {
nrows <- origRows
}
}
nrows
}
generateNormalData <- function(model, nrows, subname, empirical, returnModel, use.miss,
.backend, nrowsProportion, silent) {
origData <- findDataForSubmodel(model, subname)
origRows <- if (!is.null(origData)) { nrowMxData(origData) } else { NULL }
nrows <- calcNumRows(nrows, nrowsProportion, origRows, subname)
# Check for definition variables
if(imxHasDefinitionVariable(model[[subname]])){
if (is.null(origData)) {
stop("Definition variable(s) found, but no data are available")
}
if (origData$type != 'raw') {
stop(paste("Definition variable(s) found, but original data is type",
omxQuotes(origData$type)))
}
#Trying to use .rmvnorm() with empirical=TRUE and nrow=1 leads to an R error raised in MASS::mvrnorm():
if(empirical){
stop("argument 'empirical=TRUE' to mxGenerateData() is not compatible with definition variables")
}
origData <- origData$observed
if(nrows != nrow(origData)){
stop("Definition variable(s) found, but the number of rows in the data do not match the number of rows requested for data generation.")
}
# Generate data row by row
theCov <- imxGetExpectationComponent(model, "covariance", subname=subname)
data <- matrix(NA, nrow=nrows, ncol=ncol(theCov))
colnames(data) <- colnames(theCov)
data <- as.data.frame(data)
prevProgressLen <- 0L
for(i in 1:nrows){
theMeans <- imxGetExpectationComponent(model, "means", defvar.row=i, subname=subname)
theCov <- imxGetExpectationComponent(model, "covariance", defvar.row=i, subname=subname)
theThresh <- imxGetExpectationComponent(model, "thresholds", defvar.row=i, subname=subname)
data[i,] <- .rmvnorm(1, theMeans, theCov, empirical)
info <- paste0(subname,' ', i, '/', nrows)
if (!silent) imxReportProgress(info, prevProgressLen)
prevProgressLen <- nchar(info)
}
if (!silent) imxReportProgress('', prevProgressLen)
data <- ordinalizeDataHelper(data, theThresh, origData)
if (!is.null(origData)) {
for (dcol in setdiff(colnames(origData), colnames(data))) {
data[[dcol]] <- origData[[dcol]]
}
}
} else if (!is.null(origData) && origData$type == 'cov' && returnModel) {
theMeans <- imxGetExpectationComponent(model, "means", subname=subname)
theCov <- imxGetExpectationComponent(model, "covariance", subname=subname)
if (length(theMeans) == 0) theMeans <- NA
if (empirical) {
got <- mxModel(model[[subname]], mxData(theCov, "cov", means=theMeans, numObs=nrows))
} else {
newCov <- rWishart(1, nrows, theCov)[,,1] / nrows
dimnames(newCov) <- dimnames(theCov)
got <- mxModel(model[[subname]],
mxData(newCov, "cov", means=theMeans, numObs=nrows))
}
return(got)
} else{
theMeans <- imxGetExpectationComponent(model, "means", subname=subname)
theCov <- imxGetExpectationComponent(model, "covariance", subname=subname)
theThresh <- imxGetExpectationComponent(model, "thresholds", subname=subname)
if (length(theMeans) == 0) {
theMeans <- rep(0, nrow(theCov))
}
data <- .rmvnorm(nrows, theMeans, theCov, empirical)
colnames(data) <- colnames(theCov)
data <- as.data.frame(data)
data <- ordinalizeDataHelper(data, theThresh, origData)
}
if (use.miss && !is.null(origData) && origData$type == 'raw' &&
all(colnames(data) %in% colnames(origData$observed))) {
del <- is.na(origData$observed[,colnames(data),drop=FALSE])
if (nrows != origRows) {
del <- del[sample.int(origRows, nrows, replace=TRUE),,drop=FALSE]
}
data[del] <- NA
}
if (returnModel) {
mxModel(model[[subname]], mxData(as.data.frame(data), "raw"))
} else {
as.data.frame(data)
}
}
ordinalizeDataHelper <- function(data, thresh, origData=NULL) {
if (!is.null(origData) && is(origData, "MxData")) origData <- origData$observed
if( prod(dim(thresh)) != 0){
ordvars <- colnames(thresh)
if (length(ordvars) != ncol(thresh)) stop("thresholds missing column names")
for(avar in ordvars){
delthr <- thresh[,avar]
usethr <- 1:sum(!is.na(delthr)) # assumes NA indicates unused threshold
if (!is.null(origData) && !is.null(origData[[avar]])) {
usethr <- 1:(length(levels(origData[[avar]])) - 1L)
}
delthr <- delthr[usethr]
levthr <- 1L:(length(usethr)+1L)
if (!is.null(origData) && !is.null(origData[[avar]])) {
levthr <- levels(origData[[avar]])
}
delvar <- cut(as.vector(data[,avar]), c(-Inf, delthr, Inf), labels=levthr)
data[,avar] <- mxFactor(delvar, levels=levthr)
}
}
return(data)
}
generateRelationalData <- function(model, returnModel, .backend, subname, empirical) {
model <- model[[subname]]
# TODO add outside data reference to relational data
if (.backend) {
if (empirical) stop("mxGenerateData(..., empirical=TRUE) is not implemented for multilevel models")
plan <- mxComputeGenerateData()
model$expectation$.maxDebugGroups <- 0L
model$expectation$.optimizeMean <- 0L
modelE <- mxModel(model, plan)
modelE <- mxRun(modelE, silent=TRUE)
simData <- modelE$compute$output
names(simData) <- substr(names(simData), 1, nchar(names(simData))-5) #strip .data
for (modelName in names(simData)) {
if (modelName == model$name) {
dat <- model@data@observed
} else {
dat <- model[[modelName]]@data@observed
}
newDat <- simData[[modelName]]
toCopy <- setdiff(colnames(dat), colnames(newDat))
for (col in toCopy) newDat[[col]] <- dat[[col]]
simData[[modelName]] <- newDat
if (modelName == model$name) {
model@data@observed <- newDat
} else {
model[[modelName]]@data@observed <- newDat
}
}
if (!returnModel) {
return(simData)
} else {
return(model)
}
}
# nocov start
plan <- mxComputeSequence(list(
mxComputeOnce('expectation', 'distribution', 'flat'),
mxComputeReportExpectation()
))
modelE <- mxModel(model, plan)
modelE$expectation$.maxDebugGroups <- 100L
modelE$expectation$.rampartCycleLimit <- 0L
modelE <- mxRun(modelE, silent=TRUE)
dataEnv <- new.env()
datasets <- collectDatasets(modelE, imxGenerateNamespace(model))
for (dName in names(datasets)) {
modelName <- substr(dName, 1, nchar(dName)-5) # remove .data
assign(modelName, datasets[[ dName ]]@observed, envir=dataEnv)
}
ed <- modelE$expectation$debug
layout <- ed$layout
fmt <- paste0('g%0', ceiling(log10(ed$numGroups)), 'd')
for (gx in 1:ed$numGroups) {
groupName <- sprintf(fmt, gx)
clumpSize <- ed[[groupName]]$clumpSize
numCopies <- nrow(layout[layout$group == gx,]) %/% clumpSize
cxLength <- length(ed[[groupName]]$mean) %/% numCopies
groupTodo <- ed$layout[ed[[groupName]]$layout[,'aIndex'],]
for (cx in 1:numCopies) {
todo <- groupTodo[seq(1+(cx-1)*clumpSize, cx*clumpSize),]
repl1 <- .rmvnorm(1, ed[[groupName]]$mean[seq(1+(cx-1)*cxLength, cx*cxLength)],
sigma=as.matrix(ed[[groupName]]$covariance), empirical)
dx <- 1
for (tx in 1:nrow(todo)) {
modelName <- as.character(todo[tx,'model'])
if (modelName == modelE$name) {
submodel <- modelE
} else {
submodel <- modelE[[modelName]]
}
row <- todo[tx,'row']
manifests <- rownames(submodel$F)
beforeData <- dataEnv[[modelName]][row, manifests]
notMissing <- !is.na(beforeData)
if (sum(notMissing) > 0) {
afterData <- repl1[seq(dx, dx+sum(notMissing) - 1)]
dataEnv[[modelName]][row, manifests[notMissing] ] <- afterData
dx <- dx + sum(notMissing)
}
}
}
}
if (!returnModel) {
ret <- list()
for (n in names(dataEnv)) {
ret[[n]] <- dataEnv[[n]]
}
ret
} else {
for (modelName in names(dataEnv)) {
if (modelName == model$name) {
model@data@observed <- dataEnv[[modelName]]
} else {
model[[modelName]]@data@observed <- dataEnv[[modelName]]
}
}
model
}
# nocov end
}
simulate.MxModel <- function(object, nsim = 1, seed = NULL, ...) {
if (!is.null(seed)) {
set.seed(seed)
}
mxGenerateData(object, nsim)
}
extractData <- function(model) {
datasets <- list()
if (!is.null(model@data)) datasets <- c(datasets, list(model@data))
if (length(model@submodels)) {
datasets <- c(datasets, unlist(lapply(model@submodels, extractData), recursive=FALSE))
}
return(datasets)
}
extractObservedData <- function(model) {
lapply(extractData(model), function(mxd) {
if (mxd@type != 'raw') stop(paste("Cannot extract observed data when type=", mxd@type))
mxd@observed
})
}
mxGenerateData <- function(model, nrows=NULL, returnModel=FALSE, use.miss = TRUE,
..., .backend=TRUE, subname=NULL, empirical=FALSE, nrowsProportion=NULL, silent=FALSE) {
prohibitDotdotdot(list(...))
if (!is.null(nrows) && !is.null(nrowsProportion)) {
stop("You cannot specify both nrows and nrowsProportion")
}
if (!is.null(nrows) && nrows != round(nrows)) {
stop(paste("Cannot generate a non-integral number of rows:", nrows))
}
if (is(model, 'data.frame')) {
fake <- omxAugmentDataWithWLSSummary(mxData(model,'raw'), "ULS", "marginals",
fullWeight=FALSE, returnModel=TRUE)
obsStats <- fake$data$observedStats
fake$S$values <- obsStats$cov
fake$M$values <- obsStats$means
if (!is.null(obsStats$thresholds)) fake$thresh$values <- obsStats$thresholds
if (!is.null(nrowsProportion)) {
nrows <- round(nrow(model) * nrowsProportion)
if (nrows < 1) {
nrows <- 1
warning(paste("In model", omxQuotes(model$name),
"nrowsProportion is too small. Rounding up to 1 row"))
}
}
if(is.null(nrows)) nrows <- nrow(model)
return(mxGenerateData(fake, nrows, returnModel, empirical=empirical))
}
warnModelCreatedByOldVersion(model)
if(is.null(subname)) subname <- model$name
if (is.null(model[[subname]]$expectation) && is(model[[subname]]$fitfunction, 'MxFitFunctionMultigroup')) {
todo <- sub(".fitfunction", "", model[[subname]]$fitfunction$groups, fixed=TRUE)
for (s1 in todo) {
model <- mxModel(model, mxGenerateData(model, returnModel=TRUE, nrows=nrows,
use.miss=use.miss, .backend=.backend, subname=s1,
empirical=empirical, nrowsProportion=nrowsProportion))
}
if (!returnModel) {
return(extractObservedData(model))
} else {
return(model)
}
}
genericGenerateData(model[[subname]]$expectation, model, nrows, subname, empirical,
returnModel, use.miss, .backend, nrowsProportion, silent)
}
verifyExpectedObservedNames <- function(data, covName, flatModel, modelname, objectiveName) {
covariance <- flatModel[[covName]]
if (is(covariance, "MxMatrix") && !identical(dim(covariance), dim(data))) {
msg <- paste("The dimensions for the expected covariance matrix",
"and the observed covariance matrix",
"in the", objectiveName, "expectation function in model",
omxQuotes(modelname), "are not identical.")
stop(msg, call. = FALSE)
}
if (!identical(dimnames(covariance), dimnames(data))) {
msg <- paste0("The dimnames for the expected covariance matrix (", omxQuotes(rownames(covariance)),
") and the observed covariance matrix (", omxQuotes(rownames(data)),
") in the ", objectiveName, " expectation function in model ",
omxQuotes(modelname), " are not identical.")
stop(msg, call. = FALSE)
}
}
verifyMeans <- function(meansName, mxDataObject, flatModel, modelname) {
means <- flatModel[[meansName]]
if (!is.null(means)) {
if(any(is.na(mxDataObject@means))) {
msg <- paste("In model", omxQuotes(modelname),
"the Normal expectation function contains an expected means",
"vector but the model is missing some data",
"for the observed means.")
stop(msg, call. = FALSE)
}
meanDimnames <- dimnames(means)
}
}
setMethod("genericExpFunConvert", "MxExpectationNormal",
function(.Object, flatModel, model, labelsData, dependencies) {
.Object <- callNextMethod()
modelname <- imxReverseIdentifier(model, .Object@name)[[1]]
name <- .Object@name
if(is.na(.Object@data)) {
msg <- paste("The normal expectation function",
"does not have a dataset associated with it in model",
omxQuotes(modelname))
stop(msg, call.=FALSE)
}
mxDataObject <- flatModel@datasets[[.Object@data]]
dataName <- .Object@data
threshName <- .Object@thresholds
covName <- .Object@covariance
covariance <- flatModel[[covName]]
.Object@covariance <- imxLocateIndex(flatModel, .Object@covariance, name)
meansName <- .Object@means
.Object@means <- imxLocateIndex(flatModel, .Object@means, name)
if (inherits(mxDataObject, "MxDataDynamic")) return(.Object)
if (!(mxDataObject@type %in% c("none","raw"))) {
verifyExpectedObservedNames(mxDataObject@observed, covName, flatModel, modelname, "Normal")
verifyMeans(meansName, mxDataObject, flatModel, modelname)
}
verifyObservedNames(mxDataObject@observed, mxDataObject@means, mxDataObject@type, flatModel, modelname, "Normal")
checkNumericData(mxDataObject)
covNames <- colnames(covariance)
verifyMvnNames(covName, meansName, "expected", flatModel, modelname, class(.Object))
.Object@dataColumnNames <- covNames
.Object@dataColumns <- generateDataColumns(flatModel, covNames, dataName)
verifyThresholds(flatModel, model, labelsData, dataName, covNames, threshName)
if (single.na(.Object@dims)) {
.Object@dims <- covNames
}
if (.hasSlot(.Object, 'selectionVector')) {
.Object@selectionPlan <- prepSelectionPlan(.Object@selectionPlan, covNames)
}
return(.Object)
})
prepSelectionPlan <- function(selPlan, covNames) {
if (is.null(selPlan)) return(NULL)
ind <- match(selPlan$from, covNames)
if (any(is.na(ind))) {
stop(paste("Can't find selectionPlan from variables",
omxQuotes(selPlan$from[is.na(ind)])))
}
selPlan$from <- ind - 1L
ind <- match(selPlan$to, covNames)
if (any(is.na(ind))) {
stop(paste("Can't find selectionPlan to variables",
omxQuotes(selPlan$to[is.na(ind)])))
}
selPlan$to <- ind - 1L
selPlan
}
verifyMvnNames <- function(covName, meansName, type, flatModel, modelname, expectationName) {
if (is.na(meansName)) {
means <- NA
} else {
means <- flatModel[[meansName]]
}
covariance <- flatModel[[covName]]
if (length(covariance)) {
covDimnames <- dimnames(covariance)
if (is.null(covDimnames)) {
msg <- paste("The",type,"covariance matrix associated",
"with", expectationName, "in model",
omxQuotes(modelname), "does not contain dimnames.")
stop(msg, call. = FALSE)
}
covRows <- covDimnames[[1]]
covCols <- covDimnames[[2]]
if (is.null(covRows) || is.null(covCols) ||
(length(covRows) != length(covCols)) || !all(covRows == covCols)) {
msg <- paste("The",type,"covariance matrix associated",
"with", expectationName, "in model",
omxQuotes(modelname), "does not contain identical",
"row and column dimnames.")
stop(msg, call.=FALSE)
}
}
if (is.null(means) || (!isS4(means) && is.na(means)) || !length(means)) return()
meanDimnames <- dimnames(means)
if (is.null(meanDimnames)) {
msg <- paste("The",type,"means matrix associated",
"with", expectationName, "in model",
omxQuotes(modelname), "does not contain dimnames.")
stop(msg, call.=FALSE)
}
meanRows <- meanDimnames[[1]]
meanCols <- meanDimnames[[2]]
meanNames <- c()
if (length(meanRows) > length(meanCols)) { meanNames <- meanRows } else { meanNames <- meanCols }
if ((length(covCols) != length(meanNames)) || !all(covCols == meanNames)) {
msg <- paste("The",type,"covariance and",type,
"means matrices associated",
"with", expectationName, "in model",
omxQuotes(modelname), "do not contain identical",
"dimnames.")
stop(msg, call.=FALSE)
}
}
verifyObservedNames <- function(data, means, type, flatModel, modelname, expectationName) {
if (is.null(data)) return()
dataNames <- dimnames(data)
if(is.null(dataNames)) {
msg <- paste("The observed data associated with the",
expectationName, "expectation function in model",
omxQuotes(modelname), "does not contain dimnames.")
stop(msg, call. = FALSE)
}
if (type == "cov" || type == "cor") {
if (length(dataNames) < 2 ||
is.null(dataNames[[1]]) || is.null(dataNames[[2]]) ||
!identical(dataNames[[1]], dataNames[[2]])) {
msg <- paste("The dataset associated with the", expectationName,
"expectation function in model", omxQuotes(modelname),
"does not contain identical row and column non-NULL dimnames.")
stop(msg, call. = FALSE)
}
if (!single.na(means) && is.null(dimnames(means))) {
msg <- paste("In model", omxQuotes(modelname),
", the observed means vector does not contain column names.",
"Use the names() function to assign names to the means vector.")
stop(msg, call. = FALSE)
}
if (!single.na(means) && !identical(dataNames[[1]], dimnames(means)[[2]])) {
msg <- paste("The observed covariance or correlation matrix associated with the", expectationName,
"expectation function in model", omxQuotes(modelname),
"does not contain identical dimnames to the observed means vector.")
stop(msg, call. = FALSE)
}
} else if ((type == "raw") && (length(dataNames) < 2 || is.null(dataNames[[2]]))) {
msg <- paste("The dataset associated with the", expectationName,
"expectation function in model", omxQuotes(modelname),
"does not contain column names (use dimnames).")
stop(msg, call. = FALSE)
}
}
generateDataColumns <- function(flatModel, covNames, dataName) {
retval <- c()
if (length(covNames) == 0) return(retval)
mxd <- flatModel@datasets[[dataName]]
if (!is.null(mxd@observed)) {
dataColumnNames <- colnames(mxd@observed)
} else if (!is.null(mxd$observedStats[['cov']])) {
dataColumnNames <- colnames(mxd$observedStats[['cov']])
} else {
msg <- paste(omxQuotes(dataName), "is type='raw' but",
"raw observed data seems to be unavailable")
stop(msg, call.=FALSE)
}
retval <- match(covNames, dataColumnNames)
if (any(is.na(retval))) {
msg <- paste("The column name(s)", omxQuotes(covNames[is.na(retval)]),
"in the expected covariance matrix",
"of the expectation function",
"cannot be found in the column names of",
omxQuotes(dataName))
stop(msg, call. = FALSE)
}
return(retval - 1L)
}
updateThresholdDimnames <- function(flatExpectation, flatModel, labelsData) {
threshName <- flatExpectation@thresholds
if (is.na(threshName)) {
return(flatModel)
}
thresholds <- flatModel[[threshName]]
if (is.null(thresholds)) {
modelname <- getModelName(flatExpectation)
stop(paste("Unknown thresholds name",
omxQuotes(simplifyName(threshName, modelname)),
"detected in the expectation function",
"of model", omxQuotes(modelname)), call. = FALSE)
}
dims <- flatExpectation@threshnames
if (!is.null(colnames(thresholds)) && !single.na(dims) &&
!identical(colnames(thresholds), dims)) {
modelname <- getModelName(flatExpectation)
msg <- paste("The thresholds matrix associated",
"with the expectation function in model",
omxQuotes(modelname), "contains column names and",
"the expectation function has specified non-identical threshnames.")
stop(msg, call.=FALSE)
}
if (is.null(colnames(thresholds)) && !single.na(dims)) {
if (!flatModel@unsafe) {
tuple <- evaluateMxObject(threshName, flatModel, labelsData, new.env(parent = emptyenv()))
threshMatrix <- tuple[[1]]
if (ncol(threshMatrix) != length(dims)) {
modelname <- getModelName(flatExpectation)
msg <- paste("The thresholds matrix associated",
"with the expectation function in model",
omxQuotes(modelname), "is not of the same length as the 'threshnames'",
"argument provided by the expectation function. The 'threshnames' argument is",
"of length", length(dims), "and the expected covariance matrix",
"has", ncol(threshMatrix), "columns.")
stop(msg, call.=FALSE)
}
}
dimnames(flatModel[[threshName]]) <- list(NULL, dims)
}
return(flatModel)
}
updateExpectationDimnames <- function(flatExpectation, flatModel,
labelsData) {
unsafe <- flatModel@unsafe
covName <- flatExpectation@covariance
meansName <- flatExpectation@means
if (is.na(meansName)) {
means <- NA
} else {
means <- flatModel[[meansName]]
}
covariance <- flatModel[[covName]]
if (is.null(covariance)) {
modelname <- getModelName(flatExpectation)
stop(paste("Unknown expected covariance name",
omxQuotes(simplifyName(covName, modelname)),
"detected in the expectation function",
"of model", omxQuotes(modelname)), call. = FALSE)
}
if (is.null(means)) {
modelname <- getModelName(flatExpectation)
stop(paste("Unknown expected means name",
omxQuotes(simplifyName(meansName, modelname)),
"detected in the expectation function",
"of model", omxQuotes(modelname)), call. = FALSE)
}
dims <- flatExpectation@dims
if (!is.null(dimnames(covariance)) && !single.na(dims) &&
!identical(dimnames(covariance), list(dims, dims))) {
modelname <- getModelName(flatExpectation)
msg <- paste("The expected covariance matrix associated",
"with the expectation function in model",
omxQuotes(modelname), "contains dimnames: ",
paste(toString(dimnames(covariance)), ".", sep = ""),
"The expectation function has specified dimnames:",
paste(toString(dims), ".", sep =""))
stop(msg, call.=FALSE)
}
if (is.null(dimnames(covariance)) && !single.na(dims)) {
if (!unsafe) {
tuple <- evaluateMxObject(covName, flatModel, labelsData, new.env(parent = emptyenv()))
covMatrix <- tuple[[1]]
if (nrow(covMatrix) != ncol(covMatrix)) {
modelname <- getModelName(flatExpectation)
msg <- paste("The expected covariance matrix associated",
"with the expectation function in model",
omxQuotes(modelname), "is not a square matrix.")
stop(msg, call.=FALSE)
}
if (nrow(covMatrix) != length(dims)) {
modelname <- getModelName(flatExpectation)
msg <- paste("The expected covariance matrix associated",
"with the expectation function in model",
omxQuotes(modelname), "is not of the same length as the 'dimnames'",
"argument provided by the expectation function. The 'dimnames' argument is",
"of length", length(dims), "and the expected covariance matrix",
"has", nrow(covMatrix), "rows and columns.")
stop(msg, call.=FALSE)
}
}
dimnames(flatModel[[covName]]) <- list(dims, dims)
}
if (!isS4(means) && is.na(means)) {
return(flatModel)
}
if (!is.null(dimnames(means)) && !single.na(dims) &&
!identical(colnames(means), dims)) {
modelname <- getModelName(flatExpectation)
msg <- paste("The expected means matrix associated",
"with the expectation function in model",
omxQuotes(modelname), "contains colnames: ",
paste(toString(colnames(means)), ".", sep = ""),
"The expectation function has specified colnames:",
paste(toString(dims), ".", sep =""))
stop(msg, call.=FALSE)
}
if (is.null(dimnames(means)) && !single.na(dims)) {
if (!unsafe) {
tuple <- evaluateMxObject(meansName, flatModel, labelsData, new.env(parent = emptyenv()))
meansMatrix <- tuple[[1]]
if (nrow(meansMatrix) != 1) {
modelname <- getModelName(flatExpectation)
msg <- paste("The expected means vector associated",
"with the expectation function in model",
omxQuotes(modelname), "is not a 1 x n matrix.",
"It has dimensions", nrow(meansMatrix), "x",
paste(ncol(meansMatrix), '.', sep=''))
stop(msg, call.=FALSE)
}
if (ncol(meansMatrix) != length(dims)) {
modelname <- getModelName(flatExpectation)
msg <- paste("The expected means vector associated",
"with the expectation function in model",
omxQuotes(modelname), "is not of the same length as the 'dimnames'",
"argument provided by the expectation function. The 'dimnames' argument is",
"of length", length(dims), "and the expected means vector",
"has", ncol(meansMatrix), "columns.")
stop(msg, call.=FALSE)
}
}
dimnames(flatModel[[meansName]]) <- list(NULL, dims)
}
return(flatModel)
}
checkThreshnames <- function(threshnames) {
if (single.na(threshnames)) threshnames <- as.character(NA)
if (!is.vector(threshnames) || typeof(threshnames) != 'character') {
stop("'threshnames' argument is not a character vector")
}
if (length(threshnames) == 0) {
stop("'threshnames' argument cannot be an empty vector")
}
if (length(threshnames) > 1 && any(is.na(threshnames))) {
stop("NA values are not allowed for 'threshnames' vector")
}
tt <- table(threshnames)
if (any(tt > 1)) {
stop(paste("'threshnames' argument contains", omxQuotes(names(tt)[tt > 1]),
"more than once. \nIf you are having problems with Doppelgangers",
"perhaps you should check the basement for pods :)"))
}
return(threshnames)
}
mxExpectationNormal <-
function(covariance, means = NA, dimnames = NA, thresholds = NA,
threshnames = dimnames, ..., discrete = as.character(NA),
discreteSpec = NULL) {
prohibitDotdotdot(list(...))
if (missing(covariance) || typeof(covariance) != "character") {
stop("'covariance' argument is not a string (the name of the expected covariance matrix)")
}
if (!(single.na(means) || typeof(means) == "character")) {
stop("Means argument is not a string (the name of the expected means matrix)")
}
if (is.na(means)) means <- as.integer(NA)
if (single.na(thresholds)) thresholds <- as.character(NA)
if (single.na(dimnames)) dimnames <- as.character(NA)
if (!is.vector(dimnames) || typeof(dimnames) != 'character') {
stop("'dimnames' argument is not a character vector")
}
if (length(thresholds) != 1) {
stop("'thresholds' argument must be a single matrix or algebra name")
}
if (length(dimnames) == 0) {
stop("'dimnames' argument cannot be an empty vector")
}
if (length(dimnames) > 1 && any(is.na(dimnames))) {
stop("NA values are not allowed for 'dimnames' vector")
}
threshnames <- checkThreshnames(threshnames)
selectionVector <- NA_character_
selectionPlan <- NULL
return(new("MxExpectationNormal", covariance, means, dimnames, thresholds, threshnames,
discrete, discreteSpec, selectionVector, selectionPlan))
}
displayMxExpectationNormal <- function(expectation) {
cat("MxExpectationNormal", omxQuotes(expectation@name), '\n')
cat("$covariance :", omxQuotes(expectation@covariance), '\n')
cat("$means :", omxQuotes(expectation@means), '\n')
if (single.na(expectation@dims)) {
cat("$dims : NA \n")
} else {
cat("$dims :", omxQuotes(expectation@dims), '\n')
}
if (single.na(expectation@thresholds)) {
cat("$thresholds : NA \n")
} else {
cat("$thresholds :", omxQuotes(expectation@thresholds), '\n')
}
if (is.null(expectation@discreteSpec)) {
cat("$discreteSpec : NULL \n")
} else {
cat("$discreteSpec :", expectation@discreteSpec, '\n')
}
if (single.na(expectation@discrete)) {
cat("$discrete : NA \n")
} else {
cat("$discrete :", omxQuotes(expectation@discrete), '\n')
}
if (single.na(expectation@threshnames)) {
cat("$threshnames : NA \n")
} else {
cat("$threshnames :", omxQuotes(expectation@threshnames), '\n')
}
invisible(expectation)
}
setMethod("print", "MxExpectationNormal", function(x,...) {
displayMxExpectationNormal(x)
})
setMethod("show", "MxExpectationNormal", function(object) {
displayMxExpectationNormal(object)
})
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.