Binary Moderation

Hello, I'm trying to adapt the moderation script from the twin workshop in Boulder. At the Boulder workshop the binary moderator was divorce (0/1, and identical for twin 1 and twin 2). The moderator I'm using also is scored 0/1 but may differ for twin 1 and twin 2. Everything else is the same, and I changed very little, but instead of estimating moderation effects (aM, cM, eM), those parameters always come out as whatever I put in for the starting values. Parameter specifications verify that I am estimating them, but clearly something is breaking down before then in the script. I have a feeling it has to do with the definition variables, but I'm at a loss for how to find or fix it. Any suggestions?

Here's my script:

Thanks,
Kristine

# -----------------------------------------------------------------------
# Program: VCxAge.R
# Author: Hermine Maes
# Date: 12 01 2009
#
#
# Revision History
# Hermine Maes -- 02 01 2010 updated & reformatted
# Tim York / Danielle Dick -- 02 24 2010 modified for Mx workshop
# Marleen de Moor / Tim York -- 03 03 2010
# Modified for TCHAD puberty/depression - KPM
# -----------------------------------------------------------------------

# Univariate Heterogeneity Twin Analysis model to estimate causes of variation (ACE)
# with Heterogeneity in ACE variance decompositions for prepubertal and pubertal adolescents
# and one mean estimated
# Matrix style model input - Raw data input
#
# Datafile: TCHAD.dat
# Phenotype: depression (dep)
# Heterogeneity variable: puberty (0=pre-pubertal, 1=pubertal)
# Zygosity variable: zyg (1=MZ, 2=DZ)
#
# -----------------------------------------------------------------------
require(OpenMx)
require(psych)
source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")

#=======================================================================#
# PREPARE DATA #
#=======================================================================#

# General Family Functioning Data
data <- read.table("TCHAD.txt", header=F, na.strings=-99)
names(data) <- c("type", "sex", "dep1", "dep2", "pub1", "pub2", "pub1m", "pub2m")

mzData <- as.data.frame(subset(data, type==1&sex==2, c(dep1,dep2,pub1,pub2)))
dzData <- as.data.frame(subset(data, type==2&sex==2, c(dep1,dep2,pub1,pub2)))

head(mzData)

#CREATE DEFINITION VARIABLE FOR EACH TWIN
#mzData$pub1 <- mzData$pub1
#mzData$pub2 <- mzData$pub2
#dzData$pub1 <- dzData$pub1
#dzData$pub2 <- dzData$pub2

selVars <- c('dep1','dep2')
nv <- 1
ntv <- nv*2

TWOgroupB <- mxModel("groups2B",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=2.5, label="a11", name="a" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=2.5, label="c11", name="c" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=2.5, label="e11", name="e" ),
# Matrices a, c, and e to store moderated a, c, and e path coefficients
# These are the BetaMs that are on the path diagrams -KPM
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="aM11", name="aM" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cM11", name="cM" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="eM11", name="eM" ),
# Matrix & Algebra for expected means vector (non-moderated)
mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 38, label="mean", name="mu" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.3, label=c("l11"), name="b" ),
# Matrices A, C, and E compute non-moderated variance components
mxAlgebra( expression=a %*% t(a), name="A" ),
mxAlgebra( expression=c %*% t(c), name="C" ),
mxAlgebra( expression=e %*% t(e), name="E" ),
# Algebra to compute total variances and inverse of standard deviations (diagonal only)
mxAlgebra( expression=A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(I*V)), name="isd")
),

