| Attachment | Size |
|---|---|
 BivCon&Ord.R | 16.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)
Printer-friendly version
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).
The
is the problem here. I note also that
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.
I am a novice in OpenMx, would you please show me the syntax how to use the mxEval()?
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).I used OpenMx 2.5.2.
As for the first error you report,
, the problem is that your syntax nowhere defines an object with the "Mx name" 'rMZ'. This line,
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'.Thank you so much! I have made the change and it worked.
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, 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_Modeldifferently, say:. 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.Hello,
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
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?
Thank you so much!
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?
I think I understand 2) now. I guess the answer is simply that a correlation (or a covariance) only represents linear association, and can't really say anything about nonlinear association. I'd have to think about how the twin model you're using could be extended to incorporate the "U"-shaped relationship you describe.