rm(list=ls()) setwd("E:/thesis data") require(OpenMx) library(foreign) library(openxlsx) #=======================================================================# # PREPARE DATA # #=======================================================================# mxVersion() #mxOption(NULL, "Default optimizer", "NPSOL") mxOption( NULL, 'Default optimizer' ,'SLSQP' ) # mxOption(NULL, "Default optimizer", "CSOLNP") #mxOption( NULL, "mvnRelEps", value= mxOption(NULL, key="mvnRelEps")/5) # ----------------------------------------------------------------------- # Read Data twinData <-read.xlsx("same_sex_rear_togeter8_5.xlsx") # Select Ordinal Variables nth <- 1 # number of thresholds varso <- c('T2DM') # list of ordinal variables names nvo <- length(varso) # number of ordinal variables ntvo <- nvo*2 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") # Select Variables for Analysis mod="smk_cur" Vars <- c(mod,"T2DM") nv <- length(Vars) # number of variables nind <- 2 ntv <- nv*nind # number of total variables selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") def = c("sex", "age") ndef= length(def) defVars = paste(def,c(rep(1,ndef),rep(2,ndef)),sep="") # Select Data for Analysis twinData$age1 = scale(twinData$age1) twinData$age2 = scale(twinData$age2) mzData <- subset(twinData, zy_sex==1|zy_sex==2, c(selVars,defVars)) dzData <- subset(twinData, zy_sex==3|zy_sex==4, c(selVars,defVars)) dim(mzData) dim(dzData) mzData = mzData[complete.cases(mzData[,c(selVars,defVars)]),] dzData = dzData[complete.cases(dzData[,c(selVars,defVars)]),] mzDataF = mzData dzDataF = dzData mzDataF[,ordVars] <- mxFactor( x=mzDataF[,ordVars], levels=c(0,1) ) dzDataF[,ordVars] <- mxFactor( x=dzDataF[,ordVars], levels=c(0,1) ) dim(mzDataF) dim(dzDataF) # Raw data in OpenMx format dataMZ <- mxData(observed = mzDataF, type = "raw" ) dataDZ <- mxData(observed = dzDataF, type = "raw" ) nv <- 1 # number of phenotype variables per twin ntv <- nv*2 # total number of phenotype variables na <- 1 # number of A factors (per twin) nc <- 1 # number of C factors (per twin) ne <- 1 # number of E factors (per twin) # Set Starting Values frMV <- c(FALSE,FALSE) svTh <- 0 # start value for thresholds lbTh <- -5 # lower bounds for thresholds svMe <- c(0,0) # start value for means # Create Labels labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_") modVars <- paste(mod,c(1,2),sep='') selVars <- paste(varso,c(1,2),sep='') betaLabs_age <- paste("beta","Age",Vars, sep="_") betaLabs_sex <- paste("beta","Sex",Vars, sep="_") # Matrices to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Full", nrow=nv, ncol=na, free=TRUE, values=.5, lbound = 0.0001, label="a11", name="a" ) pathC <- mxMatrix( type="Full", nrow=nv, ncol=nc, free=TRUE, values=.5, lbound = 0.0001, label="c11", name="c" ) pathE <- mxMatrix( type="Full", nrow=nv, ncol=ne, free=TRUE, values=.5, lbound = 0.0001, label="e11", name="e" ) # Matrices to store the moderating a, c, and e Path Coefficients modPathA <- mxMatrix( type='Full', nrow=nv, ncol=na, free=TRUE, values=.1, label="aM11", name="aM" ) modPathC <- mxMatrix( type='Full', nrow=nv, ncol=nc, free=TRUE, values=.1, label="cM11", name="cM" ) modPathE <- mxMatrix( type='Full', nrow=nv, ncol=ne, free=TRUE, values=.1, label="eM11", name="eM" ) # Matrices for the expected means vector intercept <- mxMatrix( type="Zero", nrow=1, ncol=1, name="int" ) PathM <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="m11", name="m" ) threshold <-mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="Threshold" ) B_Age <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, labels=betaLabs_age, name="bAge" ) defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1","data.age2"), name="Age") B_Sex <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, labels=betaLabs_sex, name="bSex" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex") # Matrix for the moderator variable mod1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=paste("data.",modVars[1],sep=""), name="Mod1") mod2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=paste("data.",modVars[2],sep=""), name="Mod2") # Matrices to compute the moderated A, C, and E variance components covAmod1 <- mxAlgebra( expression=(a+ Mod1%*%aM) %*% t(a+ Mod1%*%aM), name="A1" ) covCmod1 <- mxAlgebra( expression=(c+ Mod1%*%cM) %*% t(c+ Mod1%*%cM), name="C1" ) covEmod1 <- mxAlgebra( expression=(e+ Mod1%*%eM) %*% t(e+ Mod1%*%eM), name="E1" ) covAmod2 <- mxAlgebra( expression=(a+ Mod2%*%aM) %*% t(a+ Mod2%*%aM), name="A2" ) covCmod2 <- mxAlgebra( expression=(c+ Mod2%*%cM) %*% t(c+ Mod2%*%cM), name="C2" ) covEmod2 <- mxAlgebra( expression=(e+ Mod2%*%eM) %*% t(e+ Mod2%*%eM), name="E2" ) covAmod12 <- mxAlgebra( expression=(a+ Mod1%*%aM) %*% t(a+ Mod2%*%aM), name="A12" ) covCmod12 <- mxAlgebra( expression=(c+ Mod1%*%cM) %*% t(c+ Mod2%*%cM), name="C12" ) covAmod21 <- mxAlgebra( expression=(a+ Mod2%*%aM) %*% t(a+ Mod1%*%aM), name="A21" ) covCmod21 <- mxAlgebra( expression=(c+ Mod2%*%cM) %*% t(c+ Mod1%*%cM), name="C21" ) # Algebra to compute the total variance per twin covPmod <- mxAlgebra( expression=A1+C1+E1, name="V" ) # Matrices to compute the non-moderated A, C, and E Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="nA" ) covC <- mxAlgebra( expression=c %*% t(c), name="nC" ) covE <- mxAlgebra( expression=e %*% t(e), name="nE" ) # Algebra to compute the total variance per twin covP <- mxAlgebra( expression = nA+nC+nE, name="nV" ) # Constrain Variance of Binary Variables matUnv <- mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv1" ) var1 <- mxConstraint( expression=diag2vec( nV )== Unv1, name="Var1" ) # Algebra for the expected mean vector and expected variance/covariance matrices covMZ <- mxAlgebra( expression= rbind (cbind(A1+C1+E1, A12+C12), cbind(A21+C21, A2+C2+E2)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind (cbind(A1+C1+E1, 0.5%x%A12+C12), cbind(0.5%x%A21+C21, A2+C2+E2)), name="expCovDZ" ) wmod1 <- mxAlgebra( expression= m %*% Mod1, name="DR1" ) wmod2 <- mxAlgebra( expression= m %*% Mod2, name="DR2" ) meanG <- mxAlgebra( expression= cbind((int + DR1+(Age %x% bAge) + (Sex %x% bSex)),(int + DR2+(Age %x% bAge) + (Sex %x% bSex))), name="expMeanG" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanG", dimnames=selVars, thresholds="Threshold", threshnames=ordVars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanG", dimnames=selVars, thresholds="Threshold", threshnames=ordVars ) fitMZ <- mxFitFunctionML() fitDZ <- mxFitFunctionML() # Combine Groups parsZ <- list( pathA, pathC, pathE, modPathA, modPathC, modPathE, PathM,B_Age, B_Sex, intercept,threshold) modelMZ <- mxModel( parsZ, meanG, covAmod1, covCmod1, covEmod1,covAmod2,covCmod2,covEmod2, covAmod12,covCmod12,covAmod21,covCmod21, covPmod, covA, covC, covE, covP, matUnv, var1, mod1, mod2, defAge, defSex, wmod1,wmod2, covMZ, dataMZ, objMZ, fitMZ, name="MZ" ) modelDZ <- mxModel( parsZ, meanG, covAmod1, covCmod1, covEmod1,covAmod2,covCmod2,covEmod2, covAmod12,covCmod12,covAmod21,covCmod21, covPmod, mod1,mod2, defAge, defSex, wmod1,wmod2, covDZ, dataDZ, objDZ, fitDZ, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective+ DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) CI <- mxCI(c("a11","c11","e11","aM11","cM11","eM11","m11")) ACEModel <- mxModel( "ModACE", parsZ, modelMZ, modelDZ, minus2ll, obj,CI ) # RUN MODEL #set.seed(1) #ACEmodFit <- mxTryHardOrdinal(ACEModel,extraTries = 10,intervals=TRUE) ACEmodFit <- mxRun( ACEModel, intervals=T ) ACEsum <- summary( ACEmodFit ) round(ACEsum$CI[, c(2, 1,3)],2) ##=======================================================================# # 2. main effect MainEffectsModel = mxModel (ACEmodFit, name='MainEffects') MainEffectsModel = omxSetParameters(MainEffectsModel, labels=c('aM11','cM11','eM11'), values=0, free=F) MainEffectsFit <- mxRun(MainEffectsModel,intervals=TRUE) Mainsum <- summary(MainEffectsFit) round(Mainsum$CI[,c(2,1,3)],2) mxCompare(ACEmodFit,MainEffectsFit) ##=======================================================================# # 3. drop ac NaMEffectsModel = mxModel (ACEmodFit, name='NaMEffects') NaMEffectsModel = omxSetParameters(NaMEffectsModel, labels=c('aM11'), values=0, free=F) NaMEffectsFit <- mxRun(NaMEffectsModel,intervals=TRUE) NaMsum <- summary(NaMEffectsFit) mxCompare(ACEmodFit,NaMEffectsFit) ##=======================================================================# # 4. drop cc NcMEffectsModel = mxModel (ACEmodFit, name='NcMEffects') NcMEffectsModel = omxSetParameters(NcMEffectsModel, labels=c('cM11'), values=0, free=F) NcMEffectsFit <- mxRun(NcMEffectsModel,intervals=TRUE) NcMsum <- summary(NcMEffectsFit) mxCompare(ACEmodFit,NcMEffectsFit) ##=======================================================================# # 5. drop ec NeMEffectsModel = mxModel (ACEmodFit, name='NeMEffects') NeMEffectsModel = omxSetParameters(NeMEffectsModel, labels=c('eM11'), values=0, free=F) NeMEffectsFit <- mxRun(NeMEffectsModel,intervals=TRUE) NeMsum <- summary(NeMEffectsFit) mxCompare(ACEmodFit,NeMEffectsFit) #6. only ac OaMEffectsModel <- omxSetParameters(model=ACEmodFit,labels = c('cM11','eM11'),free=c(F,F),values = c(0,0), name="OaMEffects") OaMEffectFit <- mxRun(OaMEffectsModel,intervals = TRUE) OaMsum <- summary(OaMEffectFit,verbose=T) round(OaMsum$CI[,c(2,1,3)],2) mxCompare(ACEmodFit, OaMEffectFit) #7. only cc OcMEffectsModel <- omxSetParameters(model=ACEmodFit,labels = c('aM11','eM11'),free=c(F,F),values = c(0,0), name="OcMEffects") OcMEffectFit <- mxRun(OcMEffectsModel,intervals = TRUE) OcMsum <- summary(OcMEffectFit,verbose=T) round(OcMsum$CI[,c(2,1,3)],2) mxCompare(ACEmodFit, OcMEffectFit) #8. only ec OeMEffectsModel <- omxSetParameters(model=ACEmodFit,labels = c('aM11','cM11'),free=c(F,F),values = c(0,0), name="OeMEffects") OeMEffectFit <- mxRun(OeMEffectsModel,intervals = TRUE) OeMsum <- summary(OeMEffectFit,verbose=T) round(OeMsum$CI[,c(2,1,3)],2) mxCompare(ACEmodFit, OeMEffectFit) com <- mxCompare(ACEmodFit,c(NaMEffectsFit,NcMEffectsFit,NeMEffectsFit,OaMEffectFit,OcMEffectFit,OeMEffectFit,MainEffectsFit)) com <- com[,-c(1,3)] com[,c(2,4,5)] <- round(com[,c(2,4,5)],2) com[,7] <- ifelse(com[,7]<0.001, "<0.001",round(com[,7],3)) ace <- rbind(round(ACEsum$CI[,c(2,1,3)],2), round(NaMsum$CI[,c(2,1,3)],2), round(NcMsum$CI[,c(2,1,3)],2), round(NeMsum$CI[,c(2,1,3)],2), round(OaMsum$CI[,c(2,1,3)],2), round(OcMsum$CI[,c(2,1,3)],2), round(OeMsum$CI[,c(2,1,3)],2), round(Mainsum$CI[,c(2,1,3)],2)) for(i in c(1:length(ace[,1]))){ ace[i,1] <- paste(ace[i,1],"(",ace[i,2],",",ace[i,3],")",sep="") } uni_smoke <- data.frame(model=c("ACE_XYZ","ACE_YZ","ACE_XZ","ACE_XY","ACE_X","ACE_Y","ACE_Z","ACE"), a11=c(ace[1,1],ace[8,1],ace[15,1],ace[22,1],ace[29,1],ace[36,1],ace[43,1],ace[50,1]), c11=c(ace[2,1],ace[9,1],ace[16,1],ace[23,1],ace[30,1],ace[37,1],ace[44,1],ace[51,1]), e11=c(ace[3,1],ace[10,1],ace[17,1],ace[24,1],ace[31,1],ace[38,1],ace[45,1],ace[52,1]), am11=c(ace[4,1],ace[11,1],ace[18,1],ace[25,1],ace[32,1],ace[39,1],ace[46,1],ace[53,1]), cm11=c(ace[5,1],ace[12,1],ace[19,1],ace[26,1],ace[33,1],ace[40,1],ace[47,1],ace[54,1]), em11=c(ace[6,1],ace[13,1],ace[20,1],ace[27,1],ace[34,1],ace[41,1],ace[48,1],ace[55,1]), m11=c(ace[7,1],ace[14,1],ace[21,1],ace[28,1],ace[35,1],ace[42,1],ace[49,1],ace[56,1])) uni_smoke <- cbind(uni_smoke,com) sink('E:/thesis data/result/moderation/univariate/uni_smoke_8.19.doc',append=T) uni_smoke sink()