mxModel("MZ",
# Matrix for moderating/interacting variable
# This is calling in the definition variables that are the moderator
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("mzData.pub1"), name="D1"), #twin1
# data.divorce1 means look in our data for the divorce1 variable.
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("mzData.pub2"), name="D2"), #twin2
# Matrices A, C, and E compute variance components
mxAlgebra( expression=(ACE.a+ D1%*%ACE.aM) %*% t(ACE.a+ D1%*%ACE.aM), name="A1" ),
mxAlgebra( expression=(ACE.c+ D1%*%ACE.cM) %*% t(ACE.c+ D1%*%ACE.cM), name="C1" ),
mxAlgebra( expression=(ACE.e+ D1%*%ACE.eM) %*% t(ACE.e+ D1%*%ACE.eM), name="E1" ),
mxAlgebra( expression=(ACE.a+ D2%*%ACE.aM) %*% t(ACE.a+ D2%*%ACE.aM), name="A2" ),
mxAlgebra( expression=(ACE.c+ D2%*%ACE.cM) %*% t(ACE.c+ D2%*%ACE.cM), name="C2" ),
mxAlgebra( expression=(ACE.e+ D2%*%ACE.eM) %*% t(ACE.e+ D2%*%ACE.eM), name="E2" ),
# Algebra for expected variance/covariance matrix and expected mean vector in MZ
mxAlgebra( expression= rbind ( cbind(A1+C1+E1 , (ACE.a+ D2%*%ACE.aM) %*% t(ACE.a+ D1%*%ACE.aM) + (ACE.c+ D2%*%ACE.cM) %*% t(ACE.c+ D1%*%ACE.cM)),
# Moderation of twin 2 times moderation of twin 1 added to c mod of t1 * c mod of twin 2 = moderation term
cbind((ACE.a+ D1%*%ACE.aM) %*% t(ACE.a+ D2%*%ACE.aM) + (ACE.c+ D1%*%ACE.cM) %*% t(ACE.c+ D2%*%ACE.cM), A2+C2+E2)), name="expCovMZ" ),
mxAlgebra( expression= ACE.b %*% D1, name="D1R"),
mxAlgebra( expression= ACE.b %*% D2, name="D2R"),
mxAlgebra( expression= cbind((ACE.mu + D1R),(ACE.mu + D2R)), name="expMean"),
# Data & Objective
mxData( observed=mzData[,c(selVars,"pub1","pub2")], type="raw" ),
mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars )
),

mxModel("DZ",
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("dzData.pub1"), name="D1"), #twin1
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("dzData.pub2"), name="D2"), #twin2
# Matrices A, C, and E compute variance components
mxAlgebra( expression=(ACE.a+ D1%*%ACE.aM) %*% t(ACE.a+ D1%*%ACE.aM), name="A1" ),
mxAlgebra( expression=(ACE.c+ D1%*%ACE.cM) %*% t(ACE.c+ D1%*%ACE.cM), name="C1" ),
mxAlgebra( expression=(ACE.e+ D1%*%ACE.eM) %*% t(ACE.e+ D1%*%ACE.eM), name="E1" ),
mxAlgebra( expression=(ACE.a+ D2%*%ACE.aM) %*% t(ACE.a+ D2%*%ACE.aM), name="A2" ),
mxAlgebra( expression=(ACE.c+ D2%*%ACE.cM) %*% t(ACE.c+ D2%*%ACE.cM), name="C2" ),
mxAlgebra( expression=(ACE.e+ D2%*%ACE.eM) %*% t(ACE.e+ D2%*%ACE.eM), name="E2" ),

# Algebra for expected variance/covariance matrix in DZ
mxAlgebra( expression= rbind ( cbind(A1+C1+E1 , 0.5%x%((ACE.a+ D2%*%ACE.aM) %*% t(ACE.a+ D1%*%ACE.aM)) + (ACE.c+ D2%*%ACE.cM) %*% t(ACE.c+ D1%*%ACE.cM)),
cbind(0.5%x%((ACE.a+ D1%*%ACE.aM) %*% t(ACE.a+ D2%*%ACE.aM)) + (ACE.c+ D1%*%ACE.cM) %*% t(ACE.c+ D2%*%ACE.cM), A2+C2+E2)), name="expCovDZ" ),
mxAlgebra( expression= ACE.b %*% D1, name="D1R"),
mxAlgebra( expression= ACE.b %*% D2, name="D2R"),
mxAlgebra( expression= cbind((ACE.mu + D1R),(ACE.mu + D2R)), name="expMean"),
# Data & Objective
mxData( observed=dzData[,c(selVars,"pub1","pub2")], type="raw" ),
mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="-2sumll" ),
mxAlgebraObjective("-2sumll")
)

TWOgroupBFIT <- mxRun(TWOgroupB)
TWOgroupBSUM <- summary(TWOgroupBFIT)

#-----------------------------------------------------------------------#
# Generate TWOgroupB Moderated Output #
#-----------------------------------------------------------------------#
parameterSpecifications(TWOgroupBFIT)
expectedMeansCovariances(TWOgroupBFIT)
tableFitStatistics(TWOgroupBFIT)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.a,ACE.c,ACE.e,ACE.aM,ACE.cM,ACE.eM), TWOgroupBFIT),4)
varComponentsACE <- round(mxEval(cbind(ACE.A/ACE.V,ACE.C/ACE.V,ACE.E/ACE.V), TWOgroupBFIT),4)
rownames(pathEstimatesACE) <- 'pathEstimates'
colnames(pathEstimatesACE) <- c('a','c','e','aM','cM','eM')
rownames(varComponentsACE) <- 'varComponents'
colnames(varComponentsACE) <- c('a^2','c^2','e^2')
pathEstimatesACE
varComponentsACE