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.
#------------------------------------------------------------------------------
# Revision History
# Fri 08 Mar 2013 15:35:29 Central Standard Time -- Michael Hunter copied file from
# OpenMx\branches\dependency-tracking\R\WLSObjective.R
# I think Tim Brick wrote the original.
# Michael Hunter edited the file to use fitfunctions
#
#------------------------------------------------------------------------------
# **DONE**
setClass(Class = "MxFitFunctionWLS",
contains = "MxBaseFitFunction",
representation = representation(
type = "character",
continuousType = "character",
fullWeight = "logical"
))
# **DONE**
setMethod("initialize", "MxFitFunctionWLS",
function(.Object, ...) {
.Object <- callNextMethod()
.Object@type <- ..1
.Object@continuousType <- ..2
.Object@fullWeight <- ..3
.Object
}
)
# **DONE**
setMethod("qualifyNames", signature("MxFitFunctionWLS"),
function(.Object, modelname, namespace) {
.Object@name <- imxIdentifier(modelname, .Object@name)
return(.Object)
})
# **DONE**
setMethod("genericFitConvertEntities", "MxFitFunctionWLS",
function(.Object, flatModel, namespace, labelsData) {
name <- .Object@name
modelname <- imxReverseIdentifier(flatModel, .Object@name)[[1]]
expectName <- paste(modelname, "expectation", sep=".")
expectation <- flatModel@expectations[[expectName]]
dataname <- expectation@data
if (flatModel@datasets[[dataname]]@type != 'raw') {
if (.Object@vector) {
modelname <- getModelName(.Object)
msg <- paste("The WLS fit function",
"in model", omxQuotes(modelname),
"does not have raw data")
stop(msg, call.=FALSE)
}
}
return(flatModel)
})
# **DONE**
setMethod("genericFitFunConvert", "MxFitFunctionWLS",
function(.Object, flatModel, model, labelsData, dependencies) {
.Object <- callNextMethod()
name <- .Object@name
modelname <- imxReverseIdentifier(model, .Object@name)[[1]]
expectName <- paste(modelname, "expectation", sep=".")
if (expectName %in% names(flatModel@expectations)) {
expectIndex <- imxLocateIndex(flatModel, expectName, name)
} else {
expectIndex <- as.integer(NA)
}
.Object@expectation <- expectIndex
eobj <- flatModel@expectations[[1L+expectIndex]]
return(.Object)
})
# **DONE**
setMethod("genericFitInitialMatrix", "MxFitFunctionWLS",
function(.Object, flatModel) {
flatFitFunction <- flatModel@fitfunctions[[.Object@name]]
modelname <- imxReverseIdentifier(flatModel, flatFitFunction@name)[[1]]
expectationName <- paste(modelname, "expectation", sep = ".")
expectation <- flatModel@expectations[[expectationName]]
if (is.null(expectation)) {
msg <- paste("The WLS fit function",
"has a missing expectation in the model",
omxQuotes(modelname))
stop(msg, call.=FALSE)
}
if (is.na(expectation@data)) {
msg <- paste("The WLS fit function has",
"an expectation function with no data in the model",
omxQuotes(modelname))
stop(msg, call.=FALSE)
}
mxDataObject <- flatModel@datasets[[expectation@data]]
if (!(mxDataObject@type %in% c('raw','none'))) {
msg <- paste("The dataset associated with the WLS expectation function",
"in model", omxQuotes(modelname), "is not raw data.")
stop(msg, call.=FALSE)
}
cols <- ncol(mxDataObject@observed)
return(matrix(as.double(NA), cols, cols))
})
# **DONE**
mxFitFunctionWLS <- function(type=c('WLS','DWLS','ULS'),
allContinuousMethod=c("cumulants", "marginals"),
fullWeight=TRUE) {
type <- match.arg(type)
allContinuousMethod <- match.arg(allContinuousMethod)
return(new("MxFitFunctionWLS", type, allContinuousMethod, as.logical(fullWeight)))
}
# **DONE**
displayMxFitFunctionWLS <- function(fitfunction) {
cat("MxFitFunctionWLS", omxQuotes(fitfunction@name), '\n')
if (length(fitfunction@result) == 0) {
cat("$result: (not yet computed) ")
} else {
cat("$result:\n")
}
print(fitfunction@result)
invisible(fitfunction)
}
# **DONE**
setMethod("print", "MxFitFunctionWLS", function(x, ...) {
displayMxFitFunctionWLS(x)
})
setMethod("show", "MxFitFunctionWLS", function(object) {
displayMxFitFunctionWLS(object)
})
# deprecated
# nocov start
imxWlsStandardErrors <- function(model){
#TODO add safety check
# Is it a WLS fit function
# Does the data have @fullWeight
isMultiGroupModel <- is.null(model$expectation) && (class(model$fitfunction) %in% "MxFitFunctionMultigroup")
fwMsg <- "Terribly sorry, master, but you cannot compute standard errors without the full weight matrix."
theParams <- omxGetParameters(model)
if( isMultiGroupModel ){
submNames <- sapply(strsplit(model$fitfunction$groups, ".", fixed=TRUE), "[", 1)
sV <- list()
sW <- list()
sD <- c()
for(amod in submNames){
sV[[amod]] <- model[[amod]]$data$acov
fullWeight <- model[[amod]]$data$fullWeight
if(single.na(fullWeight)){stop(paste(fwMsg, '\nOffending model is', amod))}
sW[[amod]] <- MASS::ginv(fullWeight)
sD[[amod]] <- single.na(model[[amod]]$data$thresholds)
}
if( !(all(sD == TRUE) || all(sD == FALSE)) ){
stop("I feel like I'm getting mixed signals. You have some ordinal data, and some continuous data, and I'm not sure what to do. Post this on the developer forums.")
}
d <- omxManifestModelByParameterJacobian(model, standardize=!any(sD))
V <- Matrix::bdiag(sV)
W <- Matrix::bdiag(sW)
} else {
d <- omxManifestModelByParameterJacobian(model, standardize=!single.na(model$data$thresholds))
V <- model$data$acov #used weight matrix
fullWeight <- model$data$fullWeight
if(single.na(fullWeight)){stop(paste(fwMsg, '\nOffending model is', model@name))}
W <- MASS::ginv(fullWeight)
}
dvd <- try(solve( t(d) %*% V %*% d ), silent=TRUE)
if (!inherits(dvd, "try-error")) {
nacov <- as.matrix(dvd %*% t(d) %*% V %*% W %*% V %*% d %*% dvd)
wls.se <- matrix(sqrt(diag(nacov)), ncol=1)
} else {
nacov <- matrix(NA, nrow=length(theParams), ncol=length(theParams))
wls.se <- matrix(NA, nrow=length(theParams), ncol=1)
warning("Standard error matrix is not invertible.\nIs your model not identified or in need of a constraint?\nCheck for identification with mxCheckIdentification.\nReturning model with all NA standard errors.")
}
dimnames(nacov) <- list(names(theParams), names(theParams))
rownames(wls.se) <- names(theParams)
#SE is the standard errors
#Cov is the analog of the Hessian for WLS
return(list(SE=wls.se, Cov=nacov, Jac=d))
}
# nocov end
# deprecated
# nocov start
imxWlsChiSquare <- function(model, J=NA){
samp.param <- mxGetExpected(model, 'standVector')
theParams <- omxGetParameters(model)
numOrdinal <- 0
numObs <- 0
isMultiGroupModel <- is.null(model$expectation) && (class(model$fitfunction) %in% "MxFitFunctionMultigroup")
fwMsg <- "Terribly sorry, master, but you cannot compute chi square without the full weight matrix."
if( isMultiGroupModel ){
submNames <- sapply(strsplit(model$fitfunction$groups, ".", fixed=TRUE), "[", 1)
sV <- list()
sW <- list()
expd.param <- c()
sD <- c()
for(amod in submNames){
cov <- model[[amod]]$data$observed
mns <- model[[amod]]$data$means
thr <- model[[amod]]$data$thresholds
numObs <- numObs + model[[amod]]$data$numObs
sD[[amod]] <- single.na(thr)
if(!single.na(thr)){
expd.param <- c(expd.param, .standardizeCovMeansThresholds(cov, mns, thr, !is.na(thr), vector=TRUE))
numOrdinal <- numOrdinal + ncol(thr)
} else {
expd.param <- c(expd.param, cov[lower.tri(cov, TRUE)], mns[!is.na(mns)], thr[!is.na(thr)])
}
sV[[amod]] <- model[[amod]]$data$acov
fullWeight <- model[[amod]]$data$fullWeight
if(single.na(fullWeight)){stop(paste(fwMsg, '\nOffending model is', amod))}
sW[[amod]] <- MASS::ginv(fullWeight)
}
V <- Matrix::bdiag(sV)
W <- Matrix::bdiag(sW)
} else {
cov <- model$data$observed
mns <- model$data$means
thr <- model$data$thresholds
numObs <- model$data$numObs
sD <- single.na(thr)
if(!single.na(thr)){
expd.param <- .standardizeCovMeansThresholds(cov, mns, thr, !is.na(thr), vector=TRUE)
numOrdinal <- numOrdinal + ncol(thr)
} else {
expd.param <- c(cov[lower.tri(cov, TRUE)], mns[!is.na(mns)], thr[!is.na(thr)])
}
V <- model$data$acov #used weight matrix
fullWeight <- model$data$fullWeight
if(single.na(fullWeight)){stop(paste(fwMsg, '\nOffending model is', model@name))}
W <- MASS::ginv(fullWeight)
}
if( !(all(sD == TRUE) || all(sD == FALSE)) ){
stop("I'm getting some mixed signals. You have some ordinal data, and some continuous data, and I'm not sure what to do. Post this on the developer forums.")
}
e <- samp.param - expd.param
if(single.na(J)){
jac <- omxManifestModelByParameterJacobian(model, standardize=!any(sD))
} else {jac <- J}
jacOC <- Null(jac)
if(prod(dim(jacOC)) > 0){
x2 <- t(e) %*% jacOC %*% ginv( as.matrix(t(jacOC) %*% W %*% jacOC) ) %*% t(jacOC) %*% e
} else {x2 <- 0}
df <- qr(jacOC)$rank
dvd <- try(solve( t(jac) %*% V %*% jac ), silent=TRUE)
if (!inherits(dvd, "try-error")) {
U <- V - V %*% jac %*% dvd %*% t(jac) %*% V
} else {
U <- matrix(NA, nrow=nrow(W), ncol=ncol(W))
}
UW <- as.matrix(U %*% W) # For multigroup case, convert Matrix type to matrix
UW2 <- UW %*% UW # unclear if this should be UW^2 i.e. elementwise power
trUW <- sum(diag(UW))
madj <- trUW/df
x2m <- as.numeric(model$fitfunction$result)/madj
dstar <- round((trUW^2) / sum(diag(UW2)))
mvadj <- trUW^2/dstar
x2mv <- as.numeric(model$fitfunction$result)/mvadj
# N.B. x2mv is off by a factor of N where N is the total number of rows in all data sets for the ULS case.
V <- as.matrix(V)
I <- diag(1, nrow=nrow(V))
x2mv <- x2mv*ifelse(all(V[V!=0] == I[V != 0]), 1/numObs, 1)
return(list(Chi=x2, ChiDoF=df, ChiM=x2m, ChiMV=x2mv, mAdjust=madj, mvAdjust=mvadj, dstar=dstar))
}
# nocov end
approveWLSIntervals <- function(flatModel, modelName) {
ff <- flatModel@fitfunctions[[ paste0(modelName, '.fitfunction') ]]
if (is(ff, "MxFitFunctionWLS")) {
if (ff@type != 'WLS') {
stop(paste("Confidence intervals are not supported for DWLS or ULS. ",
"Try mxSE or switch", omxQuotes(modelName), "to full WLS"))
}
} else if (is(ff, "MxFitFunctionMultigroup")) {
for (g1 in flatModel@fitfunction$groups) {
approveWLSIntervals(flatModel, g1)
}
}
}
#' Determine whether a dataset will have weights and summary statistics for the means if used with mxFitFunctionWLS
#'
#' Given either a data.frame or an mxData of type raw, this function determines whether \code{mxFitFunctionWLS}
#' will generate expectations for means.
#'
#' All-continuous data processed using the "cumulants" method lack means, while
#' all continuous data processed with allContinuousMethod = "marginals" will have means.
#'
#' When data are not all continuous, allContinuousMethod is ignored, and means are modelled.
#'
#' @param data the (currently raw) data being used in a \code{\link{mxFitFunctionWLS}} model.
#' @param allContinuousMethod the method used to process data when all columns are continuous.
#' @param verbose logical. Whether to report diagnostics.
#' @return - list describing the data.
#' @family Data Functions
#' @seealso - \code{\link{mxFitFunctionWLS}}, \code{\link{omxAugmentDataWithWLSSummary}}
#' @examples
#'
#' # ====================================
#' # = All continuous, data.frame input =
#' # ====================================
#'
#' tmp = mxDescribeDataWLS(mtcars, allContinuousMethod= "cumulants", verbose = TRUE)
#' tmp$hasMeans # FALSE - no means with cumulants
#' tmp = mxDescribeDataWLS(mtcars, allContinuousMethod= "marginals")
#' tmp$hasMeans # TRUE we get means with marginals
#'
#' # ==========================
#' # = mxData object as input =
#' # ==========================
#' tmp = mxData(mtcars, type="raw")
#' mxDescribeDataWLS(tmp, allContinuousMethod= "cumulants", verbose = TRUE)$hasMeans # FALSE
#' mxDescribeDataWLS(tmp, allContinuousMethod= "marginals")$hasMeans # TRUE
#'
#' # =======================================
#' # = One var is a factor: Means modelled =
#' # =======================================
#' tmp = mtcars
#' tmp$cyl = factor(tmp$cyl)
#' mxDescribeDataWLS(tmp, allContinuousMethod= "cumulants")$hasMeans # TRUE - always has means
#' mxDescribeDataWLS(tmp, allContinuousMethod= "marginals")$hasMeans # TRUE
#'
mxDescribeDataWLS <- function(data, allContinuousMethod = c("cumulants", "marginals"), verbose=FALSE){
allContinuousMethod = match.arg(allContinuousMethod)
if(inherits(data, "data.frame")){
# all good
} else if(inherits(data, "MxDataStatic") && data$type == "raw"){
data = data$observed
}else{
message("mxDescribeDataWLS currently only knows how to process dataframes and mxData of type = 'raw'.\n",
"You offered up an object of class: ", omxQuotes(class(data)))
}
if(all(sapply(data, FUN= is.numeric))){
if(verbose){ print("all continuous") }
if(allContinuousMethod == "cumulants"){
return(list(hasMeans = FALSE))
} else {
return(list(hasMeans = TRUE))
}
}else{
# Data with any non-continuous vars have means under WLS
return(list(hasMeans = TRUE))
}
}
##' imxHasWLS
##'
##' This is an internal function exported for those people who know
##' what they are doing. This function checks if a model uses a
##' fitfunction with WLS units.
##'
##' @param model model
imxHasWLS <- function(model){
wlsUnits <- c("r'Wr", "r'wr")
if(!is.null(model@output$fitUnits)){
if(model@output$fitUnits %in% wlsUnits){return(TRUE)}
}
if(is.null(model@fitfunction)){return(FALSE)}
if(is(model@fitfunction, "MxFitFunctionWLS")){return(TRUE)}
if(length(model@fitfunction$units) && model@fitfunction$units %in% wlsUnits){return(TRUE)}
if( is(model@fitfunction, "MXFitFunctionMultigroup") ){
#Just in case the user provided 'modelName.fitfunction':
submodnames <- unlist(lapply(strsplit(model@fitfunction@groups,"[.]"),function(x){x[1]}))
for(i in 1:length(model@submodels)){
if(model@submodels[[i]]@name %in% submodnames){
probe <- imxHasWLS(model@submodels[[i]])
if(probe){return(probe)}
}
}
}
return(FALSE)
}
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.