Fri, 04/16/2010 - 13:49

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

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.

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?

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.

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.