genEpi Helper Functions

I'm trying to pull out estimates in a matrix using genEpi helper functions in a trivariate model with a sibling but I'm running into some trouble with the genEpi_FormatOutputMatrices function.

Error:
Error in print(genEpi_FormatOutputMatrix(matrix = genEpi_EvalQuote(expstring = matricesList[[k]], :
error in evaluating the argument 'x' in selecting a method for function 'print': Error: The following error occurred while evaluating the expression 'iSD %*% a' in model 'CholACE' : requires numeric/complex matrix/vector arguments

Full Script is below

require(OpenMx)
require(psych)
source("myFunctions.R")
source("http://openmxhelpers.googlecode.com/svn/trunk/GenEpiHelperFunctions.R")

# --------------------------------------------------------------------
# PREPARE DATA

# Load Data
data <- read.csv('Standardized Disgust Data.csv', na.strings='-99')

head(data)
describe(data)

data$twin1age <- data$age_NEW.1
data$twin2age <- data$age_NEW.2
data$sibage <- data$sib1age
data$moral1 <- data$ZDisgustMoral.1
data$moral2 <- data$ZDisgustMoral.2
data$sibdm <- data$Zsib1dm
data$path1 <- data$ZDisgustPathogen.1
data$path2 <- data$ZDisgustPathogen.2
data$sibdp <- data$Zsib1dp
data$sexual1 <- data$ZDisgustSexual.1
data$sexual2 <- data$ZDisgustSexual.2
data$sibds <- data$Zsib1ds
data$zyg <- data$Corrected_zygosity
head(data)

# Select Variables for Analysis
Vars <- c('moral','path','sexual')
nv <- 3 # number of variables
nsib <- 3
ntv <- nv*nsib # number of total variables

selVars <- c('moral1','path1','sexual1',
'moral2','path2','sexual2',
'sibdm','sibdp','sibds')

defVars <-c('twin1age','twin2age','sibage')
useVars <- c(selVars,defVars)

#Subset data for testing
mzdata <- subset(data, zyg==3, useVars)
dzdata <- subset(data, zyg==4, useVars)
describe(mzdata, skew=F)
describe(dzdata, skew=F)
dim(mzdata)
dim(dzdata)
cov(mzdata,use="complete")
cov(dzdata,use="complete")
cor(mzdata,use="complete")
cor(dzdata,use="complete")

# ------------------------------------------------------------------------------
# Cholesky Decomposition ACE Model
# ------------------------------------------------------------------------------
# Set Starting Values
svMe <- c(.01,.01, .01) # start value for means
svPa <- valDiag(nv,.6) # start values for parameters on diagonal
lbPa <- valLUDiag(nv,.0001,-10,NA) # lower bounds for parameters on diagonal

# Matrices declared to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("a",nv), lbound=lbPa, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("c",nv), lbound=lbPa, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("e",nv), lbound=lbPa, name="e" )

# Matrices generated to hold A, C, and E computed Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )

# Algebra to compute total variances and standard deviations (diagonal only)
covP <- mxAlgebra( expression=A+C+E, name="V" )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
inSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
# Algebra for expected Mean Matrices in MZ & DZ twins
#MeanG (grand mean) also known as the intercept
#only supplied one label - will be provided to both matrix elements and will supply same value
#Want intercept to be added to beta * age - regression outcome

defAge <- mxMatrix(type="Full", nrow=1,ncol=ntv,free=FALSE,
labels=c('data.twin1age','data.twin2age','data.sibage'),name='Age')

pathB <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= .01,
label="b11", name="b")

mean <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=labFull("me",1,nv), name="mean" )

expMean <- mxAlgebra( mean + (b%x%Age), name="expMean" )

#Covariance Matricies
covMZ <- mxAlgebra( expression= rbind( cbind(V , A+C, .5%x%A+C),
cbind(A+C, V, .5%x%A+C),
cbind(.5%x%A+C, .5%x%A+C, V)),
name="expCovMZ" )

covDZ <- mxAlgebra( expression= rbind( cbind(V , .5%x%A+C, .5%x%A+C),
cbind(.5%x%A+C, V, .5%x%A+C),
cbind(.5%x%A+C, .5%x%A+C, V)),
name="expCovDZ" )

# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzdata, type="raw" )
dataDZ <- mxData( observed=dzdata, type="raw" )

# Objective objects for Multiple Groups
objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars )
objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )

# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, defAge, pathB, mean, expMean)
modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
CholAceModel <- mxModel( "CholACE", modelMZ, modelDZ, minus2ll, obj )

# ------------------------------------------------------------------------------
# RUN GENETIC MODEL

# Run Cholesky Decomposition ACE model
CholAceFit <- mxRun(CholAceModel)
CholAceSumm <- summary(CholAceFit)

mxCompare(twinSatFit,CholAceFit)
mxCompare(eqMeVaZygFit,CholAceFit)

#Generate list of parameter estimates and derived quantities using formatOutputMatrices
# ACE standardized Path Coefficients (pre-multiplied by inverse of standard deviations)
CholACEpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e")
CholACEpathLabels <- c("stPathA","stPathC","stPathE")
#genEpi_FormatOutputMatrix(matrix,label,vars,digits)
genEpi_FormatOutputMatrices(CholAceFit,CholACEpathMatrices,CholACEpathLabels,Vars,4)