Fri, 01/15/2010 - 19:04

Below is a function that prints to the console all of the means and covariances matrices in a model and all of the nested submodels. The function .collectMeansCovariances() places all of the matrices into a single list. The remaining lines of showMeansCovariances() are just pretty-printing, they are not interesting.

.collectMeansCovariances() extracts the means and covariances of the top model, and then recursively calls itself to collect the matrices from the submodels. The function does not assume a specific type of objective function. Instead it browses through the slots of the objective function, searching for the slot 'covariance' and the slot 'means'. If those slots are found, then a call to mxEval() must be constructed. It is important to use the as.symbol() when calling the substitute() function, otherwise the covariance name will be interpreted as a regular string instead of a symbol.

Finally, perform the recursive step and then collect the results. The call to unlist() is used because each call to .collectMeansCovariances() returns a list result.

showMeansCovariances <- function(model) { resultsList <- .collectMeansCovariances(model) if(length(resultsList) > 0) { resultsNames <- names(resultsList) for(i in 1:length(resultsList)) { cat(resultsNames[[i]],'\n') print(resultsList[[i]]) cat('\n') } } } .collectMeansCovariances <- function(model) { listReturn <- list() if (!is.null(model$objective)) { objective <- model$objective slots <- slotNames(objective) # extract the covariance if ('covariance' %in% slots) { covName <- objective@covariance covariance <- eval(substitute(mxEval(x, model), list(x = as.symbol(covName)))) listReturn[[covName]] <- covariance } # extract the means if ('means' %in% slots) { meansName <- objective@means means <- eval(substitute(mxEval(x, model), list(x = as.symbol(meansName)))) listReturn[[meansName]] <- means } } # Recursively collect means and covariances of submodels submodels <- lapply(model@submodels, .collectMeansCovariances) listReturn <- append(listReturn, unlist(submodels, FALSE)) return(listReturn) }

Here is a cleaned up version of the script. A couple of general comments:

I take your points. I really only wanted one function that shows statistics for

1. full model

2. nested model (if there is one)

3. list of nested models (if there is a list)

So I tried to move the old 'printFitStatistics1' down and renamed it to '.showFitStatistics' and renamed the core function 'showFitStatistics' but seem to have broken something in the process (see attached). Not familiar enough yet with lapply to know how to fix it.

Also 'showMatrixLabels' works, but maybe can be done more efficiently too?

Thanks!

h.

I wrote showFitStatistics to take two arguments. If the second argument is missing, then return statistics on the model. If the second argument is a model, then return comparative statistics. If the second argument is a list of models, then return all comparative statistics.

I also rewrote showMatrixLabels. The primary issue was to separate the matrix formatting functionality from the MxModel functionality. I wrote one function 'formatMatrix' that takes an R matrix and formats it nicely. Another function 'evalQuote' takes a character string and calls 'mxEval' on that string as if it were an expression. I can't figure out the correct column names for the matrices in this example, so you can complete that.

Hey, you're doing great with the matrix..

that is really cool! great work ^_^

Hi,

The attached job runs a saturated model and an ACE model with two submodels

and applies all the recently developed functions (specificationMatrices, showMeansCovariances, printMatrixLabels and printFitStatistics). The latter two may be improved but are working, except the third variety of printFitStatistics that would allow displaying statistics for a list of nested models. Any suggestions?

Hermine.

This is very cool. However, there is a slight problem with the attached script, which uses a data file at http://www.vipbg.vcu.edu/~vipbg/OpenMx/myData/usski.rec

The model is set up such that the expected mean matrix in the MZ submodel is actually a matrix from the parent level, indexed by ACE.expcovMZ. This type of specification generates the following error:

> showMeansCovariances(twinACEFit)

Error in eval(expr, envir, enclos) : object 'ACE.expCovMZ' not found

Although it does exist there:

> twinACEFit@algebras$expCovMZ

mxAlgebra 'expCovMZ'

@formula: rbind(cbind(A + C + E, A + C), cbind(A + C, A + C + E))

@result:

t1bic t1tri t2bic t2tri

t1bic 11.213490 9.484041 8.977883 7.908478

t1tri 9.484041 10.218682 7.908478 8.339780

t2bic 8.977883 7.908478 11.213490 9.484041

t2tri 7.908478 8.339780 9.484041 10.218682

dimnames:

[[1]]

[1] "t1bic" "t1tri" "t2bic" "t2tri"

[[2]]

[1] "t1bic" "t1tri" "t2bic" "t2tri"

Ah ok. Try this change. The names of the resulting list are long. The alternative is to generate a list of lists: the outer list for all the models, and the inner list for the covariance and means per model.

Works like a charm! Thank you!

I tried to extend the 'freeFixedLabels' function to automatically apply it to all mxMatrices in a model (ideally all matrices, except Iden,Zero,Unit matrices) - based on your 'showMeansCovariances' function but I don't know how to scroll over matrices as there isn't a type or name slot.

Here my feeble attempt that obviously has no chance of working:

specificationMatrix <- function(model) {

resultsList <- .collectMatrices(model)

if(length(resultsList) > 0) {

resultsNames <- names(resultsList)

for(i in 1:length(resultsList)) {

cat(resultsNames[[i]],'\n')

print(resultsList[[i]])

cat('\n')

}

}

}

.collectMatrices <- function(model) {

listReturn <- list()

if (!is.null(model@matrices)) {

# extract the matrices with potential free parameters

# if (model@matrices@type %in% c("Full","Lower","Diag","Symm","Stand")) {

# if(!is(matrix, "MxMatrix")) {

# stop("function argument must be MxMatrix object")

# }

retval <- mapply(.specificationHelper, matrix@labels, matrix@free, matrix@values)

retval <- matrix(retval, nrow(matrix), ncol(matrix))

dimnames(retval) <- dimnames(matrix)

retval <- print(retval, quote=FALSE)

storeName <- paste('model:', model@name,', matrix:', model@matrices@name, sep='')

listReturn[[storeName]] <- retval

# return(retval)

# }

}

# Recursively collect matrices of model

names(model@matrices) <- NULL

matrix <- lapply(model@matrices, .collectMatrices)

listReturn <- append(listReturn, unlist(matrices, FALSE))

return(listReturn)

}

.specificationHelper <- function(label, free, value) {

if(free) return(paste('[', label, ']', sep = ''))

else return(value)

}

Michael (Spiegel), can you point me in the right direction?

Thanks,

Hermine.

I've attached a script that has the functions for what I believe you were trying to do. When building complex functions, there's a very common pattern of repeating the same operation over a set of values. This pattern can be accomplished by using either a for-loop or using only of the apply statements (apply, lapply, sapply, mapply, etc.). You'll see an example of both a for-loop and an apply statement in the attached script.

The for-loop is simpler to understand. It is not "recommended" by experienced R programmers. The for-loop is said to be slower than the apply statement. When writing a function that will be used only for printing output, it's probably ok to use the for-loop. For computationally intensive repeating operations, any kind of apply is better. An apply can be replaced with a snowfall parallel apply in the future.

The apply statement uses what's known as the "functional style" of programming. The statement says "take some list or vector or matrix and some function, and apply that function to each value." I believe R is not a good language to start learning functional programming. The difficulty lies in remembering the differences between apply, lapply, mapply, sapply, tapply, etc.