rm(list=ls()) setwd("L:/课题/王冬萌小论文复核/ACE") require(OpenMx) require(psych) library(foreign) require(readstata13) source('GenEpiHelperFunctions.R') mxVersion() #=======================================================================# # PREPARE DATA # #=======================================================================# # Load Data #sink("25+ACE-drink_cur.txt") mxOption(NULL, "Default optimizer", "NPSOL") #mxOption( NULL, 'Default optimizer' ,'SLSQP' ) # mxOption(NULL, "Default optimizer", "CSOLNP") # ----------------------------------------------------------------------- # PREPARE DATA # Read Data #twinData <- read.table("D:/FHI/Sosiale forhold og helse/Analyses/IBS/SFH_17_09_2018_subset_pss_discordance_strain_lowsupport_family_format.txt", header=T, sep='\t', dec=',') twinData <-read.dta13("登记-25岁成年双胞胎-34581pair-coronary-ssex.dta") twinData <- twinData [!is.na(twinData $coronary1)&!is.na(twinData $coronary2),] twinData <- twinData [!is.na(twinData $drink_cur1)&!is.na(twinData $drink_cur2),] str(twinData) nobs = dim(twinData)[1] twinData$age1 = scale(twinData$age1) twinData$age2 = scale(twinData$age2) #twinData$Bmi1 = scale(twinData$bmi1) #twinData$Bmi2 = scale(twinData$bmi2) # Select Ordinal Variables nth <- 1 # number of thresholds varso <- c('coronary') # 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 Vars <- c("drink_cur","coronary") 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 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, "drink_cur1", "drink_cur2")]),] dzData = dzData[complete.cases(dzData[,c(defVars, "drink_cur1", "drink_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!------------------------------------ # Set up Cholesky ADE decomposition, with RawData and Matrices Input # ----------------------------------------------------------------------- ## 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) # free status for variables 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 <- .4 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .8 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix lbPa <- 0 # 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 svTh <- 1.5 # start value for thresholds pathModVal = c(0,0.1,0.1) B_AgeVal = 0.5 B_SexVal = 0.5 ## Modeling # 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" ) modPathA = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,F,F), values=c(0,0,0), labels=aModLabs, name="aMod" ) modPathC = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,F,F), values=c(0,0,0), labels=cModLabs, name="cMod" ) modPathE = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,F,F), values=c(0,0,0), labels=eModLabs, name="eMod" ) #MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('aMod21','aMod22'), values=c(0,0), free=FALSE) #MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('cMod21','cMod22'), values=c(0,0), free=FALSE) #MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('eMod21','eMod22'), values=c(0,0), free=FALSE) # Moderator mod_tw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.drink_cur1', name="Mod1") mod_tw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.drink_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(V0[2,2] == 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" ) threshold <-mxMatrix( type="Full", nrow=1, ncol=ntvo, free=T, values=svTh, labels=threshLabs, 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.drink_cur1","data.drink_cur2"), name="Age") 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") 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, matUnv, var_constraint, B_Age, B_Sex, defAge, defSex, intercept, expMean, 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, var1, var2, myVarA, myVarC, myVarE, myVar, matUnv, B_Age, B_Sex, defAge, defSex, intercept, expMean, threshold, dataDZ, expCovDZ, expDZ, funML, name = "DZ") #plan <- omxDefaultComputePlan() #plan$steps <- list( # NM=mxComputeNelderMead(centerIniSimplex=TRUE), # GD=plan$steps$GD,ND=plan$steps$ND,SE=plan$steps$SE,RD=plan$steps$RD,RE=plan$steps$RD) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) #ACEmodModel <- mxModel( "ACEmod", modelMZ, modelDZ, funML, multi MainEffectsModel<- mxModel( "MainEffects", modelMZ, modelDZ, funML, multi ) #ACEmodFit <- mxAutoStart(ACEmodModel) #ACEmodFit <- mxRun(ACEmodModel) #ACEmodFit <- mxTryHardOrdinal(ACEmodModel) #summary(ACEmodFit) MainEffectsFit <- mxRun(MainEffectsModel) set.seed(1) MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel,extraTries = 30) summary(MainEffectsFit) #model <- mxTryHard(ACEmodFit) mxGetExpected(MainEffectsModel, "covariance", "MZ") #MainEffectsModel <- mxModel (ACEmodFit, name='MainEffects') #MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('aMod21','aMod22'), values=c(0,0), free=FALSE) #MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('cMod21','cMod22'), values=c(0,0), free=FALSE) #MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('eMod21','eMod22'), values=c(0,0), free=FALSE) #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" ) ACEmodModel <- mxModel (MainEffectsModel, name='MainEffects') ACEmodModel <- omxSetParameters(ACEmodModel, labels=c('aMod21','aMod22'), values=c(pathModVal,pathModVal), free=c(F,T,T)) ACEmodModel <- omxSetParameters(ACEmodModel, labels=c('cMod21','cMod22'), values=c(pathModVal,pathModVal), free=c(F,T,T)) ACEmodModel <- omxSetParameters(ACEmodModel, labels=c('eMod21','eMod22'), values=c(pathModVal,pathModVal), free=c(F,T,T)) ACEmodModel = omxSetParameters(ACEmodModel, labels=c('aMod21','aMod22'), values=c(0.1,0.1), free=c(T,T)) ACEmodModel = omxSetParameters(ACEmodModel, labels=c('cMod21','cMod22'), values=c(0.1,0.1), free=c(T,T)) ACEmodModel = omxSetParameters(ACEmodModel, labels=c('eMod21','eMod22'), values=c(0.1,0.1), free=c(T,T)) #ACEmodFit <- mxAutoStart(MainEffectsModel) #ACEmodFit <- mxRun(MainEffectsModel) ACEmodFit <- mxRun(ACEmodModel) set.seed(1) ACEmodFit <- mxTryHardOrdinal(ACEmodModel,extraTries = 80) mxCompare(ACEmodFit,MainEffectsFit) summary(ACEmodFit) #=======================================================================# # PLOT MODERATION FULL MODEL # #=======================================================================# drink_cur <- sort(unique(scale(rbind(mzData$drink_cur1, mzData$drink_cur2, dzData$drink_cur1, dzData$drink_cur2)))) drink_cur_length=length(drink_cur) # Purcell's representation (total variances) Va=array(NA, dim=c(nv,nv,drink_cur_length)) Vc=array(NA, dim=c(nv,nv,drink_cur_length)) Ve=array(NA, dim=c(nv,nv,drink_cur_length)) Vt=array(NA, dim=c(nv,nv,drink_cur_length)) SVa=array(NA, dim=c(nv,nv,drink_cur_length)) SVc=array(NA, dim=c(nv,nv,drink_cur_length)) SVe=array(NA, dim=c(nv,nv,drink_cur_length)) corA = array(NA, dim=c(nv,nv,drink_cur_length)) corC = array(NA, dim=c(nv,nv,drink_cur_length)) corE = array(NA, dim=c(nv,nv,drink_cur_length)) corP = array(NA, dim=c(nv,nv,drink_cur_length)) for (i in 1:drink_cu_length) { Va[,,i] <- mxEval((a + aMod%x%drink_cur[i]) %*% t(a + aMod%x%drink_cur[i]), ACEmodFit$MZ) Vc[,,i] <- mxEval((c + cMod%x%drink_cur[i]) %*% t(c + cMod%x%drink_cur[i]), ACEmodFit$MZ) Ve[,,i] <- mxEval((e + eMod%x%drink_cur[i]) %*% t(e + eMod%x%drink_cur[i]), ACEmodFit$MZ) Vt[,,i] <- Va[,,i] + Vc[,,i] + Ve[,,i] SVa[,,i] <- Va[,,i]/Vt[,,i] SVc[,,i] <- Vc[,,i]/Vt[,,i] SVe[,,i] <- Ve[,,i]/Vt[,,i] corA[,,i] <- solve(sqrt(diag(2)*Va[,,i]))%*%Va[,,i]%*%solve(sqrt(diag(2)*Va[,,i])) corC[,,i] <- solve(sqrt(diag(2)*Vc[,,i]))%*%Vc[,,i]%*%solve(sqrt(diag(2)*Vc[,,i])) corE[,,i] <- solve(sqrt(diag(2)*Ve[,,i]))%*%Ve[,,i]%*%solve(sqrt(diag(2)*Ve[,,i])) corP[,,i] <- solve(sqrt(diag(2)*Vt[,,i]))%*%Vt[,,i]%*%solve(sqrt(diag(2)*Vt[,,i])) } out <- as.matrix(cbind(Va[2,2,],Vc[2,2,],Ve[2,2,],Vt[2,2,],SVa[2,2,],SVc[2,2,],SVe[2,2,], corA[2,1,],corC[2,1,],corE[2,1,], corP[2,1,])) names(out) = c('Va', 'Vc', 'Ve', 'Vt', 'SVa', 'SVc', 'SVe', 'corA', 'corC', 'corE', 'corP') head(out) par(mar=c(5.1,4.1,4.1,2.1)) par(mfrow=c(1,2)) matplot(drink_cur, out[,1:4], type="l", lty=4:1, col=4:1, xlab="drink_cur", ylab="Variance Component", main="Total - Unstandardized\nvariance components") legend(0,1, c("Va","Vc","Ve","Vt"), lty=4:1, col=4:1, cex = 1, y.intersp = .5) abline(v = 0, lty = 2, lwd = 0.5) matplot(drink_cur, out[,5:7], type="l", lty=4:2, col=4:2, xlab="drink_cur", ylab="Variance Component", main="Standardized\nvariance components") legend(0,0.4, c( "SVa", "SVc", "SVe"), lty=4:1, col=4:1, cex = 1, y.intersp = .2) abline(v = 0, lty = 2, lwd = 0.5) par(mfrow=c(1,1)) matplot(drink_cur, out[,8:11], type='l', lty=4:1, col=4:1, xlab='drink_cur', ylab="Correlation", main="Correlation") legend(4,0.4, c("rA","rC","rE","rP"), lty=4:1, col=4:1, cex = 1, y.intersp = .5) abline(v = 0, lty = 2, lwd = 0.5) par(mfrow=c(1,1)) stackpoly(drink_cur, out[,1:3], stack = T, col = c('black', 'gray', 'gray40'), axis4 = F, xat = -1:6, ylim = c(0, 13), xlab='Bmi') abline(v = 0, lty = 2, lwd = 0.5) legend(3,12, fill = c('black', 'gray', 'gray40'), legend = c("A","C","E"), border = 'black', bty = 'n') sink()