Intra Class Corelations

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