# Binary Moderation

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

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.

Log in or register to post comments

In reply to Hi Kristine I don't see by neale

## Thanks so much! That does

Log in or register to post comments

In reply to Thanks so much! That does by kpm

## I would say that neither is

Log in or register to post comments

In reply to Thanks so much! That does by kpm

## It depends on your definition

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.

Log in or register to post comments