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!