Step 1: Testing basic assumptions (equal means & variances for twin 1/twin 2 and MZ/DZ pairs)
inooradd
Joined: 03/11/2020
I am following the steps laid out in [this](https://vipbg.vcu.edu/media/course/HGEN619_2015/TwinModels.pdf ) openmx presentation and could use help with step 1.
1. Use data to test basic assumptions (equal means & variances for twin 1/twin 2 and MZ/DZ pairs)
- I am working with MZ, DZ, and siblings. If I test for differences between DZ twins and siblings, and they are not significantly different, do I just code siblings as DZ twins and proceed to compare MZ and DZ twins? I know when using ACE/ADE you just include or omit the T. How would you code the DZ twins and sibs as the same in this instance? Is the following correct, where I code sibs as dzData?
<%
conData <- twindata2
mzData <- subset(conData, zyg==1, c(selVars, covVars))
dzData <- subset(conData, zyg==2 | zyg==3, c(selVars, covVars))
sibData <- subset(conData, zyg==4 | zyg==5, c(selVars, covVars))
%>
becomes
<%
conData <- twindata2
mzData <- subset(conData, zyg==1, c(selVars, covVars))
dzData <- subset(conData, zyg==2 | zyg==3 | zyg==4 | zyg==5, c(selVars, covVars))
%>
- If this is the case, how would I proceed with multiple siblings? Would I just incorporate them into the syntax?
- Additionally, I encountered an error constraining the means and variances across twin order for MZ and DZ twins.
<%
> # Constrain expected Means and Variances to be equal across twin order
> eqMVarsTwinModel <- mxModel(eqMeansTwinFit, name="eqMVarsTwin" )
> eqMVarsTwinModel <- omxSetParameters( eqMVarsTwinModel, label=c("v1MZ","v2MZ"), free=TRUE, values=svMe, newlabels='vMZ' )
> eqMVarsTwinModel <- omxSetParameters( eqMVarsTwinModel, label=c("v1DZ","v2DZ"), free=TRUE, values=svMe, newlabels='vDZ' )
>
> eqMVarsTwinFit <- mxTryHard( eqMVarsTwinModel, intervals=F )
All fit attempts resulted in errors - check starting values or model specification
> subs <- list(eqMeansTwinFit, eqMVarsTwinFit)
> mxCompare(twinSatFit, subs)
base comparison ep minus2LL df AIC diffLL diffdf p
1 twinSat
2 twinSat eqMeansTwin 12 3811.663 1380 1051.663 2.013412 2 0.3654207
3 twinSat eqMVarsTwin 0 NA 1392 NA NA 14 NA
Warning message:
In collectStatistics1(otherStats, ref, other, bootPair) :
MxModel 'eqMVarsTwin' has no 'fitUnits' element in its output slot; has it been run?
%>
Below is my syntax.
<%
# Select Variables for Analysis
vars <- 'ideo' # 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
covVars <- c('age1',"age2", "sex1","sex2")
nc <- 1 # number of covariates
# Select Data for Analysis
conData <- twindata2
mzData <- subset(conData, zyg==1, c(selVars, covVars))
dzData <- subset(conData, zyg==2 | zyg==3, c(selVars, covVars))
sibData <- subset(conData, zyg==4 | zyg==5, c(selVars, covVars))
# Set Starting Values
svMe <- 0 # start value for means
svPa <- .6 # start value for path coefficient
svPe <- .6 # start value for path coefficient for e
lbPa <- .0000000000000000000000000000000001 # start value for lower bounds
svPas <- diag(svPa,ntv,ntv) # assign start values to diagonal1 of matrix
lbPas <- diag(lbPa,ntv,ntv) # assign lower bounds values to diagonal of matrix
laMeMZ <- c("m1MZ","m2MZ") # labels for means for MZ twins
laMeDZ <- c("m1DZ","m2DZ") # labels for means for DZ twins
laVaMZ <- c("v1MZ","c21MZ","v2MZ") # labels for (co)variances for MZ twins
laVaDZ <- c("v1DZ","c21DZ","v2DZ") # labels for (co)variances for DZ twins
# ------------------------------------------------------------------------------
# PREPARE SATURATED MODEL
# Create Matrices for Covariates and linear Regression Coefficients
defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1", "data.age2"),name="Age")
pathB1 <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=.01, label=c("b11", "b12"), 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=2, free=TRUE, values=.01, label=c("b21", "b22"), name= "b2" )
# Saturated Model
# Algebra for expected Mean Matrices in MZ & DZ twins
meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=laMeMZ, name="meanMZ" )
meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=laMeDZ, name="meanDZ" )
expMeanMZ <- mxAlgebra( expression= meanMZ + b1*Age + b2*Sex, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + b1*Age + b2*Sex, name="expMeanDZ" )
# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
expCovMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svPas, lbound=lbPas, labels=laVaMZ, name="expCovMZ" )
expCovDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svPas, lbound=lbPas, labels=laVaDZ, name="expCovDZ" )
# Create 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="expMeanMZ", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars )
funML <- mxFitFunctionML()
# Combine Groups
modelMZ <- mxModel( "MZ", defAge, defSex, pathB1, pathB2, meanMZ, expMeanMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( "DZ", defAge, defSex, pathB1, pathB2, meanDZ, expMeanDZ, expCovDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
ciCov <- mxCI( c('MZ.expCovMZ','DZ.expCovDZ' ))
ciMean <- mxCI( c('MZ.expMeanMZ','DZ.expMeanDZ' ))
twinSatModel <- mxModel( "twinSat", modelMZ, modelDZ, multi, ciCov, ciMean )
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# RUN SATURATED MODELS
# Run Saturated Model
twinSatFit <- mxRun( twinSatModel, intervals=F )
twinSatSum <- summary( twinSatFit )
twinSatSum
# Generate Saturated Model Output
twinSatFit$MZ.expMeanMZ$result
twinSatFit$DZ.expMeanDZ$result
twinSatFit$MZ.expCovMZ$values
twinSatFit$DZ.expCovDZ$values
# Print Goodness-of-Fit Statistics
twinSatSum$observedStatistics
length(twinSatSum$parameters[[1]])
twinSatFit$output$Minus2LogLikelihood
twinSatSum$degreesOfFreedom
twinSatSum$AIC
round(twinSatFit$output$estimate,4)
# Test Significance of Covariate
# Copy model, provide new name
testCovModel <- mxModel(twinSatFit, name="testCov")
# Change parameter by changing attributes for label
testCovModel <- omxSetParameters( testCovModel, label="b11", free=FALSE, values=0 )
# Fit Nested Model
testCovFit <- mxRun(testCovModel)
# Compare Nested Model with 'Full' Model
mxCompare(twinSatFit, testCovFit)
# ------------------------------------------------------------------------------
# RUN SUBMODELS
# Constrain expected Means to be equal across twin order
eqMeansTwinModel <- mxModel(twinSatFit, name="eqMeansTwin" )
eqMeansTwinModel <- omxSetParameters( eqMeansTwinModel, label=c("m1MZ","m2MZ"), free=TRUE, values=svMe, newlabels='mMZ' )
eqMeansTwinModel <- omxSetParameters( eqMeansTwinModel, label=c("m1DZ","m2DZ"), free=TRUE, values=svMe, newlabels='mDZ' )
eqMeansTwinFit <- mxRun( eqMeansTwinModel, intervals=F )
eqMeansTwinSum <- summary( eqMeansTwinFit )
eqMeansTwinLLL <- eqMeansTwinFit$output$Minus2LogLikelihood
twinSatLLL <- twinSatFit$output$Minus2LogLikelihood
chi2Sat_eqM <- eqMeansTwinLLL-twinSatLLL
pSat_eqM <- pchisq( chi2Sat_eqM, lower.tail=F, 2 )
chi2Sat_eqM; pSat_eqM
mxCompare(twinSatFit, eqMeansTwinFit)
# Constrain expected Means and Variances to be equal across twin order
eqMVarsTwinModel <- mxModel(eqMeansTwinFit, name="eqMVarsTwin" )
eqMVarsTwinModel <- omxSetParameters( eqMVarsTwinModel, label=c("v1MZ","v2MZ"), free=TRUE, values=svMe, newlabels='vMZ' )
eqMVarsTwinModel <- omxSetParameters( eqMVarsTwinModel, label=c("v1DZ","v2DZ"), free=TRUE, values=svMe, newlabels='vDZ' )
eqMVarsTwinFit <- mxRun( eqMVarsTwinModel, intervals=F )
subs <- list(eqMeansTwinFit, eqMVarsTwinFit)
mxCompare(twinSatFit, subs)
# Constrain expected Means and Variances to be equal across twin order and zygosity
eqMVarsZygModel <- mxModel(eqMVarsTwinModel, name="eqMVarsZyg" )
eqMVarsZygModel <- omxSetParameters( eqMVarsZygModel, label=c("mMZ","mDZ"), free=TRUE, values=svMe, newlabels='mZ' )
eqMVarsZygModel <- omxSetParameters( eqMVarsZygModel, label=c("vMZ","vDZ"), free=TRUE, values=svMe, newlabels='vZ' )
eqMVarsZygFit <- mxRun( eqMVarsZygModel, intervals=F )
mxCompare(eqMVarsTwinFit, eqMVarsZygFit)
# ------------------------------------------------------------------------------
# Print Comparative Fit Statistics
SatNested <- list(eqMeansTwinFit, eqMVarsTwinFit, eqMVarsZygFit)
mxCompare(twinSatFit, SatNested)
%>
Any idea what might be causing this issue?