rm(list=ls()) setwd("L:/课题/王冬萌小论文复核/ACE") require(OpenMx) require(psych) library(foreign) require(readstata13) require(ggplot) source('GenEpiHelperFunctions.R') #=======================================================================# # PREPARE DATA # #=======================================================================# #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("登记-25岁成年双胞胎-20612pair-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 Continuous Variables #varsc <- c('SL') # list of continuous variables names #nvc <- 1 # number of continuous variables #ntvc <- nvc*2 # number of total continuous variables #conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") nvc <- 0 # number of continuous variables ntvc <- nvc*2 # number of total continuous variables # 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 <- nvc+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 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$smk_cur1_num <- as.numeric(mzDataF$smk_cur1) mzDataF$smk_cur2_num <- as.numeric(mzDataF$smk_cur2) dzDataF$smk_cur1_num <- as.numeric(mzDataF$smk_cur1) dzDataF$smk_cur2_num <- as.numeric(mzDataF$smk_cur2) mzDataF[,ordVars] <- mxFactor( x=mzDataF[,ordVars], levels=c(0,1) ) dzDataF[,ordVars] <- mxFactor( x=dzDataF[,ordVars], levels=c(0,1) ) mzDataF$coronary1[is.na(mzDataF$smk_cur1_num)] <- NA mzDataF$smk_cur1_num[is.na(mzDataF$smk_cur1_num)] <- -999.0 mzDataF$coronary2[is.na(mzDataF$smk_cur2_num)] <- NA mzDataF$smk_cur2_num[is.na(mzDataF$smk_cur2_num)] <- -999.0 dzDataF$coronary1[is.na(mzDataF$smk_cur1_num)] <- NA dzDataF$smk_cur1_num[is.na(mzDataF$smk_cur1_num)] <- -999.0 dzDataF$coronary2[is.na(mzDataF$smk_cur2_num)] <- NA dzDataF$smk_cur2_num[is.na(mzDataF$smk_cur2_num)] <- -999.0 # 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="_") betaLabs_smk <- paste("beta","Smk",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(0,0) # start value for means 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 # Create Labels #labMe <- paste("mean",vars,sep="_") labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_") #labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # 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_num', name="Mod1") mod_tw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.smk_cur2_num', name="Mod2") #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)) # 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" ) ### 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=.1, labels=betaLabs_age, name="bAge" ) defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1","data.age2"), name="Age") # Create Matrices for Covariates and linear Regression Coefficients B_Sex <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_sex, name="bSex" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex") B_Smk <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_smk, name="bSmk" ) defSmk <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.smk_cur1","data.smk_cur2"), name="Smk") expMean <- mxAlgebra( intercept + (Age %x% bAge) + (Sex %x% bSex) + (Smk %x% bSmk) , name="expMean") inclusions <- list (B_Age, B_Sex, B_Smk,defAge, defSex, defSmk,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, B_Age, B_Sex, B_Smk, defAge, defSex,defSmk, 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, matUnv, var_constraint, # B_Age, B_Sex, defAge, defSex, intercept, expMean, threshold, # dataDZ, expCovDZ, expDZ, funML, name = "DZ") modelDZ <- mxModel(pathA, pathC, pathE, modPathA, modPathC, modPathE, mod_tw1, mod_tw2, varA1, varC1, varE1, covA12, covC12, covE12, covA21, covC21, covE21, varA2, varC2, varE2, myVarA, myVarC, myVarE, myVar, matUnv, var_constraint, B_Age, B_Sex, B_Smk, defAge, defSex,defSmk, 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)