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)
What specifically, goes wrong
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()
.Log in or register to post comments