Binary Moderation

4 replies [Last post]
kpm's picture
kpm
Offline
Joined: 04/16/2010

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

neale's picture
Offline
Joined: 07/31/2009
Hi Kristine I don't see

Hi Kristine

I don't see anything obviously wrong with your script. It is not possible for me to debug it by trying it out because the data files were not attached. From experience, parameters that do not change the -2lnL end up not changing from their starting values. You could establish this by fixing all the parameters and just changing the parameter in question to see if it does change the -2lnL. I would expect that it does not and that is why optimization has failed. Possibly, the data are the culprit. If, for example, all the the definition variables were zero, then whatever the values of ACE.aM, ACE.cM, ACE.eM and ACE.b then no change to the -2lnL would occur.

kpm's picture
kpm
Offline
Joined: 04/16/2010
Thanks so much! That does

Thanks so much! That does seem to be the problem. The definition variable is very skewed. Quick follow up question: Does the definition variable have to be relatively even (0's and 1's) within each sibling type? Or does it only have to be relatively even across the whole sample in order for optimization to succeed?

Steve's picture
Offline
Joined: 07/30/2009
It depends on your definition

It depends on your definition of "succeed". As Mike noted, if all of the values in your definition variable were zero (or all 1 for that matter), there is no variance that can be used to predict the effect of moderation.

Now suppose you have 500 zeros in your definition variable and you change one of those numbers from 0 to 1. You have variance to make a prediction, but the variance is very small. While from a strictly mathematical standpoint the optimization would succeed, you would find very little change from the estimated parameters.

As the number of 0s and 1s becomes more equal, you have greater power to detect a moderation effect because you have more variance in your predictor. This is not a property of the optimizer, it's a property of the numbers.

neale's picture
Offline
Joined: 07/31/2009
I would say that neither is

I would say that neither is needed. Also, please please note that non-zero status code 6 does not necessarily mean that optimization has failed. It might have, but it might be just fine.