You are here

Bivariate Cholesky - Continuous & Ordinal some errors

8 posts / 0 new
Last post
Fengxia's picture
Offline
Joined: 11/02/2017 - 03:33
Bivariate Cholesky - Continuous & Ordinal some errors
AttachmentSize
Binary Data BivCon&Ord.R16.48 KB

Hello,

I am trying run bivariate cholesky model analyses with my data in OpenMx2.5.0 recently. But I've got a couple of questions:

1)When I ran the code "CorFit <-mxRun(CorModel, intervals=F)", I got error message"Error in eval(result, envir = env) : object 'MZ.rMZ' not found".

2)I would also like to try models that include e.g. CE for variable 1 and ACE for variable 2,I used
ACE_noa1_Cholesky_Model <- omxSetParameters(model=AceFit,labels=c('a11','a21'),free=F,values=0,name="ACE_noa1_Cholesky").
But I got error message"Error: The following error occurred while evaluating the subexpression 'solve(sqrt(ACE_noa1_Cholesky.ACE_noa1_Cholesky.I * ACE_noa1_Cholesky.ACE_noa1_Cholesky.A))' during the evaluation of 'ACE_noa1_Cholesky.Ra' in model 'ACE_noa1_Cholesky' : Lapack routine dgesv: system is exactly singular: U[1,1] = 0".

Thanks for some help!

Cheers,
Fengxia

 rm(list=ls())
#ls()

require(OpenMx)
require(psych)
 
Bivdata <- read.csv ('bivtwindata.csv', header=T, sep=",", na.strings="NA")
 
names(Bivdata)
summary(Bivdata)
str(Bivdata)
describe(Bivdata)
 
BivdataOrd      <- data.frame(Bivdata) 
 
#Make the integer variables ordered factors (Note:because of the 'cuting' above, this is not necessary here)
BivdataOrd$sleepda1     <-mxFactor(BivdataOrd$sleepda1, levels=c(0:2) )
BivdataOrd$sleepda2    <-mxFactor(BivdataOrd$sleepda2, levels=c(0:2) )
 
describe(BivdataOrd)
 
vars        <-c('Hb','sleepda')
selVars <-c('Hb1', 'sleepda1', 'Hb2', 'sleepda2' )
useVars <-c('Hb1', 'sleepda1', 'Hb2', 'sleepda2', 'age1', 'age2')
 
# Select Data for Analysis
mzData <- subset(BivdataOrd, zygosity==1, useVars)
dzData <- subset(BivdataOrd, zygosity==2, useVars)
describe(mzData)
describe(dzData)
 
# To get the CTs for ADHD
table(mzData$sleepda1, mzData$sleepda2)
table(dzData$sleepda1, dzData$sleepda2)
 
nv  <- 2       # number of variables per twin
nvo  <- 1       # number of ordinal variables per twin
nvc  <- nv-nvo  # number of continuous variables per twin
poso <- nvc+1  # position where ordinal variables start
ntv  <- nv*2   # number of variables per pair
nth  <- 2       # number of max thresholds
ninc <- nth-1  # number of max increments

## -----------------------------------------------------------------------------------------
# Part 1:
# Constrained Polyserial correlation model 
# MZ and DZ Thresholds, but constrained across twins
# Age effect is different acros variables, but same across thresholds within ordinal variable(s)if c>2
# There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric
## ------------------------------------------------------------------------------------------------------------------------------

# CREATE LABELS in objects(to ease specification)

LabThMZ <-c('Tmz_11','imz_11')         # THs for ordinal var for a twin individual (mz)
LabThDZ  <-c('Tdz_11','idz_11')         # THs for ordinal var for a twin individual (dz)

LabCorMZ    <-c('r21','rMZ1','MZxtxt','MZxtxt','rMZ2','r21')
LabCorDZ    <-c('r21','rDZ1','DZxtxt','DZxtxt','rDZ2','r21')
 
(LabCovM   <-c( paste("Bcv",1:nvc,sep=""), rep(NA, nvo)))
(LabCovTH  <-rep( paste("Bov",1:nvo,sep=""), nth)) 
 
(LabM      <- c( paste("m",1:nv,sep="")))
(LabSD <- c( paste("sd",1:nv,sep="")))
(PatSD <- c( rep(T,nvc), rep(F, nvo)))
(PatM      <- c( rep(T,nvc), rep(F, nvo)))
(PatBetaM  <- c( rep(T,nvc), rep(F, nvo)))
 
