You are here

Cholesky decomposition between a binary and a continuous variable?

2 posts / 0 new
Last post
RFrank's picture
Joined: 06/08/2012 - 14:22
Cholesky decomposition between a binary and a continuous variable?

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 ( and listed it below, but it's not running.

source( "")

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")

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(IV)), 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 Cholesky Decomposition ACE model

CholAceFit <- mxRun(CholAceModel,intervals=T)
CholAceSum <- summary(CholAceFit)

RobK's picture
Joined: 04/19/2011 - 21:00
What specifically, goes wrong

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, 

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().