#
# Copyright 2007-2020 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.
#-------------------------------------------------------------------------------------
# Author: Michael D. Hunter
# Date: 2012.10.23
# Filename: mxSaturatedModel.R
# Purpose: This is a helper function for making and fitting a saturated model.
# It takes an mxModel and gives back the fitted saturated model.
# Example usage
# amod <- mxModel(blah)
# smod <- mxSaturatedModel(amod)
# summary(amod, SaturatedLikelihood=smod)
#-------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------
# Revision History
# Tue Oct 23 00:55:48 Central Daylight Time 2012 -- Michael Hunter copied file from personal work
# Wed 19 Jun 2013 12:23:45 Central Daylight Time -- Michael Hunter added ordinal support etc with aid of Mike Neale and Ryne Estabrook
#
#-------------------------------------------------------------------------------------
# 1. Add ability to fit multiple groups
# 2. If given a model whose data set has 10 variables but the model only uses 2
# variables, adjust function to make a saturated model of only the 2 used variables.
# TODO Check that the above are done reasonably correct.
# TODO Improve interaction between 1 & 2.
# Added ability to do Independence model.
# Added Fix variances of binary variables to 1.0
# via mxConstraint on filtered expected cov.
#-------------------------------------------------------------------------------------
# Saturated Model function definition
generateNormalReferenceModels <- function(modelName, obsdata, datatype, withMeans=FALSE, numObs, means=NA,
distribution, equateThresholds, weight = NULL) {
if (!is.null(weight) && !is.na(weight)) {
datasource <- mxData(observed=obsdata, type=datatype, numObs=numObs, means=means, weight=weight)
obsdata <- obsdata[,!names(obsdata) %in% weight]
numVar <- ncol(obsdata)
varnam <- colnames(obsdata)
} else{
datasource <- mxData(observed=obsdata, type=datatype, numObs=numObs, means=means)
numVar <- ncol(obsdata)
varnam <- colnames(obsdata)
}
if(is.null(varnam)) {
varnam <- paste("V", 1:numVar, sep="")
dimnames(obsdata) <- list(varnam, varnam)
}
if(datatype == "raw") {
if (is.data.frame(obsdata)) {
ordinalCols <- sapply(obsdata, is.ordered)
} else {
ordinalCols <- rep(FALSE, numVar)
}
if(!any(ordinalCols)){
sampcov <- cov(obsdata, use="pairwise.complete.obs")
startcov <- try(t(chol(sampcov)))
# if the cholesky fails, just use the diagonal elements
if("try-error" %in% class(startcov)){
startcov <- t(chol(diag(diag(sampcov), nrow=nrow(sampcov))))
}
startcov <- startcov[lower.tri(startcov, TRUE)]
indepcov <- diag(sampcov)
startmea <- colMeans(obsdata, na.rm=TRUE)
}
else {
ordnam <- names(obsdata[,ordinalCols])
ordinalLevels <- lapply(obsdata[,ordinalCols], levels)
numOrdinal <- sum(ordinalCols)
numOrdinalLevels <- sapply(ordinalLevels, length)
maxLevels <- max(numOrdinalLevels)
numThresholds <- maxLevels-1
startcov <- t(chol(diag(1, numVar)))
startcov <- startcov[lower.tri(startcov, TRUE)]
indepcov <- diag(1, numVar)
startmea <- rep(0, numVar)
}
} else {
startcov <- 0.3
indepcov <- 0.3
startmea <- 3.0
ordinalCols <- rep(FALSE, numVar)
}
# For all continuous data, use the Cholesky decomposition
# but for joint and all-ordinal use a symmetric matrix
# This allows one to "contrain" the total variance of ordinal variables
# to 1 without using mxConstraint.
if(!any(ordinalCols)){
ltCov <- mxMatrix(type="Lower", nrow=numVar, ncol=numVar,
values=startcov, free=TRUE, name="ltCov")
diag(ltCov$lbound) <- 0
satCov <- mxAlgebra(name="satCov", expression= ltCov %*% t(ltCov), dimnames=list(varnam, varnam))
} else {
ltCov <- NULL
satFre <- matrix(as.logical(diag(!ordinalCols, numVar)), numVar, numVar)
satFre[lower.tri(satFre, diag=FALSE)] <- TRUE
satCov <- mxMatrix(type="Symm", nrow=numVar, ncol=numVar,
values=startcov, free=satFre[lower.tri(satFre, diag=TRUE)], name="satCov",
dimnames=list(varnam, varnam))
}
saturatedModel <- mxModel(name=paste("Saturated", modelName),
datasource,
ltCov,
satCov,
mxExpectationNormal("satCov"),
mxFitFunctionML())
indCov <- mxMatrix(type="Diag", nrow=numVar, ncol=numVar, values=indepcov, free=!ordinalCols,
lbound=0, name="indCov", dimnames=list(varnam, varnam))
independenceModel <- mxModel(name=paste("Independence", modelName),
datasource, indCov,
mxExpectationNormal("indCov"), mxFitFunctionML())
if(datatype == "raw" || withMeans) {
saturatedModel <- mxModel(saturatedModel,
mxMatrix(nrow=1, ncol=numVar, values=startmea, free=TRUE, name="satMea", dimnames=list(NA, varnam)),
mxExpectationNormal("satCov", "satMea")
)
independenceModel <- mxModel(independenceModel,
mxMatrix(nrow=1, ncol=numVar, values=startmea, free=TRUE, name="satMea", dimnames=list(NA, varnam)),
mxExpectationNormal("indCov", "satMea")
)
if(any(ordinalCols)) {
thrdnam <- paste(rep(ordnam, each=numThresholds), 'ThrDev', 1:numThresholds, sep='')
unitLower <- mxMatrix("Lower", numThresholds, numThresholds, values=1, free=FALSE, name="unitLower")
thrdM <- rbind(numOrdinalLevels-1, numThresholds - numOrdinalLevels+1)
thrdfre <- apply(thrdM, 2, rep, x=c(TRUE, FALSE))
thresholdDeviations <- mxMatrix("Full",
name="thresholdDeviations", nrow=numThresholds, ncol=numOrdinal,
values=.2,
free = thrdfre,
lbound = rep( c(-Inf,rep(.01, (numThresholds-1))) , numOrdinal), # TODO adjust increment value
dimnames = list(c(), varnam[ordinalCols]),
)
if (equateThresholds) {
thresholdDeviations$labels <- thrdnam
}
saturatedMeans <- mxMatrix(nrow=1, ncol=numVar,
values=startmea, free=c(!ordinalCols), name="satMea", dimnames=list(NA, varnam))
saturatedThresholds <- mxAlgebra(unitLower %*% thresholdDeviations, name="thresholdMatrix")
saturatedModel <- mxModel(saturatedModel,
saturatedMeans, thresholdDeviations, unitLower, saturatedThresholds,
mxExpectationNormal("satCov", "satMea", thresholds="thresholdMatrix")
)
independenceModel <- mxModel(independenceModel,
saturatedMeans, thresholdDeviations, unitLower, saturatedThresholds,
mxExpectationNormal("indCov", "satMea", thresholds="thresholdMatrix")
)
}
}
if (all(ordinalCols)) {
if (distribution == 'multinomial') {
if (any(is.na(obsdata))) {
stop(paste("Saturated model for the multinomial",
"distribution is not implemented"))
}
weights <- ordinalRowWeights(saturatedModel)
saturatedModel <- list(fit=-2 * sum(weights * log(weights / sum(weights))), df=0)
} else if (distribution == 'default') {
if (length(ordinalCols) >= 12) {
message(paste("Your model has many ordinal variables.",
"It will take a long time to estimate the multivariate normal saturated model.",
"Consider using distribution='multinomial'"))
}
} else {
stop(paste("Don't know how to build reference models for the",
distribution, "distribution"))
}
} else {
if (distribution != 'default') {
stop(paste("Don't know how to build reference models for the",
distribution, "distribution"))
}
}
return(list(Saturated=saturatedModel, Independence=independenceModel))
}
ordinalRowWeights <- function(model) {
obs <- model$data$observed
weightColumn <- model$expectation$weightColumn # old API
if (is.na(weightColumn)) weightColumn <- model$data$weight # new API
if (!is.na(weightColumn)) {
weights <- obs[weightColumn]
} else {
obs <- obs[rpf::orderCompletely(obs),]
weights <- as.numeric(rpf::tabulateRows(obs))
}
weights
}
generateIFAReferenceModels <- function(model, distribution) {
if (distribution != 'default') {
stop(paste("Don't know how to build reference models for the",
distribution, "distribution"))
}
modelName <- model@name
expectation <- model@expectation
spec <- expectation$ItemSpec
nullspec <- lapply(spec, rpf::rpf.modify, 0)
data <- model$data$observed
itemName <- expectation$item
item <- model[[itemName]]
nullitem <- mxMatrix(name="item", values=mxSimplify2Array(lapply(nullspec, rpf::rpf.rparam)), condenseSlots=FALSE)
if (is.null(item)) {
stop(paste("Cannot find matrix", omxQuotes(itemName),"in model",
omxQuotes(modelName),"to create independence model"))
}
pmap <- matrix(NA, nrow(nullitem), ncol(nullitem))
for (cx in 1:ncol(item)) {
map1 <- match(names(rpf::rpf.rparam(nullspec[[cx]])),
names(rpf::rpf.rparam(spec[[cx]])))
if (!length(map1)) next
pmap[1:length(map1),cx] <- item$labels[map1,cx]
}
nullitem$labels[,] <- pmap
ind <- mxModel(name=paste("Independence", modelName),
nullitem, model$data,
mxExpectationBA81(ItemSpec=nullspec,
qpoints = expectation$qpoints,
qwidth = expectation$qwidth),
mxFitFunctionML(),
# Only need 1 iteration, but allow 2 to avoid code BLUE warning.
mxComputeEM(estep=mxComputeOnce('expectation', 'scores'),
mstep=mxComputeSequence(list(
mxComputeNewtonRaphson(),
mxComputeOnce('expectation'))),
maxIter = 2L))
dimnames(ind$item) = list(paste('p', 1:nrow(ind$item), sep=""), colnames(item))
ind$item$free <- !is.na(ind$item$values)
weights <- ordinalRowWeights(model)
saturated <- NA
if (!any(is.na(data[1,]))) { # Not sure how to handle missingness
saturated <- -2 * sum(weights * log(weights / sum(weights)))
}
return(list(Saturated=list(fit=saturated, df=0),
Independence=ind))
}
ReferenceModelHelper <- function(x, distribution, equateThresholds) {
if ( (!(isS4(x) && is(x, "MxModel"))) && !is.data.frame(x) && !(is.matrix(x) && is.numeric(x)) ) {
stop("The 'x' argument must be (1) an MxModel object, (2) a raw data frame, or (3) a raw data matrix.")
}
if ( is(x, "MxModel") ) {
if (is.null(x$fitfunction)) {
stop("Model", omxQuotes(x$name), "has no fitfunction")
}
generateReferenceModels(x$fitfunction, x, distribution, equateThresholds)
} else {
obsdata <- x
if(ncol(obsdata) != nrow(obsdata)) {
datatype <- "raw"
}
else {datatype <- "cov"}
generateNormalReferenceModels("Data Model", obsdata, datatype,
distribution=distribution, equateThresholds=equateThresholds)
}
}
mxRefModels <- function(x, run=FALSE, ..., distribution="default", equateThresholds = TRUE) {
prohibitDotdotdot(list(...))
if(is(x,"MxModel")){
warnModelCreatedByOldVersion(x)
if(imxHasDefinitionVariable(x)){
warning(
"argument 'x' is an MxModel that contains definition variables, but mxRefModels() ignores definition variables, and therefore may not do what you expect")
}
if(imxIsMultilevel(x)){
warning("The right reference models for the multilevel case are not yet known.\nI made reference models for level 1.\nI hope you know what you're doing because I don't.")
}
}
models <- lapply(ReferenceModelHelper(x, distribution, equateThresholds), function(model) {
if (!isS4(model)) return(model)
model <- omxAssignFirstParameters(model)
model <- mxOption(model, "Standard Errors", "No")
model <- mxOption(model, "Calculate Hessian", "No")
if (run) {
model <- mxRun(model, silent=FALSE)
}
model
})
models
}
omxSaturatedModel <- mxRefModels # old name
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.