You are here

Intra Class Corelations

13 posts / 0 new
Last post
JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Intra Class Corelations

Hi I am doing some univariate analyses but I have one doubt when about how getting the intraclass correlations for MZ and DZ. I have tried with cov2cor but I can not make that work. Here is the script:

# Select Variables for Analysis
vars      <- "AG_Ln"                     # list of variables names
nv        <- 1                         # number of variables
ntv       <- nv*2                      # number of total variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
 
# Select Covariates for Analysis
#twinData[,'age'] <- twinData[,'age']/100
#twinData  <- twinData[-which(is.na(twinData$age)),]
covVars   <- c("age","Sex1","Sex2")
 
# Select Data for Analysis
mzData    <- subset(twinData, Zyg==1, c(selVars, covVars))
dzData    <- subset(twinData, Zyg==2, c(selVars, covVars))
 
# Set Starting Values
svMe      <- 4                       # start value for means
svPa      <- .4                        # start value for path coefficient
svPe      <- .8                        # start value for path coefficient for e
lbPa      <- .01                     # start value for lower bounds
 
# ------------------------------------------------------------------------------
# PREPARE MODEL
 
# ACE Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" )
pathB1     <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="b11", name="b1" )
defSex    <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.Sex1", "data.Sex2"), name="Sex" )
pathB2    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="b12", name="b2" )
 
# Create Algebra for expected Mean Matrices
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels="xbmi", name="meanG" )
expMean   <- mxAlgebra( expression= meanG + cbind(b1%*%Age,b1%*%Age)+b2%*%Sex, name="expMean" )
 
# Create Matrices for Path Coefficients
pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="a11", lbound=lbPa, name="a" ) 
pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="c11", lbound=lbPa, name="c" )
pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="e11", lbound=lbPa, name="e" )
 
# Create Algebra for 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" )
 
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP      <- mxAlgebra( expression= A+C+E, name="V" )
covMZ     <- mxAlgebra( expression= A+C, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" )
expCovMZ  <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ  <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )
 
# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )
 
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()
 
# Create Model Objects for Multiple Groups
pars      <- list(pathB1,pathB2, meanG, pathA, pathC, pathE, covA, covC, covE, covP)
defs      <- list(defAge, defSex)
modelMZ   <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
 
# Create Algebra for Variance Components
rowVC     <- rep('VC',nv)
colVC     <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estVC     <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC))
 
# Create Confidence Interval Objects
ciACE     <- mxCI( "VC" )
 
# Build Model with Confidence Intervals
modelACE  <- mxModel( "oneACEca", pars, modelMZ, modelDZ, multi, estVC, ciACE )
 
# ------------------------------------------------------------------------------
# RUN MODEL
 
# Run ACE Model
fitACE    <- mxTryHard( modelACE, intervals=T, extraTries = 10, finetuneGradient = T, exhaustive = T )
sumACE    <- summary( fitACE )
summary(fitACE)
# Compare with Saturated Model
mxCompare( fit, fitACE )
#lrtSAT(fitACE,-2ll,df)
 
# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEsts(fitACE)
 
# ------------------------------------------------------------------------------
# RUN SUBMODELS
 
# Test Significance of Covariate
modelCov  <- mxModel( fitACE, name="testCov" )
modelCov  <- omxSetParameters( modelCov, label=c("b11","b12"), free=FALSE, values=0 )
fitCov    <- mxRun( modelCov )
mxCompare( fitACE, fitCov )
 
# Run AE model
modelAE   <- mxModel( fitACE, name="oneAEca" )
modelAE   <- omxSetParameters( modelAE, labels="c11", free=FALSE, values=0 )
fitAE     <- mxRun( modelAE, intervals=T )
mxCompare( fitACE, fitAE )
fitGofs(fitAE)
fitEsts(fitAE)
summary(fitAE)
# Run CE model
modelCE   <- mxModel( fitACE, name="oneCEca" )
modelCE   <- omxSetParameters( modelCE, labels="a11", free=FALSE, values=0 )
fitCE     <- mxRun( modelCE, intervals=T )
mxCompare( fitACE, fitCE )
fitGofs(fitCE)
fitEsts(fitCE)
 
