Cholesky decomposition between a binary and a continuous variable?

Posted on
No user picture. RFrank Joined: 06/08/2012
Is anyone aware of a script to do a Cholesky decomposition between a binary and a continuous variable?

I attempted to adapt a script I found here (http://www.genepi.qimr.edu.au/staff/sarahMe/workshop14/Day3.zip) and listed it below, but it's not running.

require(OpenMx)
require(psych)
source( "http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")

nl <- read.delim("Cholesky.dat", header=T, sep="\t", na.strings=".")

nl$family1 <- nl$fam
nl$family2 <- nl$famb
nl$happy1 <- nl$max
nl$happy2 <- nl$maxb

# Select Variables for Analysis
Vars <- c('family','happy')
nv <- 2 # number of variables
nind <- 2
ntv <- nv*2 # number of total variables
selVars <- c("family1", "happy1", "family2", "happy2")
#paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Data for Analysis

mzData <- subset(nl, ZYGOSITY==1, selVars)
dzData <- subset(nl, ZYGOSITY==0, selVars)

# Define threshold and fix factor ordering

mzData$family1 <- mxFactor(mzData$family1,levels=c(0,1))
mzData$family2 <- mxFactor(mzData$family2,levels=c(0,1))
dzData$family1 <- mxFactor(dzData$family1,levels=c(0,1))
dzData$family2 <- mxFactor(dzData$family2,levels=c(0,1))
threshold <-mxMatrix( type="Full", nrow=1, ncol=2, free=T, values=1, labels="thresh", name="Threshold" )
inclusions <- list (threshold)

# Matrices declared to store a, c, and e Path Coefficients

pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=.6, labels=c("a11","a21","a22"), name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=.6, labels=c("c11","c21","c22"), name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=.6, labels=c("e11","e21","e22"), name="e" )

# Matrices generated to hold A, C, and E computed Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )

# Algebra to compute total variances and standard deviations (diagonal only)
covP <- mxAlgebra( expression=A+C+E, name="V" )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")
StPathA <- mxAlgebra(iSD %*% a, name="sta")
StPathC <- mxAlgebra(iSD %*% c, name="stc")
StPathE <- mxAlgebra(iSD %*% e, name="ste")

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=.6, labels=c("family","happy"), name="expMean" )
covMZ <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZ" )
covDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZ" )

# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="Threshold", threshnames=c('family1','family2') )
objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="Threshold", threshnames=c('family1','family2') )

# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, meanG,StPathA,StPathC,StPathE)
modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj )

# ------------------------------------------------------------------------------
# RUN GENETIC MODEL
# Run Cholesky Decomposition ACE model
CholAceFit <- mxRun(CholAceModel,intervals=T)
CholAceSum <- summary(CholAceFit)

Replied on Fri, 10/17/2014 - 11:45
Picture of user. RobK Joined: 04/19/2011

What specifically, goes wrong when you try to run it?

One thing I notice immediately is that the mean, variance, and threshold of the binary variable are all free parameters. The model isn't identified. You'll need to fix its variance (or its nonshared-environmental variance) to some value, say, 1. You'll also need to either fix the the mean or the threshold to some value, say, 0. To these ends, let me make a few suggestions.

You currently have

meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,values=.6, labels=c("family","happy"), name="expMean")

Try replacing it with

meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv,
free=c(F,T,F,T),
values=c(0,.6,0,.6),
labels=c("family1","happy1","family2","happy2"),
name="expMean")

Create a constraint,

vconst <- mxConstraint(V[1,1]==1, name="Vconst")

and be sure to include it in

pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, meanG,StPathA,StPathC,StPathE,vconst)

See if these changes make any improvement. If not, please reply with the error message you see, as well as the output of mxVersion().