# CREATE START VALUES in objects
#(StM      <-c( colMeans(mzData[,1:nvc],na.rm=TRUE), rep(0,nvo)))
StM      <-c( 100, 0)
(StSD  <-c( sd(mzData[,1:nvc],na.rm=TRUE), rep(1,nvo)) )
(StBetaM   <-c( rep(.2,nvc), rep(0,nvo)))
 
StCorMZ <-c(-.3, .7, -.3,-.3,   .7, -.3)
StCorDZ <-c(-.3, .4, -.15,-.15, .3, -.3)
StTH        <-c(-1.4, 1.7 )
 
# Define definition variables to hold the Covariates
obsAge1  <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age1"), name="Age1")
obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2")
 
# Matrix & Algebra for expected observed means, Thresholds, effect of Age on Th and correlations
Mean <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatM, values=StM, labels=LabM, name="M" )
 
betaAm  <-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaM, values=StBetaM, labels=LabCovM, name="BageM" )
betaAth <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTH, name="BageTH" )
 
Mu1 <-mxAlgebra( expression= M + t(BageM%*%Age1), name="MeanCtw1" )
Mu2 <-mxAlgebra( expression= M + t(BageM%*%Age2), name="MeanCtw2" )
Mu12    <-mxAlgebra( expression= cbind(MeanCtw1,MeanCtw2), name="ExpMean" )
 
inc     <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")
Tmz     <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThMZ, name="ThMZ")
ThresMZ <-mxAlgebra( expression= cbind(Low%*%ThMZ + BageTH%x%Age1, Low%*%ThMZ + BageTH%x%Age2), name="ExpThresMZ")
 
Tdz     <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThDZ, name="ThDZ")
ThresDZ <-mxAlgebra( expression= cbind(Low%*%ThDZ + BageTH%x%Age1, Low%*%ThDZ + BageTH%x%Age2), name="ExpThresDZ")
 
SD  <-mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=c(PatSD,PatSD), values=c(StSD,StSD), labels=c(LabSD,LabSD), name="sd") 
rMZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZ, lbound=-.99, ubound=.99, name="CorMZ") 
rDZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZ, lbound=-.99, ubound=.99, name="CorDZ") 
 
# Matrices for the Covariance model        
MZCov    <-mxAlgebra( expression=sd %&% CorMZ, name="ExpCovMZ" )
DZCov   <-mxAlgebra( expression=sd %&% CorDZ, name="ExpCovDZ" )
 
# Data objects 
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="ExpThresMZ", threshnames=c("sleepda1","sleepda2"))
objdz   <-mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThresDZ", threshnames=c("sleepda1","sleepda2"))
 
# Combine Groups
pars     <-list (obsAge1,obsAge2,Mean,betaAm,betaAth,Mu1,Mu2,Mu12,inc,SD)
modelMZ <-mxModel(pars, Tmz, ThresMZ, rMZ, MZCov, DataMZ, objmz, name="MZ" )
modelDZ <-mxModel(pars, Tdz, ThresDZ, rDZ, DZCov, DataDZ, objdz, name="DZ" )
minus2ll    <-mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj     <-mxAlgebraObjective("m2LL" )
Conf1       <-mxCI (c  ('MZ.rMZ[2,1]' )  )                     # IQ-ADHD phenotypic correlations
Conf2        <-mxCI (c  ('MZ.rMZ[3,1]', 'MZ.rMZ[4,2]','MZ.rMZ[3,2]')  )     # MZ twin cor var1, var 2 and xtr-xtwin
Conf3        <-mxCI (c  ('DZ.rDZ[3,1]', 'DZ.rDZ[4,2]','DZ.rDZ[3,2]')  )     # DZ twin cor var1, var 2 and xtr-xtwin
CorModel <-mxModel( "Cor", modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2, Conf3) 
 
# RUN Correlation MODEL 

CorFit  <-mxRun(CorModel, intervals=F)
(CorSumm   <-summary(CorFit))
 
round(CorFit@output$estimate,4)
CorFit$Cor$MZ.Rph
CorFit$Cor$MZ.Rbmz
CorFit$Cor$DZ.Rbdz
 
 
# ********************************************************************************************************

# -------------------------------------------------------------------------------------------------------------------
# (2) Specify Bivariate ACE Model using Cholesky Decomposition, 
# One set of THRESHOLDS (same across twins and zyg group)
# -------------------------------------------------------------------------------------------------------------------

LabTh       <-c('T_11','i_11')         # THs and inc for ordinal var for a twin individual 

# CREATE LABELS FOR Lower Triangular Matrices
(aLabs  <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep=""))
(cLabs     <- paste("c", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep=""))
(eLabs     <- paste("e", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep=""))
 
