Hello,
First a thank you to this great forum! I am a beginner at twin model and have found it very helpful. I recently ran the Univariate Sex Limitation Model script http://www.vipbg.vcu.edu/HGEN619_2014/twinHet5AceCon.R At the end one runs a Homogeneity ACE model. My understanding is this should produce the same results as running the OpenMx ACE tutorial limited to just the young twins from the psych package, and excluding the opposite sex pairs. https://openmx.ssri.psu.edu/docs/openmx/latest/GeneticEpi_Matrix.html
However, when I compare the results while the estimates for A and E are similar, the estimate for C is not. Am I wrong in expecting these two models to match?
I have copied below the OpenMx tutorial script that I altered to use just the young MZ and DZ(same sex) twins from the psych package. I also included the results the two model compared.
Forgive me if this is has an obvious answer, I am very new to twin modeling. :)
Thank you!!
Brittany
------------------------------------------------------------------------------
MODEL COMPARISONS
Homogeneity ACE Model
round(HomoAceFit@output$estimate,4)
a11 c11 e11 meanf meanm
0.7510 0.0000 0.4083 21.3820 21.7235
OpenMx Univariate ACE Model
round(AceFit@output$estimate,4)
a11 c11 e11 mean
0.6036 0.4895 0.3989 21.4973
------------------------------------------------------------------------------
OPENMX TUTORIAL WITH YOUNG TWINS
------------------------------------------------------------------------------
require(OpenMx)
Select Variables for Analysis
Vars <- 'bmi' #this way the next expected number is either 1 or two for the selection below
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('T1tw1','T1tw2')
selVars
Select Data for Analysis
mzData <- subset(twinData, zyg==1 | zyg==2, selVars)
dzData <- subset(twinData, zyg==3 | zyg==4, selVars)
nrow(mzData) #569+273=842
nrow(dzData) #351+206=557
Set Starting Values
svMe <- 20
svPa <- 0.6 # start value for path coefficients (sqrt(variance/#ofpaths))
ACE Model
Matrices declared to store a, d, and e Path Coefficients
pathA <- mxMatrix( type="Full", nrow=nv, ncol=nv,
free=TRUE, values=svPa, label="a11", name="a" )
pathC <- mxMatrix( type="Full", nrow=nv, ncol=nv,
free=TRUE, values=svPa, label="c11", name="c" )
pathE <- mxMatrix( type="Full", nrow=nv, ncol=nv,
free=TRUE, values=svPa, label="e11", 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 total variances
covP <- mxAlgebra( expression=A+C+E, name="V" )
Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv,
free=TRUE, values=svMe, label="mean", name="expMean" )
free=TRUE, values=svMedz, label="meandz", name="expMeandz" )
covMZ <- mxAlgebra( expression=rbind( cbind(V, A+C),
cbind(A+C, V)), name="expCovMZ" )
covDZ <- mxAlgebra( expression=rbind( cbind(V, 0.5%x%A+ 0.25%x%C),
cbind(0.5%x%A+ 0.25%x%C , V)), name="expCovDZ" )
Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
Objective objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean",
dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean",
dimnames=selVars )
funML <- mxFitFunctionML()
Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP )
modelMZ <- mxModel( pars, meanG, covMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, meanG, covDZ, dataDZ, expDZ, funML, name="DZ" )
fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction") )
AceModel <- mxModel( "ACE", pars, modelMZ, modelDZ, fitML )
Run ADE model
AceFit <- mxRun(AceModel, intervals=T)
AceSum <- summary(AceFit)
AceSum
round(AceFit@output$estimate,4)
Thanks!! That is a huge help! :) Great catch that the tutorial was modeling ADE not ACE.
I re-ran it and it makes sense. I now realize that the difference between the tutorial and the final homogeneity model in the sex-limitation script is that the tutorial is only estimating 3paths and 1mean while the homogeneity model is doing 3 paths and 2 means.
I still have a general question. If my goal is to rule out sex differences for my estimate of A what is the most appropriate test? Do I want to compare QuanAceFit to HomoAceFit (both have different means for males and females)? OR should I be constraining the means and the pathways?
Thanks!
Hey!
There is an error in http://openmx.psyc.virginia.edu/docs/openmx/latest/GeneticEpi_Matrix.html . The DZ covariance structure is actually set up as if you were running and ADE model (just think of all C estimates you are getting as D). To change it back, switch:
covDZ <- mxAlgebra( expression=rbind( cbind(V, 0.5%x%A+ 0.25%x%C),
cbind(0.5%x%A+ 0.25%x%C , V)), name="expCovDZ" )
to:
covDZ <- mxAlgebra( expression=rbind( cbind(V, 0.5%x%A+ C),
cbind(0.5%x%A+C , V)), name="expCovDZ" )
It took me awhile to spot, I bet this has tripped up others too!
Brooke
Good eye, Brooke!