# Program: UnivSat&ACE_Bin.R (H. Maes) #require(foreign) #require(psych) require(OpenMx) citation("OpenMx") #stut <- read.table("twinstut.dat") data("twinstut") #read.table("twinstut.dat") stut <- twinstut str(stut) head(stut) # Females stutf <- subset(stut, sex=="female") stutf$sex <- stutf$sex[,drop=TRUE] stutf <- subset(stutf, zyg=="mz" | zyg=="dz") stutf$zyg <- stutf$zyg[,drop=TRUE] with(stutf, table(sex,zyg)) stutfwide <- reshape(stutf,idvar="tvparnr",timevar="nr", v.names=c("stutter","age","sex"), direction="wide") head(stutfwide) stutfw <- data.frame(cbind(stutfwide$zyg,stutfwide$stutter.1,stutfwide$stutter.2)) colnames(stutfw) <- c( 'zyg', 'stut1', 'stut2') str(stutfw) head(stutfw) Vars <- 'stut' nv <- 1 # number of variables per twin ntv <- nv*2 # number of variables per pair nthresh1<- 1 # number of max thresholds selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('stut1','stut2') #mz and dzss with(stutfw, table(stut1,stut2,zyg)) # make an mxFactor stutfw$stut1 <- mxFactor(stutfw$stut1, levels=c(1:2) ) stutfw$stut2 <- mxFactor(stutfw$stut2, levels=c(1:2) ) # Select Data for Analysis mzDataBin <- data.frame(subset(stutfw, zyg=="2", selVars)) dzDataBin <- data.frame(subset(stutfw, zyg=="1", selVars)) # Generate Descriptive Statistics table(mzDataBin$stut1,mzDataBin$stut2) table(dzDataBin$stut1,dzDataBin$stut2) # MxFactor enforces ordinal data to be sorted (for OpenMx) # the levels parameter is NOT optional, so that missing levels are not skipped, # and leaves the list of levels complete mzData <- mzDataBin dzData <- dzDataBin # Set Starting Values thVals <- 1.8 # start value for thresholds corValsM <-.8 # start value for MZ correlations corValsD <-.2 # start value for DZ correlations # Print Descriptive Statistics # --------------------------------------------------------------------- summary(mzData) summary(dzData) table(mzData$stut1, mzData$stut2 ) table(dzData$stut1, dzData$stut2) # ___________________________________ # PREPARE SATURATED MODEL # ___________________________________ # Matrices for expected Means & Thresholds (on liabilities) in MZ & DZ twin pairs meanG <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" ) threMZ <-mxMatrix(type="Full", nrow=1, ncol=ntv, free=TRUE, values=thVals, labels=c("tMZ1","tMZ2"), name="expThreMZ" ) threDZ <-mxMatrix(type="Full", nrow=1, ncol=ntv, free=TRUE, values=thVals, labels=c("tDZ1","tDZ2"), name="expThreDZ" ) corMZ <-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=corValsM, lbound=-.99, ubound=.99, labels=c("rMZ"), name="expCorMZ") corDZ <-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=corValsD, lbound=-.99, ubound=.99, labels=c("rDZ"), name="expCorDZ") # Data objects for Multiple Groups dataMZ <-mxData(mzData, type="raw") dataDZ <-mxData(dzData, type="raw") # Objective objects for Multiple Groups objMZ <-mxFIMLObjective( covariance="expCorMZ", means="expMean", dimnames=selVars, thresholds="expThreMZ" ) objDZ <-mxFIMLObjective( covariance="expCorDZ", means="expMean", dimnames=selVars, thresholds="expThreDZ" ) # Combine Groups groupMZ <-mxModel("MZ", corMZ, meanG, threMZ, dataMZ, objMZ ) groupDZ <-mxModel("DZ", corDZ, meanG, threDZ, dataDZ, objDZ ) minus2ll<-mxAlgebra( MZ.objective + DZ.objective, name="minus2sumloglikelihood" ) obj <-mxAlgebraObjective("minus2sumloglikelihood") ciCor <-mxCI(c('MZ.expCorMZ[2,1]', 'DZ.expCorDZ[2,1]')) ciThre <-mxCI( c('MZ.expThreMZ','DZ.expThreDZ' )) twinSatModel <-mxModel( "twinSat", minus2ll, obj, groupMZ, groupDZ, ciCor, ciThre ) # ----------------------------------------------------------------------- # RUN SATURATED MODEL (Tetrachoric correlations) # ----------------------------------------------------------------------- twinSatFit <- mxRun(twinSatModel, intervals=FALSE) twinSatSumm <- summary(twinSatFit) round(twinSatFit@output$estimate,4) twinSatSumm # Generate Saturated Model Output #-------------------------------------------------------------------------- rMZ<- twinSatFit$MZ.expCorMZ@values[2,1] rDZ<- twinSatFit$DZ.expCorDZ@values[2,1] tMZ<- twinSatFit$MZ.expThreMZ@values tDZ<- twinSatFit$DZ.expThreDZ@values twinSatLLL <-twinSatFit@output$Minus2LogLikelihood twinSatAIC <-twinSatSumm$AIC twinSatOS <-twinSatSumm$observedStatistics twinSatDF <-twinSatSumm$degreesOfFreedom twinSatNP <-length(twinSatSumm$parameters[[1]]) # Use Helper Functions source("GenEpiHelperFunctions.R") expectedMeansCovariances(twinSatFit) tableFitStatistics(twinSatFit) # ------------------------------------------------------------------------- # RUN SUBMODELS # ------------------------------------------------------------------------- # SubModel 1: constraining Thresholds across Twin 1 and Twin 2 within zyg group to be equal # ---------------------------------------------------------------------------------------------------------------------------------------- eqThresholdsTwinModel <- twinSatFit eqThresholdsTwinModel <- omxSetParameters( eqThresholdsTwinModel, label="tMZ1", free=TRUE, values=thVals, newlabels='tMZ' ) eqThresholdsTwinModel <- omxSetParameters( eqThresholdsTwinModel, label="tMZ2", free=TRUE, values=thVals, newlabels='tMZ' ) eqThresholdsTwinModel <- omxSetParameters( eqThresholdsTwinModel, label="tDZ1", free=TRUE, values=thVals, newlabels='tDZ' ) eqThresholdsTwinModel <- omxSetParameters( eqThresholdsTwinModel, label="tDZ2", free=TRUE, values=thVals, newlabels='tDZ' ) eqThresholdsTwinFit <- mxRun( eqThresholdsTwinModel, intervals=T ) eqThresholdsTwinSumm <- summary( eqThresholdsTwinFit ) eqThresholdsTwinLLL <- eqThresholdsTwinFit@output$Minus2LogLikelihood tableFitStatistics(twinSatFit, eqThresholdsTwinFit) eqThresholdsTwinSumm # SubModel 2: constraining Thresholds across Twin 1 and Twin 2 and across zyg group to be equal # ---------------------------------------------------------------------------------------------------------------------------------------------- eqThresholdsZygModel <- eqThresholdsTwinModel eqThresholdsZygModel <- omxSetParameters( eqThresholdsZygModel, label="tMZ", free=TRUE, values=thVals, newlabels='tZ' ) eqThresholdsZygModel <- omxSetParameters( eqThresholdsZygModel, label="tDZ", free=TRUE, values=thVals, newlabels='tZ' ) eqThresholdsZygFit <- mxRun( eqThresholdsZygModel, intervals=F ) eqThresholdsZygSumm <- summary( eqThresholdsZygFit ) eqThresholdsZygLLL <- eqThresholdsZygFit@output$Minus2LogLikelihood tableFitStatistics(eqThresholdsTwinFit, eqThresholdsZygFit) eqThresholdsZygSumm # Print Comparative Fit Statistics for Saturated Models # ----------------------------------------------------------------------------- SatNested <- list(eqThresholdsTwinFit, eqThresholdsZygFit) tableFitStatistics(twinSatFit, SatNested) # ___________________________________ # PREPARE GENETIC MODEL # ___________________________________ # ACE Model with one overall Threshold # Matrices to store a, c, and e Path Coefficients pathA <-mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a11", name="a" ) pathC <-mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c11", name="c" ) pathE <-mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e11", name="e" ) # Algebra to generate Matrices 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" ) # Matrices for expected Means & Thresholds (on liabilities) meanG <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" ) threT <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=thVals, label="thre", name="expThre" ) # Algebra to compute Total Variance covP <-mxAlgebra( expression=A+C+E, name="V" ) # Algebra to generate Matrices to hold Parameter Estimates and Derived Variance Components rowVars <-rep('vars',nv) colVars <-rep(c('A','C','E','SA','SC','SE'),each=nv) estVars <-mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="Vars", dimnames=list(rowVars,colVars)) # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covMZ <-mxAlgebra( expression= rbind( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) covDZ <-mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) # Constraint on variance of the liability of Binary variables (assumed to have a SND) matUnv <-mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv1" ) var1 <-mxConstraint( expression=diag2vec(V)==Unv1, name="Var1" ) # Data objects for Multiple Groups dataMZ <-mxData( observed=mzData, type="raw" ) dataDZ <-mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <-mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="expThre" ) objDZ <-mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="expThre" ) # Combine Groups pars <-list( pathA, pathC, pathE, covA, covC, covE, covP, matUnv, estVars ) modelMZ <-mxModel( pars, meanG, threT, covMZ, dataMZ, objMZ, name="MZ" ) modelDZ <-mxModel( pars, meanG, threT, covDZ, dataDZ, objDZ, name="DZ" ) minus2ll<-mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <-mxAlgebraObjective( "m2LL" ) ciVC <-mxCI(c('A', 'C', 'E')) AceModel <-mxModel( "ACE", pars, var1, modelMZ, modelDZ, minus2ll, obj, ciVC) # ------------------------------------------------- # RUN GENETIC MODEL # ------------------------------------------------- # Run ACE Model AceFit <- mxRun(AceModel, intervals=T ) AceSumm <- summary(AceFit) AceSumm round(AceFit@output$estimate,4) round(AceFit$Vars@result,4) # Generate ACE Model Output (using helper function GenEpiHelperFunctions.R) parameterSpecifications(AceFit) expectedMeansCovariances(AceFit) tableFitStatistics(AceFit) ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covs_A","covs_C","covs_E","Var","prop_A","prop_C","prop_E") formatOutputMatrices(AceFit,ACEcovMatrices,ACEcovLabels,Vars,4) # ----------------------------------------------------------- # RUN SUBMODELS # ----------------------------------------------------------- # Fit AE model # ----------------------------------------------------------------- AeModel <- mxModel( AceFit, name="AE") AeModel <- omxSetParameters( AeModel, labels="c11", free=FALSE, values=0 ) AeFit <- mxRun(AeModel, intervals=T) AeSumm <- summary(AeFit) AeSumm round(AeFit@output$estimate,4) round(AeFit$Vars@result,4) # Run CE model # ------------------------------------------------------------------ CeModel <- mxModel( AceFit, name="CE") CeModel <- omxSetParameters( CeModel, labels="a11", free=FALSE, values=0 ) CeFit <- mxRun(CeModel) round(CeFit@output$estimate,4) round(CeFit$Vars@result,4) # ------------------------------------------------------------------------------ # Run E model # ------------------------------------------------------------------ eModel <- mxModel( AeFit, name="E") eModel <- omxSetParameters( eModel, labels="a11", free=FALSE, values=0 ) eFit <- mxRun(eModel) round(eFit@output$estimate,4) round(eFit$Vars@result,4) # ------------------------------------------------------------------------------ # Print Comparative Fit Statistics AceNested <- list(AeFit, CeFit, eFit) tableFitStatistics(AceFit,AceNested) round(rbind(AceFit@output$estimate,AeFit@output$estimate,CeFit@output$estimate,eFit@output$estimate),4) round(rbind(AceFit$Vars@result,AeFit$Vars@result,CeFit$Vars@result,eFit$Vars@result),4)