# Run E model
modelE    <- mxModel( fitAE, name="oneEca" )
modelE    <- omxSetParameters( modelE, labels="a11", free=FALSE, values=0 )
fitE      <- mxRun( modelE, intervals=T )
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)
 
# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitCov, fitAE, fitCE, fitE) )
round(rbind(fitACE$VC$result,fitAE$VC$result,fitCE$VC$result,fitE$VC$result),4)

Thank you so much.

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

So, what's your syntax involving cov2cor(), and what results does it give you?

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Hi,

Hi,

I have tried with this:

mxEval(cov2cor(V), fitACE, T) the result is: [,1] [1,] 1

and with: mxEval(cov2cor(expCovDZ), fitACE, T) the result is:
Error: The following error occurred while evaluating the expression 'cov2cor(expCovDZ)' in model 'oneACEca' : 'V' is not a square numeric matrix

Thank you so much!

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Ahhh

The first one is not an error. According to the script you posted, V is a 1x1 matrix. If you run cov2cor() on it, it should produce a 1x1 matrix with a 1 in it.

The second one is a namespace issue. mxEval is confused about the R object called expCovDZ which is an MxAlgebra and the OpenMx object given the same name which can be interpreted as a numeric matrix. Name the R object something else and you should get a 2x2 correlation matrix.

# change this line 
expCovDZ  <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )
 
# to something like
rexpCovDZ  <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )
 
# then put 'rexpCovDZ' into the model.
JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much for your

Thank you so much for your answer. I did the changes but now how can I get the value? sorry but I am very new

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Are you having a new problem?

Are you having a new problem? A new error message? Please describe what you are trying in as much detail as possible and what problem you are encountering, preferably with example code that demonstrates the issue.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
my guess
I did the changes but now how can I get the value?

I'm going to take a wild guess that all three of the following lines of code will do what you want for the DZ twins (and with the obvious changes, what you want for the MZ twins):

mxEval(cov2cor(expCovDZ), fitACE$DZ, T)
cov2cor(mxGetExpected(fitACE$DZ, "covariance"))
mxEval(cov2cor(DZ.expCovDZ), fitACE, T)

Am I correct?

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Yeah!!! you are right. Thank

Yeah!!! you are right. Thank you so much I had tried with that function but not properly.

Thanks for your help I really appreciate it.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Just one thing more I did the

Just one thing more I did the same in the trivariate model (you helped me in other thread) but I don't know which value is for each variable. Here is the result:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00000000 0.66182596 0.37644256 0.62216985 0.54028297 0.42317649
[2,] 0.66182596 1.00000000 0.32231799 0.54028297 0.70316889 0.32361090
[3,] 0.37644256 0.32231799 1.00000000 0.42317649 0.32361090 0.78676051
[4,] 0.62216985 0.54028297 0.42317649 1.00000000 0.66182596 0.37644256
[5,] 0.54028297 0.70316889 0.32361090 0.66182596 1.00000000 0.32231799
[6,] 0.42317649 0.32361090 0.78676051 0.37644256 0.32231799 1.00000000

So I should have 3 intra-class correlations but I don't know what value is for each variable.

Also, I would like to know if the intra-class correlations should be the same in the univariate analyses and trivariate analysis?

Thanks in advance

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
between-twin, within-trait?

You probably want the between-twin, within-trait correlations, which would be on the diagonal of the upper-right and lower-left 3x3 submatrix "blocks." BTW, the variables corresponding to each column and row of the correlation matrix are named by object 'selVars' in your script.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much!!! and the

Thank you so much!!! and the other 3x3 matrices are the phenotypic correlations right?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
Right.

Right.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much!!!

Thank you so much!!!