You are here

Bivariate Cholesky - Continuous & Ordinal some errors

12 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.

Fengxia's picture
Offline
Joined: 11/02/2017 - 03:33
more about bivariate SEM

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
some answers
1)Wether opposite sex twins can be used in model?

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.

3)The difference among Cholesky model, Common Pathways model and Independent Pathways model, and How to choose which model is the best?

See my post in another thread. I recommend selecting the model with the smallest AIC.

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.

Now THAT is an interesting question. Could you say more about what you're trying to accomplish with this analysis?

Fengxia's picture
Offline
Joined: 11/02/2017 - 03:33
More about 2)& A NEW question

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?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
I think I understand 2) now.

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.