You are here

genEpi Helper Functions

4 posts / 0 new
Last post
James Sherlock's picture
Joined: 06/03/2014 - 01:45
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 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


# --------------------------------------------------------------------

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


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

# 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',

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)

# ------------------------------------------------------------------------------
# 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,

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 Cholesky Decomposition ACE model
CholAceFit <- mxRun(CholAceModel)
CholAceSumm <- summary(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")

RobK's picture
Joined: 04/19/2011 - 21:00
The objects you are asking to

The objects you are asking to be multiplied together are not in the fitted supermodel, CholAceFit, but in the submodels. There are a few small changes you could make to get the result you want.

One possibility is to put certain objects from pars back into the supermodel. I think they would be pathA, pathC, pathE, covA, covC, covE, covP, matI, and invSD. Then, your call to genEpi_FormatOutputMatrices() should work as-is.

Another possibility is to replace
CholACEpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e")
with (say)
CholACEpathMatrices <- c("MZ.iSD %*% MZ.a","MZ.iSD %*% MZ.c","MZ.iSD %*% MZ.e")

A third possibility, and probably the easiest, is to change
to (say)

James Sherlock's picture
Joined: 06/03/2014 - 01:45
problem solved

Thanks so much for your help Rob!

RobK's picture
Joined: 04/19/2011 - 21:00
You're welcome

No problem! It's partly my fault you got this error in the first place: I told you to take pars out of the supermodel entirely in that other thread, which was a bit of overkill, since the definition-variables problem could have been solved in any way that kept defAge out of the supermodel.

Log in or register to post comments