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(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 GENETIC MODEL
Run Cholesky Decomposition ACE model
CholAceFit <- mxRun(CholAceModel,intervals=T)
CholAceSum <- summary(CholAceFit)
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
Try replacing it with
Create a constraint,
and be sure to include it in
See if these changes make any improvement. If not, please reply with the error message you see, as well as the output of
mxVersion()
.