You are here

Binary Moderation

5 posts / 0 new
Last post
kpm's picture
kpm
Offline
Joined: 04/16/2010 - 13:32
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 &amp; 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 &amp; 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 - 15:14
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 - 13:32
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?

neale's picture
Offline
Joined: 07/31/2009 - 15:14
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.

Steve's picture
Offline
Joined: 07/30/2009 - 14:03
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.