You are here

Revision of HOWTO Make an Objective Function from Sun, 10/11/2009 - 17:15

Revisions allow you to track differences between multiple versions of your content, and revert back to older versions.

Making An Objective Function

This howto provides step-by-step instructions on writing your own objective function for OpenMx. You will first want to prototype your objective function using either an MxAlgebraObjective or an MxRObjective. An algebra objective consists of an arbitrary number of matrix algebra statements that ultimately evaluate to a 1 x 1 matrix. An R objective consists of an arbitrary R function that evaluates a model to return a numeric value. Algebra objectives are computationally faster than R objectives, but algebra objectives are more limited in their expressive power. When you are satisfied prototyping your new objective function, we will guide you through the process of writing the front-end for this objective function in R, and writing the back-end for this objective function in C.

Step 1: Prototyping

Using MxAlgebraObjective

Can you express your objective function using the matrix operators and functions provided by OpenMx? If so, then create an MxAlgebraObjective. The following example is from the SimpleAlgebra demo.

A <- mxMatrix("Full", values = c(1, 2), nrow=1, ncol=2, name="A")
B <- mxMatrix("Full", values = c(3, 4), nrow = 2, ncol = 1, free = c(TRUE, FALSE), name="B")
algebra <- mxAlgebra(A %*% B - 3, "algebra")
algSquared <- mxAlgebra(algebra %*% algebra, "algebraSquared")
 
model <- mxModel("themodel")
model <- mxModel(model, A, B)
model <- mxModel(model, algebra, algSquared)
model <- mxModel(model, mxAlgebraObjective("algebraSquared"))
 
# Test 1: Algebra is a multiply and subtract
 
modelOut <- mxRun(model)

Using MxRObjective

Otherwise, express your objective function in the R front-end. Write a function in R that accepts two arguments. The first argument is the mxModel that should be evaluated, and the second argument is some persistent state information that can be stored between one iteration of optimization to the next iteration. It is valid for the function to simply ignore the second argument. The function must return either a single numeric value, or a list of exactly two elements. If the function returns a list, the first argument must be a single numeric value and the second element will be the new persistent state information to be passed into this function at the next iteration. The single numeric value will be used by the optimizer to perform optimization. For more help on the slots in the MxModel class, see the help provided with the OpenMx package (type '?MxModel' in R). The following example is from the SimpleRObjective demo.

A <- mxMatrix(nrow = 2, ncol = 2, values = c(1:4), free = TRUE, name = 'A')
 
squared <- function(x) { x ^ 2 }
 
objFunction <- function(model, state) {
    values <- mxEval(A, model) 
    return(squared(values[1,1] - 4) + squared(values[1,2] - 3) +
        squared(values[2,1] - 2) + squared(values[2,2] - 1))
}
objective <- mxRObjective(objFunction)
 
model <- mxModel('model', A, objective)
 
modelOut <- mxRun(model)

Step 2: Write Your Objective Function

Step 2a: Implement the front-end

In order to write your own objective function, you will need to create your own objective function class in R. For teaching purposes, we are going to use the MxAlgebraObjective class as a guide. First, your class needs to be a sub-class of the MxBaseObjective class. The MxBaseObjective class contains three slots: 'name' which is a character string, 'data' which is either a character string or a number, and 'result' which is a 1x1 matrix that stores the result of optimization. Your class should contain any additional slots that you need to pass along to the optimizer. In the MxAlgebraObjective, we have an additional slot that will inform the optimizer as to which MxAlgebra to use for optimization.

setClass(Class = "MxAlgebraObjective",
    representation = representation(
        algebra = "MxCharOrNumber"),
    contains = "MxBaseObjective")

Next, you will need to overwrite three generic methods for your new objective function class. These three methods are initialize, omxObjFunConvert, and omxObjFunNamespace.

initialize

The initialize method is a constructor method. It will never be directly called by the user, but it will be critical for making new instances of your class. When writing your initialize method, it is required that the default value for the data slot is as.integer(NA) and the default value for the name slot is 'objective'. In fact, the data slot will be modified only be the omxObjFunConvert method, and the name slot will only be modified by the omxObjFunNamespace method.

setMethod("initialize", "MxAlgebraObjective",
    function(.Object, algebra, 
        data = as.integer(NA), name = 'objective') {
        .Object@name <- name
        .Object@algebra <- algebra
        .Object@data <- data
        return(.Object)
    }
)

omxObjModelConvert

The omxObjModelConvert method is used in the event that you need to perform transformations on the model that are specific to your objective function. One example use-case of this method is to transform RAM objective functions to FIML objective functions, if raw data is provided to the RAM objective function. During that particular transformation, several additional matrices and algebras are added to the model, and the RAM objective function is replaced with a FIML objective function. For now, we will show you the default behavior of this function, which is the most common behavior. You do not need to copy/paste the default behavior into your objective function, it is inherited from "MxBaseObjective".The function omxObjModelConvert accepts four parameters: '.Object' is the current objective function, 'job' is the entire job that was passed to the mxRun() function, 'model' is the individual model associated with the current objective function, and 'flatJob' is a flattened version of 'job'.  This function MUST return the job.  In order to transform the current model, use the following idiom after updating the model: job[model@name] &lt;- model.

