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.
cov2cor()?
So, what's your syntax involving
cov2cor()
, and what results does it give you?Log in or register to post comments
In reply to cov2cor()? by AdminRobK
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!
Log in or register to post comments
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.
Log in or register to post comments
In reply to Ahhh by mhunter
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
Log in or register to post comments
In reply to Thank you so much for your by JuanJMV
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.
Log in or register to post comments
In reply to Thank you so much for your by JuanJMV
my guess
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?
Log in or register to post comments
In reply to my guess by AdminRobK
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.
Log in or register to post comments
In reply to my guess by AdminRobK
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
Log in or register to post comments
In reply to Just one thing more I did the by JuanJMV
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.
Log in or register to post comments
In reply to between-twin, within-trait? by AdminRobK
Thank you so much!!! and the
Thank you so much!!! and the other 3x3 matrices are the phenotypic correlations right?
Log in or register to post comments
In reply to Thank you so much!!! and the by JuanJMV
Right.
Right.
Log in or register to post comments
In reply to Right. by AdminRobK
Thank you so much!!!
Thank you so much!!!
Log in or register to post comments