# Step 1: Testing basic assumptions (equal means & variances for twin 1/twin 2 and MZ/DZ pairs)

1 post / 0 new
Offline
Joined: 03/11/2020 - 11:42
Step 1: Testing basic assumptions (equal means & variances for twin 1/twin 2 and MZ/DZ pairs)

Hello!

I am following the steps laid out in this 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        <NA> 14 3809.650 1378 1053.650       NA     NA        NA
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?