Login

News

Bivariate Cholesky - Continuous & Ordinal some errors

Posted on
No user picture. Fengxia Joined: 11/02/2017
Forums
Behavioral Genetics Models
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

#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)?>

Replied on Sun, 11/05/2017 - 11:21
Picture of user. AdminNeale Joined: 03/01/2013

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.

Log in or register to post comments
Replied on Sun, 11/05/2017 - 13:09
Picture of user. AdminRobK Joined: 01/24/2014

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).
Log in or register to post comments
Replied on Sun, 11/05/2017 - 13:22
Picture of user. AdminRobK Joined: 01/24/2014

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

Log in or register to post comments
Replied on Tue, 11/07/2017 - 10:43
Picture of user. AdminRobK Joined: 01/24/2014

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
Replied on Wed, 02/07/2018 - 09:18
No user picture. Fengxia Joined: 11/02/2017

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

Log in or register to post comments
Replied on Mon, 02/19/2018 - 15:12
Picture of user. AdminRobK Joined: 01/24/2014

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?

Log in or register to post comments
Replied on Sun, 03/04/2018 - 11:44
No user picture. Fengxia Joined: 11/02/2017

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?

Log in or register to post comments
Replied on Tue, 03/20/2018 - 16:35
Picture of user. AdminRobK Joined: 01/24/2014

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.
Log in or register to post comments
  Login

News