# Intra Class Corelations

13 posts / 0 new
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. Offline Joined: 01/24/2014 - 12:15 cov2cor()? So, what's your syntax involving cov2cor(), and what results does it give you? 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! 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. 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 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. 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?

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.

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?

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.

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?

Offline
Joined: 01/24/2014 - 12:15
Right.

Right.

Offline
Joined: 07/20/2016 - 13:13
Thank you so much!!!

Thank you so much!!!