Extracting means and covariances from a model

10 replies [Last post]
mspiegel's picture
Offline
Joined: 07/31/2009

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)
}

mspiegel's picture
Offline
Joined: 07/31/2009
Here is a cleaned up version

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

  • Avoid eval() and parse() if possible. Don't make a vector of variable names (as character strings) and use eval and parse to get the variables. Instead make a list of the variables themselves.
  • Avoid multiple calls to the summary() function. For some important reason I can't remember, summary needs to flatten a model. This is OK as a one time cost, but is expensive if done repeatedly. Store the results of summary() in a variable.
  • The function printFitStatistics1() takes two arguments. If the second argument is missing, then print only the statistics for the one model. If the second argument is available, then print the comparative statistics. This eliminates the need for printFitStatistics2().
  • I recommend renaming the functions printFitStatistics1() and printFitStatistics(). The names are misleading; the functions don't print to the terminal.
AttachmentSize
BivariateTwin_MatrixRaw.R 13.38 KB
Hermine's picture
Offline
Joined: 07/31/2009
I take your points. I really

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.

AttachmentSize
BivariateTwin_MatrixRaw2.R 12.78 KB
mspiegel's picture
Offline
Joined: 07/31/2009
I wrote showFitStatistics to

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.

AttachmentSize
BivariateTwin_MatrixRaw2.R 12.94 KB
ronjacobson's picture
Offline
Joined: 02/01/2010
Hey, you're doing great with

Hey, you're doing great with the matrix..
that is really cool! great work ^_^

Hermine's picture
Offline
Joined: 07/31/2009
Hi, The attached job runs a

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.

AttachmentSize
BivariateTwin_MatrixRaw.R 14.57 KB
neale's picture
Offline
Joined: 07/31/2009
This is very cool. However,

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"

AttachmentSize
BivariateTwinAnalysis_MatrixRaw.R 4.35 KB
mspiegel's picture
Offline
Joined: 07/31/2009
Ah ok. Try this change. The

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.

showMeansCovariances <- function(model) {
   resultsList <- .collectMeansCovariances(model, 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, topModel) {
   listReturn <- list()
   if (!is.null(model$objective)) {
      objective <- model$objective
      slots <- slotNames(objective)
 
      # extract the covariance
      if ('covariance' %in% slots) {
         covName <- objective@covariance
         if (length(grep('.', covName, fixed=TRUE)) == 1) {
            covariance <- eval(substitute(mxEval(x, topModel),
               list(x = as.symbol(covName))))
         } else {
            covariance <- eval(substitute(mxEval(x, model),
               list(x = as.symbol(covName))))
         }
         storeName <- paste('model:', model@name,
            ', covariance:', covName, sep='')
         listReturn[[storeName]] <- covariance
      }
 
      # extract the means
      if ('means' %in% slots) {
         meansName <- objective@means
         if (length(grep('.', meansName, fixed=TRUE)) == 1) {
            means <- eval(substitute(mxEval(x, topModel),
               list(x = as.symbol(meansName))))
         } else {
            means <- eval(substitute(mxEval(x, model),
               list(x = as.symbol(meansName))))
         }
         storeName <- paste('model:', model@name,
            ', means:', meansName, sep='')
         listReturn[[storeName]] <- means
      }      
   }
 
   # Recursively collect means and covariances of submodels
   names(model@submodels) <- NULL
   submodels <- lapply(model@submodels, .collectMeansCovariances, topModel)
   listReturn <- append(listReturn, unlist(submodels, FALSE))
   return(listReturn)
}

neale's picture
Offline
Joined: 07/31/2009
Works like a charm! Thank

Works like a charm! Thank you!

Hermine's picture
Offline
Joined: 07/31/2009
I tried to extend the

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.

mspiegel's picture
Offline
Joined: 07/31/2009
I've attached a script that

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.

AttachmentSize
specificationMatrix.R 1.27 KB