rm(list=ls()) require(OpenMx) require(psych) library(foreign) require(readstata13) require(ggplot) #mxOption(NULL, "Default optimizer", "NPSOL") mxOption( NULL, 'Default optimizer' ,'SLSQP' ) #mxOption(NULL, "Default optimizer", "CSOLNP") mxOption( NULL, "mvnRelEps", value= mxOption(NULL, key="mvnRelEps")/5) # ----------------------------------------------------------------------- # PREPARE DATA # Read Data twinData <-read.dta13("20477pair-coronary.dta") dim(twinData) describe(twinData, skew=F) twinData <- twinData [!is.na(twinData $coronary1)&!is.na(twinData $coronary2),] twinData <- twinData [!is.na(twinData $smk_cur1)&!is.na(twinData $smk_cur2),] head(twinData) nobs = dim(twinData)[1] twinData$age1 = scale(twinData$age1) twinData$age2 = scale(twinData$age2) # Select Ordinal Variables nth <- 1 # number of thresholds varso <- c("smk_cur","coronary") # list of ordinal variables names nvo <- 2 # number of ordinal variables ntvo <- nvo*2 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") ordData <- twinData # Select Variables for Analysis vars <- c("smk_cur","coronary") # list of variables names nv <- nvo # number of variables ntv <- nv*2 # number of total variables oc <- c(0,1,1) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis nc <- 2 # number of covariates def = c("sex", "age") ndef= length(def) defVars = paste(def,c(rep(1,ndef),rep(2,ndef)),sep="") # Select Data for Analysis mzData <- subset(twinData, zy_sex1==1|2, c(selVars,defVars)) dzData <- subset(twinData, zy_sex1==3|4, c(selVars,defVars)) mzData = mzData[complete.cases(mzData[,c(defVars, "smk_cur1", "smk_cur2")]),] dzData = dzData[complete.cases(dzData[,c(defVars, "smk_cur1", "smk_cur2")]),] mzDataF = mzData dzDataF = dzData mzDataF[,ordVars] <- mxFactor( x=mzDataF[,ordVars], levels=c(0,1) ) dzDataF[,ordVars] <- mxFactor( x=dzDataF[,ordVars], levels=c(0,1) ) # Raw data in OpenMx format dataMZ <- mxData(observed = mzDataF, type = "raw" ) dataDZ <- mxData(observed = dzDataF, type = "raw" ) # ---------------------Cholesky part!------------------------------------ ## Labeling aLabs <- c("aM","aC","aU") cLabs <- c("cM","cC","cU") eLabs <- c("eM","eC","eU") meanLabs <- c("meanM","meanP") aModLabs <- c("aMod11","aMod21","aMod22") cModLabs <- c("cMod11","cMod21","cMod22") eModLabs <- c("eMod11","eMod21","eMod22") threshLabs <- paste(varso,"thresh",sep="_") betaLabs_age <- paste("beta","Age",vars, sep="_") betaLabs_sex <- paste("beta","Sex",vars, sep="_") # Set Starting Values #frMV <- c(TRUE,FALSE,FALSE) # free status for variables frMV <- c(FALSE,FALSE) frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements frCv <- matrix(as.logical(frCvD),4) #svMe <- c(3,1.5,1.5) # start value for means svMe <- c(1.5,1.5) svPa <- .3 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .4 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix lbPa <- .0001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements svLTh <- 0.5 # start value for first threshold svITh <- 0.2 # start value for increments svTh <- -1 # start value for thresholds lbTh <- -1 # lower bounds for thresholds pathModVal = c(0,0.1,0.1) B_AgeVal = 0.5 B_SexVal = 0.5 labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_") # ------------------------------------------------------------------------------ # PREPARE MODEL # Matrices a, c, and e to store a, c, and e Path Coefficients pathA <- mxMatrix(name = "a", type = "Lower", nrow = nv, ncol = nv, free=T, labels = aLabs, values=svPaD, lbound=lbPaD) pathC <- mxMatrix(name = "c", type = "Lower", nrow = nv, ncol = nv, free=T, labels = cLabs, values=svPaD, lbound=lbPaD) pathE <- mxMatrix(name = "e", type = "Lower", nrow = nv, ncol = nv, free=T, labels = eLabs, values=svPeD, lbound=lbPaD) modPathA = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=aModLabs, name="aMod" ) modPathC = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=cModLabs, name="cMod" ) modPathE = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=eModLabs, name="eMod" ) # Moderator mod_tw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.smk_cur1', name="Mod1") mod_tw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.smk_cur2', name="Mod2") # Matrices generated to hold A, D, and E computed Variance Components varA1 <- mxAlgebra(name = "A1", expression = (a + Mod1%x%aMod) %*% t(a+ Mod1%x%aMod)) varC1 <- mxAlgebra(name = "C1", expression = (c + Mod1%x%cMod) %*% t(c+ Mod1%x%cMod)) varE1 <- mxAlgebra(name = "E1", expression = (e + Mod1%x%eMod) %*% t(e+ Mod1%x%eMod)) covA12 = mxAlgebra(name = "A12", expression = (a + Mod1%x%aMod) %*% t(a+ Mod2%x%aMod)) covC12 = mxAlgebra(name = "C12", expression = (c + Mod1%x%cMod) %*% t(c+ Mod2%x%cMod)) covE12 = mxAlgebra(name = "E12", expression = (e + Mod1%x%eMod) %*% t(e+ Mod2%x%eMod)) covA21 = mxAlgebra(name = "A21", expression = (a + Mod2%x%aMod) %*% t(a+ Mod1%x%aMod)) covC21 = mxAlgebra(name = "C21", expression = (c + Mod2%x%cMod) %*% t(c+ Mod1%x%cMod)) covE21 = mxAlgebra(name = "E21", expression = (e + Mod2%x%eMod) %*% t(e+ Mod1%x%eMod)) varA2 <- mxAlgebra(name = "A2", expression = (a + Mod2%x%aMod) %*% t(a+ Mod2%x%aMod)) varC2 <- mxAlgebra(name = "C2", expression = (c + Mod2%x%cMod) %*% t(c+ Mod2%x%cMod)) varE2 <- mxAlgebra(name = "E2", expression = (e + Mod2%x%eMod) %*% t(e+ Mod2%x%eMod)) myVarA <- mxAlgebra(name = "A0", expression = a %*% t(a)) myVarC <- mxAlgebra(name = "C0", expression = c %*% t(c)) #myVarE <- mxAlgebra(name = "E0", expression = e %*% t(e)) myVarE <- mxMatrix(type="Symm",nrow=nv, ncol=nv, free=c(F,T,F), values=c(1, 0.1, 1), labels=c("e_1_1", "e_2_1", "e_2_2"), lbound=c(0.0001,rep(NA,2)), name="E0") #mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(F,T,F), values=c(1, 0.1, 1), label=labLower("E",nv), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) per twin var1 <- mxAlgebra( A1+C1+E1, name="V1" ) var2 <- mxAlgebra( A2+C2+E2, name="V2" ) myVar <- mxAlgebra( A0+C0+E0, name="V0" ) # Constraint on variance of Binary variables #matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" ) #var_constraint <- mxConstraint( expression=diag2vec(V0)==Unv1, name="Var1" ) #matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" ) #var_constraint <- mxConstraint(expression=rbind(V0[2,2],V0[1,1]) == Unv1, name="Var1") ### Algebra for expected variance/covariance matrices expCovMZ <- mxAlgebra(name = "expCovMZ", expression = rbind (cbind(A1+C1+E1, A12+C12), cbind(A21+C21, A2+C2+E2))) expCovDZ <- mxAlgebra(name = "expCovDZ", expression = rbind (cbind(A1+C1+E1, 0.5%x%A12+C12), cbind(0.5%x%A21+C21, A2+C2+E2))) # Matrices for expected Means for females and males #setting up the regression intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=meanLabs, name="intercept" ) thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threshold <- mxAlgebra( expression= inc %*% thinG, name="Threshold" ) # Regression effects B_Age <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, labels=betaLabs_age, name="bAge" ) defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1","data.age2"), name="Age") B_Sex <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, labels=betaLabs_sex, name="bSex" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex") expMean <- mxAlgebra( intercept + (Age %x% bAge) + (Sex %x% bSex) , name="expMean") inclusions <- list (B_Age, B_Sex, defAge, defSex, intercept, expMean, threshold) # Objective objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="Threshold", threshnames=ordVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ,thresholds="Threshold", threshnames=ordVars ) funML <- mxFitFunctionML() # MZ and DZ models modelMZ <- mxModel(pathA, pathC, pathE, modPathA, modPathC, modPathE, mod_tw1, mod_tw2, varA1, varC1, varE1, covA12, covC12, covE12, covA21, covC21, covE21, varA2, varC2, varE2, var1, var2, myVarA, myVarC, myVarE, myVar, B_Age, B_Sex, defAge, defSex, intercept, expMean, thinG,inc,threshold, dataMZ, expCovMZ, expMZ, funML, name = "MZ") modelDZ <- mxModel(pathA, pathC, pathE, modPathA, modPathC, modPathE, mod_tw1, mod_tw2, varA1, varC1, varE1, covA12, covC12, covE12, covA21, covC21, covE21, varA2, varC2, varE2, var1, var2, myVarA, myVarC, myVarE, myVar, B_Age, B_Sex, defAge, defSex, intercept, expMean,thinG,inc, threshold, dataDZ, expCovDZ, expDZ, funML, name = "DZ") multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) ACEmodModel <- mxModel( "ACEmod", modelMZ, modelDZ, funML, multi ,mxCI(c("aMod21","aMod22","cMod21","cMod22","eMod21","eMod22"))) ACEmodFit <- mxRun(ACEmodModel,intervals=TRUE) summary(ACEmodFit) set.seed(1) TACEmodFit <- mxTryHardOrdinal(ACEmodModel,extraTries = 30,intervals=TRUE, bestInitsOutput=TRUE , silent=FALSE,OKstatuscodes=c(0,1,5)) summary(TACEmodFit) ##=======================================================================# # 2. main effect MainEffectsModel = mxModel (TACEmodFit, name='MainEffects') MainEffectsModel = omxSetParameters(MainEffectsModel, labels=c('aMod21','aMod22'), values=0, free=F) MainEffectsModel = omxSetParameters(MainEffectsModel, labels=c('cMod21','cMod22'), values=0, free=F) MainEffectsModel = omxSetParameters(MainEffectsModel, labels=c('eMod21','eMod22'), values=0, free=F) set.seed(1) TMainEffectsFit <- mxTryHardOrdinal(MainEffectsModel,extraTries = 30,intervals=TRUE, bestInitsOutput=TRUE) summary(TMainEffectsFit) mxCompare(TACEmodFit,TMainEffectsFit) ##=======================================================================# # 3. drop ac NacEffectsModel = mxModel (TACEmodFit, name='NacEffects') NacEffectsModel = omxSetParameters(NacEffectsModel, labels=c('aMod21'), values=0, free=F) set.seed(1) TNacEffectsFit <- mxTryHardOrdinal(NacEffectsModel,extraTries = 30,intervals=TRUE, bestInitsOutput=TRUE) summary(TNacEffectsFit) mxCompare(TACEmodFit,TNacEffectsFit) ##=======================================================================# # 4. drop cc NccEffectsModel = mxModel (TACEmodFit, name='NccEffects') NccEffectsModel = omxSetParameters(NccEffectsModel, labels=c('cMod21'), values=0, free=F) set.seed(1) TNccEffectsFit <- mxTryHardOrdinal(NccEffectsModel,extraTries = 30,intervals=TRUE) summary(TNccEffectsFit) mxCompare(TACEmodFit,TNccEffectsFit) ##=======================================================================# # 5. drop ec NecEffectsModel = mxModel (TACEmodFit, name='NecEffects') NecEffectsModel = omxSetParameters(NecEffectsModel, labels=c('eMod21'), values=0, free=F) set.seed(1) TNecEffectsFit <- mxTryHardOrdinal(NecEffectsModel,extraTries = 30,intervals=TRUE) summary(TNecEffectsFit) mxCompare(TACEmodFit,TNecEffectsFit) ##=======================================================================# # 6 drop au NauEffectsModel = mxModel (TACEmodFit, name='NauEffects') NauEffectsModel = omxSetParameters(NauEffectsModel, labels=c('aMod22'), values=0, free=F) set.seed(1) TNauEffectsFit <- mxTryHardOrdinal(NauEffectsModel,extraTries = 30,intervals=TRUE) summary(TNauEffectsFit) mxCompare(TACEmodFit,TNauEffectsFit) ##=======================================================================# # 7. drop cu NcuEffectsModel = mxModel (TACEmodFit, name='NcuEffects') NcuEffectsModel = omxSetParameters(NcuEffectsModel, labels=c('cMod22'), values=0, free=F) set.seed(1) TNcuEffectsFit <- mxTryHardOrdinal(NcuEffectsModel,extraTries = 30,intervals=TRUE) summary(TNcuEffectsFit) mxCompare(TACEmodFit,TNcuEffectsFit) ##=======================================================================# # 8. drop eu NeuEffectsModel = mxModel (TACEmodFit, name='NeuEffects') NeuEffectsModel = omxSetParameters(NeuEffectsModel, labels=c('eMod22'), values=0, free=F) set.seed(1) TNeuEffectsFit <- mxTryHardOrdinal(NeuEffectsModel,extraTries = 30,intervals=TRUE) summary(TNeuEffectsFit) mxCompare(TACEmodFit,TNeuEffectsFit) sink()