Bivariate Cholesky - Continuous & Ordinal some errors

Attachment | Size |
---|---|
BivCon&Ord.R | 16.48 KB |
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
#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)?>
Problems with algebras included when they aren't needed
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.
Log in or register to post comments
In reply to Problems with algebras included when they aren't needed by AdminNeale
I am a novice in OpenMx,
Log in or register to post comments
mxVersion() ?
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).Log in or register to post comments
In reply to mxVersion() ? by AdminRobK
OpenMx 2.5.2
Log in or register to post comments
As for the first error
, 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'.Log in or register to post comments
In reply to As for the first error by AdminRobK
I have made the change
And I still need some help about the second error.
Log in or register to post comments
In reply to I have made the change by Fengxia
CIs for 'Ra' won't work.
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, makingsqrt(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 MxModelACE_noa1_Cholesky_Model
. Off the top of my head, the simplest solution is probably to defineACE_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.
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
In reply to CIs for 'Ra' won't work. by AdminRobK
more about bivariate SEM
Recently I have some questions about bivariate model in OpenMx:
1)Wether opposite sex twins can be used in model?
2)The hypothesis of my study is that there is a ‘‘U’’ shaped relationship between sleep and HbA1c,that is short sleep duration and long sleep duration would all lead higher HbA1c. When I only have one rP and rG using bivariate model, how can I explain the relationship of sleep and HbA1c? Or there are some other methods? In my study, the sleep duration is divided into three groups.
3)The difference among Cholesky model, Common Pathways model and Independent Pathways model, and How to choose which model is the best?
Thanks for some help!
Cheers,
Fengxia
Log in or register to post comments
In reply to more about bivariate SEM by Fengxia
some answers
Yes. I see that you aren't adjusting for the main effect of sex in your script, nor are you separating groups by sex (only by zygosity). If you are sure that sex is irrelevant to the phenotypes (which I rather doubt), then you can include opposite-sex DZ twins without modification. Otherwise, you'll need to at least adjust for the main effect of sex, similarly to what you're currently doing with age.
See my post in another thread. I recommend selecting the model with the smallest AIC.
Now THAT is an interesting question. Could you say more about what you're trying to accomplish with this analysis?
Log in or register to post comments
In reply to some answers by AdminRobK
More about 2)& A NEW question
More about 2):
I have found that compared to sleeping 6~8h, both short sleep duration (<6h) and long sleep duration (>8h) were associated with a higher prevalence and of diabetes and higher HbA1c using logistic regression model and mixed linear model, that's the " ‘‘U’’ shaped relationship". Now I want to learn the traits (sleep duration and diabetes, sleep duration and HbA1c) may be correlated due to shared genetic factors or shared environmental factors. So I did this twin study using bivariate model. BUT the model only give one rP and rG with one direction, disagree with the "U" shaped relationship". I have no idea how to explain the result. OR could I dive it in to TWO models, that's one divariate model for diabetes and sleep duration (<6h and 6~8h) when another for diabetes and sleep duration (6~8h and >8h)? IF this is a feasible method, how about the syntax?
A NEW question:
In my study, the beast fitting univariate model of diabetes was ACE, while the sleep duration was ADE, than the Cholesky model using ACE-ADE at first right?
Log in or register to post comments
In reply to More about 2)& A NEW question by Fengxia
I think I understand 2) now.
Log in or register to post comments