Hello!
I am following the steps laid out in this openmx presentation and could use help with step 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?
Printer-friendly version