R/MxSE.R

In OpenMx: Extended Structural Equation Modelling

#   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: 2016-08-12 09:28:02 by mhunter
# Date: 2016-08-13 add roxygen by tbates
# Filename: MxSE.R
# Purpose: Define the mxSE() function.		
#  This function lets the user obtain standard errors for arbitrary		
#  expressions, named entities, and algebras.		
#  It is a frontend-only file that works much like mxEval.
#------------------------------------------------------------------------------


##' Compute standard errors in OpenMx
##' 
##' @description
##' This function allows you to obtain standard errors for arbitrary
##' expressions, named entities, and algebras.
##' 
##' @param x the parameter to get SEs on (reference or expression)
##' @param model the \code{\link{mxModel}} to use.
##' @param details logical. Whether to provide further details, e.g. the full
##' sampling covariance matrix of x.
##' @param cov optional matrix of covariances among the free parameters. If 
##' missing, \code{vcov(model)} is used.
##' @param forceName logical; defaults to \code{FALSE}.  Set to \code{TRUE}
##' if \code{x} is an R symbol that refers to a character string.
##' @param silent logical; defaults to \code{FALSE}.  If \code{TRUE},
##' message-printing is suppressed.
##' @param ... further named arguments passed to \code{\link{mxEval}}
##' @param defvar.row which row to load for any definition variables
##' @param data name of data from which to load definition variables
##' 
##' @details
##' x can be the name of an algebra, a bracket address, named entity
##' or arbitrary expression.
##' When the \code{details} argument is TRUE, the full
##' sampling covariance matrix of \code{x} is also returned as part of a list.
##' The square root of the diagonals of this sampling covariance matrix are
##' the standard errors.
##' 
##' When supplying the \code{cov} argument, take care that the free parameter
##' covariance matrix is given, not the information matrix.  These 
##' two are inverses of one another.
##' 
##' This function uses the delta method to compute the standard error of arbitrary
##' and possibly nonlinear functions of the free parameters.  The delta method
##' makes a first-order Taylor approximation of the nonlinear function.  The
##' nonlinear function is a map from all the free parameters to some transformed
##' subset of parameters: the linearization of this map is given by the Jacobian
##' \eqn{J}.  In equation form, the delta method computes standard errors by the following:
##' 
##' \deqn{J^T C J}
##' 
##' where \eqn{J} is the Jacobian of the nonlinear parameter transformation
##' and \eqn{C} is the covariance matrix of the free parameters (e.g., two
##' times the inverse of the Hessian of the minus two log likelihood function).
##' 
##' @return SE value(s) returned as a matrix when \code{details} is FALSE.
##' When \code{details} is TRUE, a list of the SE value(s) and the full 
##' sampling covariance matrix.
##' @seealso - \code{\link{mxCI}}
##' @references - \url{https://en.wikipedia.org/wiki/Standard_error}
##' @examples
##' library(OpenMx)
##' data(demoOneFactor)
##' # ===============================
##' # = Make and run a 1-factor CFA =
##' # ===============================
##' 
##' latents  = c("G") # the latent factor
##' manifests = names(demoOneFactor) # manifest variables to be modeled
##' # ===========================
##' # = Make and run the model! =
##' # ===========================
##' m1 <- mxModel("One Factor", type = "RAM", 
##' 	manifestVars = manifests, latentVars = latents, 
##' 	mxPath(from = latents, to = manifests, labels=paste0('lambda', 1:5)),
##' 	mxPath(from = manifests, arrows = 2),
##' 	mxPath(from = latents, arrows = 2, free = FALSE, values = 1),
##' 	mxData(cov(demoOneFactor), type = "cov", numObs = 500)
##' )
##' m1 = mxRun(m1)
##' mxSE(lambda5, model = m1)
##' mxSE(lambda1^2, model = m1)
mxSE <- function(x, model, details=FALSE, cov, forceName=FALSE, silent=FALSE, ...,
		 defvar.row=as.integer(NA), data='data'){
  warnModelCreatedByOldVersion(model)
	if(length(model@output) > 0 && missing(cov)){
		ParamsCov <- try(vcov(model))
		if(is(ParamsCov,"try-error")){
			msg <- "Model does not have a reasonable vcov matrix or standard errors."
			stop(msg)
		}
		# if(length(model@output$infoDefinite) && !single.na(model@output$infoDefinite)){
		# 	# An indefinite Hessian usually means some SEs will be NaN:
		# 	ParamsCov <- 2*solve(model@output$hessian)
		# 	dimnames(ParamsCov) <- dimnames(model@output$hessian)
		# } else {
	} else if (missing(cov)){
		stop("Model does not have output and 'cov' argument is missing.  I'm a doctor, not a bricklayer!\nWas this model run with mxRun?")
	} else {
		ParamsCov <- cov
		if(is.null(dimnames(ParamsCov))){
			if(length(paramnames) == nrow(ParamsCov)){
				dimnames(ParamsCov) <- list(paramnames, paramnames)
			}else{
				stop(paste0("dimnames of user-supplied parameter covariance matrix are null\nand the number of rows (",  nrow(ParamsCov), ") do not match the number of free parameters (", length(paramnames), ")."))
			}
		}
	}
	
	xorig <- "x" #<--Initialize as something that will always be understandable in an error message.
	isCallEtc <- any(c('call', 'language', 'MxAlgebraFormula') %in% is(match.call()$x))
	ex <- try(eval(x), silent=TRUE)
	isChar <- !('try-error' %in% is(ex)) && is.character(ex)
	if(isCallEtc && !forceName && !isChar){
		if(!silent){message('Treating first argument as an expression')}
		xalg <- mxAlgebraFromString(Reduce(paste, deparse(match.call()$x)), name='onTheFlyAlgebra')
		xorig <- Reduce(paste, deparse(match.call()$x))
		x <- "onTheFlyAlgebra"
		model <- mxModel(model, xalg)
	} else if ('character' %in% is(x) && !isCallEtc) {
		if(!silent){message('Treating first argument as a character')}
		xalg <- mxAlgebraFromString(Reduce(paste, match.call()$x), name='onTheFlyAlgebra')
		xorig <- x
		x <- "onTheFlyAlgebra"
		model <- mxModel(model, xalg)
	} else if(isChar){
		if(!silent){message('Treating first argument as an object that stores a character')}
		xalg <- mxAlgebraFromString(ex, name='onTheFlyAlgebra')
		xorig <- ex
		x <- "onTheFlyAlgebra"
		model <- mxModel(model, xalg)
	} else {
		stop("Please, sir.  'x' must be either the name of an entity in the model, or an expression for an MxAlgebra.")
	}
	
	# Get current algebra/matrix values:
	freeparams <- omxGetParameters(model)
	paramnames <- names(freeparams)
	zoutMat <- try(mxEvalByName(x, model, compute=TRUE),silent=silent)
	if(is(zoutMat, "try-error")) {
		stop(paste0("Couldn't evaluate expression ", omxQuotes(xorig), ". Might help to check if it works in mxEval.\n",
		"Recall also that elements of submodels are addressed using **dot** syntax (not $), e.g., submodelName.objectName\n",
		"For example, to refer to an object called 'bob' in submodel 'sub1', you would say 'sub1.bob'.",
		"Also, you don't need to refer to values: so top.a_cp[1,1] NOT top.a_cp.values[1,1]"))
	}
	
	covParam <- ParamsCov
	jModel <- mxModel(model, mxComputeJacobian(of=x, defvar.row=defvar.row, data=data))
	jModel <- mxRun(jModel, silent=TRUE)
	jacTrans <- jModel$compute$output$jacobian
	covSparam <- jacTrans %*% covParam %*% t(jacTrans)
	# dimnames(covSparam) <- list(rownames(zoutMat), colnames(zoutMat))
	if(any(diag(covSparam) < 0) || any(is.na(diag(covSparam)))){
		warning("Some diagonal elements of the repeated-sampling covariance matrix of the point estimates are less than zero or NA.\nI know, right? Set details=TRUE and check the 'Cov' element of this object.")
	}
	SEs <- suppressWarnings(sqrt(diag(covSparam)))
	SEsMat <- matrix(SEs, nrow = nrow(zoutMat), ncol = ncol(zoutMat))
	if(details==TRUE){
		return(list(SE=SEsMat, Cov=covSparam))
	} else{
		return(SEsMat)
	}
}