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)