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) }
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")) {
}
}
# 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.
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.
Here is a cleaned up version of the script. A couple of general comments:
<
ul>
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 ^_^