setMethod("omxObjModelConvert", "MxBaseObjective",
    function(.Object, job, model, flatJob) {
        return(job)
})

omxObjFunNamespace The omxObjFunNamespace method is used to transform an objective function from its local model namespace, into the flattened model namespace. Basically you must use the function omxConvertIdentifier on any slot that is storing name information. Remember to always invoke omxConvertIdentifier on the 'name' slot and the 'data' slot. Shown below is the implementation for the FIML objective function.

setMethod("omxObjFunNamespace", signature("MxFIMLObjective"), 
    function(.Object, modelname, namespace) {
        .Object@name <- omxIdentifier(modelname, .Object@name)
        .Object@covariance <- omxConvertIdentifier(.Object@covariance, 
            modelname, namespace)
        .Object@means <- omxConvertIdentifier(.Object@means, 
            modelname, namespace)
        .Object@data <- omxConvertIdentifier(.Object@data, 
            modelname, namespace)
        return(.Object)
})

omxObjFunConvert

The omxObjFunConvert method is invoked to transform an objective function from a form that interacts with the R front-end, to a form that is suitable for processing by the C back-end. Any character strings in the objective function are converted into matrix/algebra pointers for the back-end to process. You should read the documentation for the function omxLocateIndex to familiarise yourself with matrix/algebra pointers. Two arguments are passed into this function, a flattened version of the model, and the original model itself. When at all possible, it is recommended that you interact with the flattened version of the model. Below we have included the omxObjFunConvert method for both the algebra objective and the FIML objective. The FIML objective has been included so that you may see how the data slot is processed. NOTE: .Object@algebra is a string on input, and is converted into a matrix/algebra pointer.

setMethod("omxObjFunConvert", signature("MxAlgebraObjective"), 
    function(.Object, flatModel, model) {
        name <- .Object@name
        algebra <- .Object@algebra
        if (is.na(algebra)) {
            modelname <- omxReverseIdentifier(model, .Object@name)[[1]]
            msg <- paste("The algebra name cannot be NA",
            "for the algebra objective of model", omxQuotes(modelname))
            stop(msg, call. = FALSE)
        }
        algebraIndex <- omxLocateIndex(flatModel, algebra, name)
        .Object@algebra <- algebraIndex
        return(.Object)
})

NOTE: .Object@data, .Object@means, and .Object@covariance are all strings at the beginning of this function, and are converted into data pointers or matrix/algebra pointers.

setMethod("omxObjFunConvert", signature("MxFIMLObjective"), 
    function(.Object, flatModel, model) {
        modelname <- omxReverseIdentifier(model, .Object@name)[[1]]
        name <- .Object@name
        if(is.na(.Object@data)) {
            msg <- paste("The FIML objective",
                "does not have a dataset associated with it in model",
                omxQuotes(modelname))
            stop(msg, call.=FALSE)
        }
        mxDataObject <- flatModel@datasets[[.Object@data]]
        if (mxDataObject@type != 'raw') {
            msg <- paste("The dataset associated with the FIML objective", 
                "in model", omxQuotes(modelname), "is not raw data.")
            stop(msg, call.=FALSE)
        }
        verifyObservedNames(mxDataObject@observed, mxDataObject@type, flatModel, modelname, "FIML")
        checkNumericData(mxDataObject)
        meansName <- .Object@means
        covName <- .Object@covariance
        dataName <- .Object@data
        threshNames <- .Object@thresholds
        .Object@means <- omxLocateIndex(flatModel, .Object@means, name)
        .Object@covariance <- omxLocateIndex(flatModel, .Object@covariance, name)
        .Object@data <- omxLocateIndex(flatModel, .Object@data, name)
        verifyExpectedNames(covName, meansName, flatModel, modelname, "FIML")
        .Object@definitionVars <- generateDefinitionList(flatModel)
        .Object@dataColumns <- generateDataColumns(flatModel, covName, dataName)
        .Object@thresholdColumns <- generateThresholdColumns(flatModel, threshNames, dataName)
        if (length(mxDataObject@observed) == 0) {
            .Object@data <- as.integer(NA)
        }
        return(.Object)
})

User-interface

Finally, you should write your own method that creates instances of your new objective function class. This interface will invoke the 'initialize' constructor that you wrote earlier, but should also include some basic sanity-checking before invoking the constructor. Here is the algebra objective function interface. Since this method does not have access to any parts of the model, you will not be able to perform full error-checking at this stage. Full error-checking should be performed in the omxObjFunConvert method.

mxAlgebraObjective <- function(algebra) {
    if (missing(algebra) || typeof(algebra) != "character") {
        stop("Algebra argument is not a string (the name of the algebra)")
    }
    return(new("MxAlgebraObjective", algebra))
}

mxRun() Workflow

The following is the workflow of the mxRun() function, showing where each objective function method occurs in the workflow.

  1. check namespace for errors
  2. translate objectives (omxObjModelConvert)
  3. run independent sub-models
  4. flatten model (omxObjFunNamespace)
  5. check algebras for errors
  6. convert global variables and constants in algebra expressions
  7. convert model pieces to back-end (omxObjFunConvert)
  8. invoke the back-end
  9. repopulate the model results from the optimizer

Step 2(b): Implement the back-end

TODO: write this section