# CREATE START VALUES FOR Lower Triangular Matrices
(SDcon  <-(sqrt(vech(cov(mzData[,1:nvc],mzData[,1:nvc],use="complete")))))/3
Stpaths <- c(5,-0.2,3)
 
# Define definition variables to hold the Covariates
obsAge1  <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age1"), name="Age1")
obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2")
 
# Matrix & Algebra for expected observed means, Thresholds, effect of Age on Th and correlations
Mean     <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatM, values=StM, labels=LabM, name="M" )
 
betaAm  <-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaM, values=StBetaM, labels=LabCovM, name="BageM" )
betaAth <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTH, name="BageTH" )
 
Mu1     <-mxAlgebra( expression= M + t(BageM%*%Age1), name="MeanCtw1" )
Mu2     <-mxAlgebra( expression= M + t(BageM%*%Age2), name="MeanCtw2" )
Mu12        <-mxAlgebra( expression= cbind(MeanCtw1,MeanCtw2), name="ExpMean" )
 
inc     <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")
T       <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabTh, name="Th")
Thres       <-mxAlgebra( expression= cbind(Low%*%Th + BageTH%x%Age1, Low%*%Th + BageTH%x%Age2), name="ExpThres")
 
# MATRICES FOR THE ACE COV MODEL        
# Matrices to store a, c, and e Path Coefficients
pathA    <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=aLabs, name="a" ) 
pathC   <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=cLabs, name="c" )
pathE   <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=eLabs, 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 standardized variance components
covP <- mxAlgebra( expression=A+C+E, name="V" )
StA <- mxAlgebra( expression=A/V, name="h2")
StC <- mxAlgebra( expression=C/V, name="c2")
StE <- mxAlgebra( expression=E/V, name="e2")
 
# Algebra to compute Phenotypic, A, C & E correlations
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
rph <- mxAlgebra( expression= solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)), name="Rph")
rA  <- mxAlgebra( expression= solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)), name="Ra" )
rC  <- mxAlgebra( expression= solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)), name="Rc" )
rE  <- mxAlgebra( expression= solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)), name="Re" )
 
# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covMZ    <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C),
                                        cbind(A+C  , A+C+E))  , name="ExpCovMZ" )
covDZ   <- mxAlgebra( expression= rbind( cbind(A+C+E      , 0.5%x%A+C),
                                        cbind(0.5%x%A+C , A+C+E)),   name="ExpCovDZ" )
 
# Unity constraint on total variance of Ordinal variable(s)
mUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv" )
varL    <- mxConstraint( expression=diag2vec(V[(poso:nv),(poso:nv)])==Unv, name="VarL" )
 
# Data objects 
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="ExpThres", threshnames=c("sleepda1","sleepda2"))
objdz   <-mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThres", threshnames=c("sleepda1","sleepda2"))
 
# Combine Groups
pars1        <-list(obsAge1,obsAge2,Mean,betaAm,betaAth,Mu1,Mu2,Mu12,inc,T,Thres,mUnv,varL)
pars2       <-list(pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, rph, rA, rC, rE )
modelMZ <- mxModel( pars1, pars2, covMZ, DataMZ, objmz, name="MZ" )
modelDZ <- mxModel( pars1, pars2, covDZ, DataDZ, objdz, name="DZ" )
minus2ll    <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj     <- mxAlgebraObjective( "m2LL" )
Conf1       <- mxCI (c ('h2[1,1]', 'h2[2,2]', 'c2[1,1]', 'c2[2,2]', 'e2[1,1]', 'e2[2,2]') )
Conf2       <- mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Rc[2,1]', 'Re[2,1]') )
AceModel    <- mxModel( "ACE", pars2, modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2)
 
 
# ------------------------------------------------------------------------------
# RUN Constrained Bivariate ACE Model

AceFit   <- mxRun(AceModel, intervals=TRUE)
(AceSumm  <- summary(AceFit))
round(AceFit@output$estimate,4)
round(AceFit$Est@result,4)
 
AceFit$ACE.h2
AceFit$ACE.c2
AceFit$ACE.e2
AceFit$ACE.Rph
AceFit$ACE.Ra
AceFit$ACE.Rc
AceFit$ACE.Re
 
# Generate AceModel Output

# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
mxCompare(CorFit,AceFit)
 
# Fit ACE_noa1 model
# -----------------------------------------------------------------------

