# Load Libraries & Options library(OpenMx) library(psych) library('frequency') library("gee") library(frequency) library(dplyr) library(haven) source("miFunctions3.R") # I normally read in the data here # Select Data for Analysis mzfData <- subset(WB_pers, zyg5==3, c(selVars,covVars)) dzfData <- subset(WB_pers, zyg5==4, c(selVars,covVars)) mzmData <- subset(WB_pers, zyg5==1, c(selVars,covVars)) dzmData <- subset(WB_pers, zyg5==2, c(selVars,covVars)) dzoData <- subset(WB_pers, zyg5==5, c(selVars,covVars)) mzfsibData$zyg5<-3 dzfsibData$zyg5<-4 mzmsibData$zyg5<-1 dzmsibData$zyg5<-2 dzosibData$zyg5<-5 WB_pers_final<-rbind(mzfsibData,dzfsibData,mzmsibData,dzmsibData,dzosibData) # Select Continuous Variables varsc <- c('contvar1','contvar2','contvar3','contvar4','contvar5', 'contvar6','contvar7','contvar8') # list of variables names nvc <- 8 # number of variables ntvc <- nvc*4 # number of total continuous variables conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc),rep(3,nvc),rep(4,nvc)),sep="_") # Select Ordinal Variables nth <- 2 # number of thresholds varso <- c('ordvar1','ordvar2') # list of ordinal variables names nvo <- 2 # number of ordinal variables ntvo <- nvo*4 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo),rep(3,nvo),rep(4,nvo)),sep="_") ordData <- WB_pers_final for (i in c('ordvar2_1','ordvar2_2','ordvar2_3','ordvar2_4')) { ordData[,i] <- ifelse(ordData[,i]>2.5,ifelse(ordData[,i]>3.5,2,1),0) } for (i in c('ordvar1_1','ordvar1_2','ordvar1_3','ordvar1_4')) { ordData[,i] <- ifelse(ordData[,i]>0.50,ifelse(ordData[,i]>1.10,2,1),0) } # Select Variables for Analysis vars <- c('ordvar1','ordvar2','contvar1','contvar2','contvar3','contvar4','contvar5', 'contvar6','contvar7','contvar8') # list of variables names nv <- nvc+nvo # number of variables ntv <- nv*4 # number of total variables oc <- c(1,1,rep(0,8)) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv),rep(3,nv),rep(4,nv)),sep="_") varscov <- c('age') # list of covariate names covVars <- paste(varscov,c(rep(1,1),rep(2,1),rep(3,1), rep(4,1)),sep="_") # Select Data for Analysis mzfsibData <- subset(ordData, zyg5==3, c(selVars,covVars)) dzfsibData <- subset(ordData, zyg5==4, c(selVars,covVars)) mzmsibData <- subset(ordData, zyg5==1, c(selVars,covVars)) dzmsibData <- subset(ordData, zyg5==2, c(selVars,covVars)) dzosibData <- subset(ordData, zyg5==5, c(selVars,covVars)) mzmsibDataF <- cbind(mxFactor( x=mzmsibData[,c('ordvar1_1','ordvar2_1')], levels=c(0:nth)), mzmsibData[,c('contvar1_1','contvar2_1','contvar3_1','contvar4_1','contvar5_1','contvar6_1','contvar7_1','contvar8_1')], mxFactor( x=mzmsibData[,c('ordvar1_2','ordvar2_2')], levels=c(0:nth)), mzmsibData[,c('contvar1_2','contvar2_2','contvar3_2','contvar4_2','contvar5_2','contvar6_2','contvar7_2','contvar8_2')], mxFactor( x=mzmsibData[,c('ordvar1_3','ordvar2_3')], levels=c(0:nth)), mzmsibData[,c('contvar1_3','contvar2_3','contvar3_3','contvar4_3','contvar5_3','contvar6_3','contvar7_3','contvar8_3')], mxFactor( x=mzmsibData[,c('ordvar1_4','ordvar2_4')], levels=c(0:nth)), mzmsibData[,c('contvar1_4','contvar2_4','contvar3_4','contvar4_4','contvar5_4','contvar6_4','contvar7_4','contvar8_4')], mzmsibData[,covVars]) mzfsibDataF <- cbind(mxFactor( x=mzfsibData[,c('ordvar1_1','ordvar2_1')], levels=c(0:nth)), mzfsibData[,c('contvar1_1','contvar2_1','contvar3_1','contvar4_1','contvar5_1','contvar6_1','contvar7_1','contvar8_1')], mxFactor( x=mzfsibData[,c('ordvar1_2','ordvar2_2')], levels=c(0:nth)), mzfsibData[,c('contvar1_2','contvar2_2','contvar3_2','contvar4_2','contvar5_2','contvar6_2','contvar7_2','contvar8_2')], mxFactor( x=mzfsibData[,c('ordvar1_3','ordvar2_3')], levels=c(0:nth)), mzfsibData[,c('contvar1_3','contvar2_3','contvar3_3','contvar4_3','contvar5_3','contvar6_3','contvar7_3','contvar8_3')], mxFactor( x=mzfsibData[,c('ordvar1_4','ordvar2_4')], levels=c(0:nth)), mzfsibData[,c('contvar1_4','contvar2_4','contvar3_4','contvar4_4','contvar5_4','contvar6_4','contvar7_4','contvar8_4')], mzfsibData[,covVars]) dzmsibDataF <- cbind(mxFactor( x=dzmsibData[,c('ordvar1_1','ordvar2_1')], levels=c(0:nth)), dzmsibData[,c('contvar1_1','contvar2_1','contvar3_1','contvar4_1','contvar5_1','contvar6_1','contvar7_1','contvar8_1')], mxFactor( x=dzmsibData[,c('ordvar1_2','ordvar2_2')], levels=c(0:nth)), dzmsibData[,c('contvar1_2','contvar2_2','contvar3_2','contvar4_2','contvar5_2','contvar6_2','contvar7_2','contvar8_2')], mxFactor( x=dzmsibData[,c('ordvar1_3','ordvar2_3')], levels=c(0:nth)), dzmsibData[,c('contvar1_3','contvar2_3','contvar3_3','contvar4_3','contvar5_3','contvar6_3','contvar7_3','contvar8_3')], mxFactor( x=dzmsibData[,c('ordvar1_4','ordvar2_4')], levels=c(0:nth)), dzmsibData[,c('contvar1_4','contvar2_4','contvar3_4','contvar4_4','contvar5_4','contvar6_4','contvar7_4','contvar8_4')], dzmsibData[,covVars]) dzfsibDataF <- cbind(mxFactor( x=dzfsibData[,c('ordvar1_1','ordvar2_1')], levels=c(0:nth)), dzfsibData[,c('contvar1_1','contvar2_1','contvar3_1','contvar4_1','contvar5_1','contvar6_1','contvar7_1','contvar8_1')], mxFactor( x=dzfsibData[,c('ordvar1_2','ordvar2_2')], levels=c(0:nth)), dzfsibData[,c('contvar1_2','contvar2_2','contvar3_2','contvar4_2','contvar5_2','contvar6_2','contvar7_2','contvar8_2')], mxFactor( x=dzfsibData[,c('ordvar1_3','ordvar2_3')], levels=c(0:nth)), dzfsibData[,c('contvar1_3','contvar2_3','contvar3_3','contvar4_3','contvar5_3','contvar6_3','contvar7_3','contvar8_3')], mxFactor( x=dzfsibData[,c('ordvar1_4','ordvar2_4')], levels=c(0:nth)), dzfsibData[,c('contvar1_4','contvar2_4','contvar3_4','contvar4_4','contvar5_4','contvar6_4','contvar7_4','contvar8_4')], dzfsibData[,covVars]) dzosibDataF <- cbind(mxFactor( x=dzosibData[,c('ordvar1_1','ordvar2_1')], levels=c(0:nth)), dzosibData[,c('contvar1_1','contvar2_1','contvar3_1','contvar4_1','contvar5_1','contvar6_1','contvar7_1','contvar8_1')], mxFactor( x=dzosibData[,c('ordvar1_2','ordvar2_2')], levels=c(0:nth)), dzosibData[,c('contvar1_2','contvar2_2','contvar3_2','contvar4_2','contvar5_2','contvar6_2','contvar7_2','contvar8_2')], mxFactor( x=dzosibData[,c('ordvar1_3','ordvar2_3')], levels=c(0:nth)), dzosibData[,c('contvar1_3','contvar2_3','contvar3_3','contvar4_3','contvar5_3','contvar6_3','contvar7_3','contvar8_3')], mxFactor( x=dzosibData[,c('ordvar1_4','ordvar2_4')], levels=c(0:nth)), dzosibData[,c('contvar1_4','contvar2_4','contvar3_4','contvar4_4','contvar5_4','contvar6_4','contvar7_4','contvar8_4')], dzosibData[,covVars]) # Set Starting Values frMV <- c(FALSE,FALSE,rep(TRUE,8)) # free status for variables frCvD <- valDiagLU(frMV,T,T,ntv) # free status for diagonal, lower & upper elements of covariance matrix frCv <- matrix(as.logical(frCvD),ntv) svMe <- c(0,0,3,6,5,8,7,4,6,3) # start value for means svVa <- c(1,1,rep(0.9,8)) # start value for variances lbVa <- -100 # lower bound for variances svLTh <- c(-0.2,0.5) # start value for first threshold svITh <- c(0,0.7) # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=ntvo) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),ntvo),nrow=nth,ncol=ntvo) # lower bound for thresholds # Create Labels labMeMZf <- labVars("meanMZf",selVars) labMeDZf <- labVars("meanDZf",selVars) labMeMZm <- labVars("meanMZm",selVars) labMeDZm <- labVars("meanDZm",selVars) labMeDZo <- labVars("meanDZo",selVars) labMeMZ <- labVars("meanMZ",selVars) labMeDZ <- labVars("meanDZ",selVars) labMeZ <- labVars("meanZ",selVars) labCvMZf <- labLower("covMZf",ntv) labCvDZf <- labLower("covDZf",ntv) labCvMZm <- labLower("covMZm",ntv) labCvDZm <- labLower("covDZm",ntv) labCvDZo <- labLower("covDZo",ntv) labCvSIB <- labLower("covSIB",ntv) labCvMZ <- labLower("covMZ",ntv) labCvDZ <- labLower("covDZ",ntv) labCvZ <- labLower("covZ",ntv) labVaMZf <- labDiag("covMZf",ntv) labVaDZf <- labDiag("covDZf",ntv) labVaMZm <- labDiag("covMZm",ntv) labVaDZm <- labDiag("covDZm",ntv) labVaDZo <- labDiag("covDZo",ntv) labVaSIB <- labDiag("covSIB",ntv) labVaMZ <- labDiag("covMZ",ntv) labVaDZ <- labDiag("covDZ",ntv) labVaZ <- labDiag("covZ",ntv) labThMZm <- labTh("MZm",ordVars,nth) labThMZf <- labTh("MZf",ordVars,nth) labThDZm <- labTh("DZm",ordVars,nth) labThDZf <- labTh("DZf",ordVars,nth) labThDZo <- labTh("DZo",ordVars,nth) labThMZ <- labTh("MZ",ordVars,nth) labThDZ <- labTh("DZ",ordVars,nth) labThZ <- labTh("Z",ordVars,nth) # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=4, free=FALSE, labels=c("data.age_1","data.age_2","data.age_3","data.age_4"), name="age" ) pathB <- mxMatrix( type="Full", nrow=1, ncol=10, free=frMV, values=.01, label=c("beta11","beta12","beta13","beta14","beta15", "beta16","beta17","beta18","beta19","beta110"), name="b" ) pathB_Ord <- mxMatrix( type="Full", nrow=nth, ncol=nvo, free=TRUE, values= 0, labels=c("betaOrdv1","betaOrdv1","betaOrdv2", "betaOrdv2"), name="bOrd") # Create Algebra for expected Mean Matrices meanMZf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZf, name="meanMZf" ) meanDZf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZf, name="meanDZf" ) meanMZm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZm, name="meanMZm" ) meanDZm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZm, name="meanDZm" ) meanDZo <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZo, name="meanDZo" ) expMeanMZf <- mxAlgebra( expression= meanMZf + age%x%b, name="expMeanMZf" ) expMeanDZf <- mxAlgebra( expression= meanDZf + age%x%b, name="expMeanDZf" ) expMeanMZm <- mxAlgebra( expression= meanMZm + age%x%b, name="expMeanMZm" ) expMeanDZm <- mxAlgebra( expression= meanDZm + age%x%b, name="expMeanDZm" ) expMeanDZo <- mxAlgebra( expression= meanDZo + age%x%b, name="expMeanDZo" ) #expMeanSIB <- mxAlgebra( expression= meanDZo + b*age, name="expMeanSIB" ) # Create thresholds thinMZm <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labThMZm, name="thinMZm" ) thinMZf <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labThMZf, name="thinMZf" ) thinDZm <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labThDZm, name="thinDZm" ) thinDZf <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labThDZf, name="thinDZf" ) thinDZo <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labThDZo, name="thinDZo" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threMZm <- mxAlgebra( expression= (inc %*% thinMZm) + (age%x%bOrd), name="threMZm" ) threMZf <- mxAlgebra( expression= (inc %*% thinMZf) + (age%x%bOrd), name="threMZf" ) threDZm <- mxAlgebra( expression= (inc %*% thinDZm) + (age%x%bOrd), name="threDZm" ) threDZf <- mxAlgebra( expression= (inc %*% thinDZf) + (age%x%bOrd), name="threDZf" ) threDZo <- mxAlgebra( expression= (inc %*% thinDZo) + (age%x%bOrd), name="threDZo" ) covMZf <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=valDiag(svVa,ntv), labels=labCvMZf, name="covMZf" ) covDZf <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=valDiag(svVa,ntv), labels=labCvDZf, name="covDZf" ) covMZm <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=valDiag(svVa,ntv), labels=labCvMZm, name="covMZm" ) covDZm <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=valDiag(svVa,ntv), labels=labCvDZm, name="covDZm" ) covDZo <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=valDiag(svVa,ntv), labels=labCvDZo, name="covDZo" ) # Create correlation matrices cov2cor corMZf <- mxAlgebra(cov2cor(covMZf), name="corMZf") corDZf <- mxAlgebra(cov2cor(covDZf), name="corDZf") corMZm <- mxAlgebra(cov2cor(covMZm), name="corMZm") corDZm <- mxAlgebra(cov2cor(covDZm), name="corDZm") corDZo <- mxAlgebra(cov2cor(covDZo), name="corDZo") corSIB <- mxAlgebra(cov2cor(covSIB), name="corDZo") # Create Data Objects for Multiple Groups dataMZf <- mxData( observed=mzfsibDataF, type="raw" ) dataDZf <- mxData( observed=dzfsibDataF, type="raw" ) dataMZm <- mxData( observed=mzmsibDataF, type="raw" ) dataDZm <- mxData( observed=dzmsibDataF, type="raw" ) dataDZo <- mxData( observed=dzosibDataF, type="raw" ) # Create Expectation Objects for Multiple Groups expMZf <- mxExpectationNormal( covariance="covMZf", means="expMeanMZf", dimnames=selVars, thresholds="threMZf", threshnames=c('ordvar1_1','ordvar2_1','ordvar1_2','ordvar2_2','ordvar1_3','ordvar2_3','ordvar1_4','ordvar2_4') ) expDZf <- mxExpectationNormal( covariance="covDZf", means="expMeanDZf", dimnames=selVars, thresholds="threDZf", threshnames=c('ordvar1_1','ordvar2_1','ordvar1_2','ordvar2_2','ordvar1_3','ordvar2_3','ordvar1_4','ordvar2_4') ) expMZm <- mxExpectationNormal( covariance="covMZm", means="expMeanMZm", dimnames=selVars, thresholds="threMZm", threshnames=c('ordvar1_1','ordvar2_1','ordvar1_2','ordvar2_2','ordvar1_3','ordvar2_3','ordvar1_4','ordvar2_4') ) expDZm <- mxExpectationNormal( covariance="covDZm", means="expMeanDZm", dimnames=selVars, thresholds="threDZm", threshnames=c('ordvar1_1','ordvar2_1','ordvar1_2','ordvar2_2','ordvar1_3','ordvar2_3','ordvar1_4','ordvar2_4') ) expDZo <- mxExpectationNormal( covariance="covDZo", means="expMeanDZo", dimnames=selVars, thresholds="threDZo", threshnames=c('ordvar1_1','ordvar2_1','ordvar1_2','ordvar2_2','ordvar1_3','ordvar2_3','ordvar1_4','ordvar2_4') ) #expSIB <- mxExpectationNormal( covariance="covSIB", means="expMeanSIB", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups modelMZf <- mxModel( meanMZf, defAge, pathB, pathB_Ord, covMZf, expMeanMZf, corMZf, dataMZf, expMZf, funML, thinMZf, inc,threMZf, name="MZf" ) modelDZf <- mxModel( meanDZf, defAge, pathB, pathB_Ord, covDZf, expMeanDZf, corDZf, dataDZf, expDZf, funML, thinDZf, inc,threDZf, name="DZf" ) modelMZm <- mxModel( meanMZm, defAge, pathB, pathB_Ord, covMZm, expMeanMZm, corMZm, dataMZm, expMZm, funML, thinMZm, inc,threMZm, name="MZm" ) modelDZm <- mxModel( meanDZm, defAge, pathB, pathB_Ord, covDZm, expMeanDZm, corDZm, dataDZm, expDZm, funML, thinDZm, inc,threDZm, name="DZm" ) modelDZo <- mxModel( meanDZo, defAge, pathB, pathB_Ord, covDZo, expMeanDZo, corDZo, dataDZo, expDZo, funML, thinDZo, inc,threDZo, name="DZo" ) multi <- mxFitFunctionMultigroup( c("MZf","DZf","MZm","DZm","DZo") ) # Create Confidence Interval Objects ciMean <- mxCI( c('MZf.meanMZf','DZf.meanDZf','MZm.meanMZm','DZm.meanDZm','DZo.meanDZo' )) ciCor <- mxCI( c("MZf.corMZf", "MZm.corMZm","DZf.corDZf","DZm.corDZm","DZo.corDZo")) ciVar <- mxCI( c("MZf.covMZf", "MZm.covMZm","DZf.covDZf","DZm.covDZm","DZo.covDZo")) # Build Saturated Model with Confidence Intervals modelSAT <- mxModel( "twoSATc", modelMZf, modelDZf,modelMZm, modelDZm, modelDZo, multi, ciMean, ciCor, ciVar ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run Saturated Model tick <- proc.time()[3] gc() modelSAT<-mxOption(model= modelSAT, key="Number of Threads", value= omxDetectCores() ) modelSAT <- mxOption(modelSAT, "Standard Errors", "No") modelSAT <- mxOption(modelSAT, "Calculate Hessian", "No") mxOption(NULL,"Default optimizer","CSOLNP") fitSAT2 <- mxTryHardOrdinal(modelSAT, extraTries = 30,bestInitsOutput=T, intervals=FALSE) tock <- proc.time()[3] tock-tick fitSAT2$output$matrices #sumSAT <- summary( fitSAT2 ) #Print correlation matrices #round(fitSAT2$output$algebras$MZf.corMZf,3) #round(fitSAT2$output$algebras$DZf.corDZf,3) #round(fitSAT2$output$algebras$MZm.corMZm,3) #round(fitSAT2$output$algebras$DZm.corDZm,3) #round(fitSAT2$output$algebras$DZo.corDZo,3)