ACE_noa1_Cholesky_Model <- omxSetParameters(model=AceFit,labels=c('a11','a21'),free=F,values=0,name="ACE_noa1_Cholesky")
 
 
ACE_noa1_Cholesky_Fit <- mxRun(ACE_noa1_Cholesky_Model, intervals=TRUE)
ACE_noa1_Cholesky_Summ <- summary(ACE_noa1_Cholesky_Fit)
ACE_noa1_Cholesky_Summ
 
mxCompare(AceFit,ACE_noa1_Cholesky_Fit)
tableFitStatistics(AceFit,ACE_noa1_Cholesky_Fit)
AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Problems with algebras included when they aren't needed

Hi

There's an issue with including expressions like the following in the model itself. I recommend using mxEval() to calculate them after the model has been run. That way, if, e.g., A is set to zero the script still runs, and the only thing that fails is the mxEval() step (which fails because the A matrix is not positive definite when it has zeroes in it).

rph <- mxAlgebra( expression= solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)), name="Rph")
rA  <- mxAlgebra( expression= solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)), name="Ra" )
rC  <- mxAlgebra( expression= solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)), name="Rc" )
rE  <- mxAlgebra( expression= solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)), name="Re" )

The

solve(sqrt(I *A))

is the problem here. I note also that

solve(sqrt(I *A)) %&% A 

would be about half the computational effort (only do the inverse once and use the quadratic operator), but would still break when some diagonal elements of A are set to zero. Putting it in an mxEval() after the optimization would save MUCH more cpu effort because it isn't needed during optimization.

Fengxia's picture
Offline
Joined: 11/02/2017 - 03:33
I am a novice in OpenMx,

I am a novice in OpenMx, would you please show me the syntax how to use the mxEval()?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
mxVersion() ?

What version of OpenMx are you running (use the mxVersion() function to see)? Your syntax and the error message you report both make me strongly suspect that you're running version 1.x. I believe this issue is resolved in recent minor versions of v2.x , such that, as long as the offending MxAlgebra isn't strictly necessary, no error is thrown and the algebra is not computed (e.g., a previous post of mine).

Fengxia's picture
Offline
Joined: 11/02/2017 - 03:33
OpenMx 2.5.2

I used OpenMx 2.5.2.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
As for the first error

As for the first error you report,

1)When I ran the code "CorFit <-mxRun(CorModel, intervals=F)", I got error message"Error in eval(result, envir = env) : object 'MZ.rMZ' not found".

, the problem is that your syntax nowhere defines an object with the "Mx name" 'rMZ'. This line,

rMZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZ, lbound=-.99, ubound=.99, name="CorMZ")

creates an object with R symbol rMZ, but Mx name 'CorMZ'. I think all you need to do is change its Mx name to 'rMZ', and also change the Mx name of the DZ correlation matrix to 'rDZ'.

Fengxia's picture
Offline
Joined: 11/02/2017 - 03:33
I have made the change

Thank you so much! I have made the change and it worked.

And I still need some help about the second error.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
CIs for 'Ra' won't work.
And I still need some help about the second error.

With regard to the second error: your situation is very similar to that of JuanJMV in a different thread. Juan wanted to drop shared-environmental variance from only one of three traits. Consequently, the C covariance matrix had a diagonal element of zero, making sqrt(I C) uninvertible, and thereby making the shared-environmental correlation matrix for all three traits non-computable. In your case, you are similarly dropping additive-genetic variance from only one of two traits, making sqrt(I A) uninvertible. I think the real difference between your case and his is that you're requesting confidence intervals for the MxAlgebra named 'Ra', which is non-computable in your MxModel ACE_noa1_Cholesky_Model. Off the top of my head, the simplest solution is probably to define ACE_noa1_Cholesky_Model differently, say:

ACE_noa1_Cholesky_Model <- omxSetParameters(model=AceFit,labels=c('a11','a21'),free=F,values=0,name="ACE_noa1_Cholesky")
ACE_noa1_Cholesky_Model@intervals <- list()
ACE_noa1_Cholesky_Model <- mxModel(ACE_noa1_Cholesky_Model, Conf1, mxCI (c ('Rph[2,1]', 'Rc[2,1]', 'Re[2,1]') ))

. It wouldn't hurt to update to the latest version of OpenMx, either, but I don't anymore think that's the real issue.

I am a novice in OpenMx, would you please show me the syntax how to use the mxEval()?

See the documentation for mxEval(). There's some example syntax here. However, in your case, mxEval() is not going to be useful--you'd be requesting a standardized 1x1 matrix, which of course will contain a single 1.0.

Log in or register to post comments