#trivariate #full ACE #ACE with c21, c22, c32 = 0 (males); C for intellect and its covariances = 0 setwd('/Users/iza/iza/project') # Load libraries source('http://openmx.psyc.virginia.edu/getOpenMx.R') source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R") require(psych) require(OpenMx) # PREPARE DATA Data <- read.table (file.choose (), header=T, sep="\t",na=c(-1, -99)) describe(Data, skew=F) # Select Variables for Analysis Vars <- c('CITOsc_','intellect_','open_') nv <- 3 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Data for Analysis mzmData <- subset(Data, zyg5==1, selVars) dzmData <- subset(Data, zyg5==2, selVars) mzfData <- subset(Data, zyg5==3, selVars) dzfData <- subset(Data, zyg5==4, selVars) dosmfData <- subset(Data, zyg5==5, selVars) # ---------------------Cholesky part!------------------------------------ # Set up Cholesky ACE decomposition, with RawData and Matrices Input # ----------------------------------------------------------------------- ## Labeling aLabsm <- c("a11m","a21m","a22m","a31m","a32m","a33m") cLabsm <- c("c11m","c21m","c22m","c31m","c32m","c33m") eLabsm <- c("e11m","e21m","e22m","e31m","e32m","e33m") aLabsf <- c("a11f","a21f","a22f","a31f","a32f","a33f") cLabsf <- c("c11f","c21f","c22f","c31f","c32f","c33f") eLabsf <- c("e11f","e21f","e22f","e31f","e32f","e33f") ## Modeling # Matrices a, c, and e to store a, c, and e Path Coefficients pathAm <- mxMatrix(name = "am", type = "Lower", nrow = nv, ncol = nv, free = TRUE, values = c(1, .8, 1, .8, .8, 1),labels = c("a11m","a21m","a22m","a31m","a32m","a33m")) pathCm <- mxMatrix(name = "cm", type = "Lower", nrow = nv, ncol = nv, free = TRUE, values = c(1, .8, 1, .8, .8, 1),labels = c("c11m","c21m","c22m","c31m","c32m","c33m")) pathEm <- mxMatrix(name = "em", type = "Lower", nrow = nv, ncol = nv, free = TRUE, values = c(1, .8, 1, .8, .8, 1),labels = c("e11m","e21m","e22m","e31m","e32m","e33m")) pathAf <- mxMatrix(name = "af", type = "Lower", nrow = nv, ncol = nv, free = TRUE, values = c(1, .8, 1, .8, .8, 1),labels = c("a11f","a21f","a22f","a31f","a32f","a33f")) pathCf <- mxMatrix(name = "cf", type = "Lower", nrow = nv, ncol = nv, free = TRUE, values = c(1, .8, 1, .8, .8, 1),labels = c("c11f","c21f","c22f","c31f","c32f","c33f")) pathEf <- mxMatrix(name = "ef", type = "Lower", nrow = nv, ncol = nv, free = TRUE, values = c(1, .8, 1, .8, .8, 1),labels = c("e11f","e21f","e22f","e31f","e32f","e33f")) # Matrices generated to hold A, C, and E computed Variance Components covAm <- mxAlgebra(expression = am %*% t(am), name = "Am") covAf <- mxAlgebra(expression = af %*% t(af), name = "Af") covCm <- mxAlgebra(expression = cm %*% t(cm), name = "Cm") covCf <- mxAlgebra(expression = cf %*% t(cf), name = "Cf") covEm <- mxAlgebra(expression = em %*% t(em), name = "Em") covEf <- mxAlgebra(expression = ef %*% t(ef), name = "Ef") # Algebra to compute total variances and standard deviations (diagonal only) covVm <- mxAlgebra(expression = Am+Cm+Em, name = "Vm") covVf <- mxAlgebra(expression = Af+Cf+Ef, name = "Vf") # Algebra to compute total variances and standard deviations (diagonal only) matI <- mxMatrix(name= "I", type="Iden", nrow = nv, ncol = nv) iSDm <- mxAlgebra(name ="iSDm", expression = solve(sqrt(I*Vm))) iSDf <- mxAlgebra(name ="iSDf", expression = solve(sqrt(I*Vf))) corPhm <- mxAlgebra(name ="rPhm", expression = iSDm%*%Vm%*%iSDm) corPhf <- mxAlgebra(name ="rPhf", expression = iSDf%*%Vf%*%iSDf) # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins # Mean structure, Algebra M to store Expected means grandMeanm <- mxMatrix(type="Full", nrow=1, ncol=nv, free=TRUE, labels=c("mean1m","mean2m", "mean3m"), values=c(539, 16, 19), name="Mm") grandMeanf <- mxMatrix(type="Full", nrow=1, ncol=nv, free=TRUE, labels=c("mean1f","mean2f","mean3f"), values=c(538, 16, 19), name="Mf") expMeanMZM <- mxAlgebra(expression= cbind(Mm,Mm) , name="expMeanMZM") expMeanDZM <- mxAlgebra(expression= cbind(Mm,Mm) , name="expMeanDZM") expMeanMZF <- mxAlgebra(expression= cbind(Mf,Mf) , name="expMeanMZF") expMeanDZF <- mxAlgebra(expression= cbind(Mf,Mf) , name="expMeanDZF") expMeanDOSmf <- mxAlgebra(expression= cbind(Mm,Mf) , name="expMeanDOSmf") expCovMZM <- mxAlgebra( expression= rbind( cbind(Vm, Am+Cm), cbind(Am+Cm, Vm)), name="expCovMZM" ) expCovDZM <- mxAlgebra( expression= rbind( cbind(Vm, 0.5%x%Am+Cm), cbind(0.5%x%Am+Cm, Vm)), name="expCovDZM" ) expCovMZF <- mxAlgebra( expression= rbind( cbind(Vf, Af+Cf), cbind(Af+Cf, Vf)), name="expCovMZF" ) expCovDZF <- mxAlgebra( expression= rbind( cbind(Vf, 0.5%x%Af+Cf), cbind(0.5%x%Af+Cf, Vf)), name="expCovDZF" ) expCovDOSmf <- mxAlgebra( expression= rbind( cbind(Vm, 0.5%x%Am+Cm), cbind(0.5%x%Af+Cf, Vf)), name="expCovDOSmf" ) amCITOsc_ <- mxAlgebra(expression=Am[1,1]/Vm[1,1], name="amCITOsc_") cmCITOsc_ <- mxAlgebra(expression=Cm[1,1]/Vm[1,1], name="cmCITOsc_") emCITOsc_ <- mxAlgebra(expression=Em[1,1]/Vm[1,1], name="emCITOsc_") afCITOsc_ <- mxAlgebra(expression=Af[1,1]/Vf[1,1], name="afCITOsc_") cfCITOsc_ <- mxAlgebra(expression=Cf[1,1]/Vf[1,1], name="cfCITOsc_") efCITOsc_ <- mxAlgebra(expression=Ef[1,1]/Vf[1,1], name="efCITOsc_") amintellect_ <- mxAlgebra(expression=Am[2,2]/Vm[2,2], name="amintellect_") cmintellect_ <- mxAlgebra(expression=Cm[2,2]/Vm[2,2], name="cmintellect_") emintellect_ <- mxAlgebra(expression=Em[2,2]/Vm[2,2], name="emintellect_") afintellect_ <- mxAlgebra(expression=Af[2,2]/Vf[2,2], name="afintellect_") cfintellect_ <- mxAlgebra(expression=Cf[2,2]/Vf[2,2], name="cfintellect_") efintellect_ <- mxAlgebra(expression=Ef[2,2]/Vf[2,2], name="efintellect_") amopen_ <- mxAlgebra(expression=Am[3,3]/Vm[3,3], name="amopen_") cmopen_ <- mxAlgebra(expression=Cm[3,3]/Vm[3,3], name="cmopen_") emopen_ <- mxAlgebra(expression=Em[3,3]/Vm[3,3], name="emopen_") afopen_ <- mxAlgebra(expression=Af[3,3]/Vf[3,3], name="afopen_") cfopen_ <- mxAlgebra(expression=Cf[3,3]/Vf[3,3], name="cfopen_") efopen_ <- mxAlgebra(expression=Ef[3,3]/Vf[3,3], name="efopen_") amCint <- mxAlgebra(expression=Am[2,1]/Vm[2,1], name="amCint") cmCint <- mxAlgebra(expression=Cm[2,1]/Vm[2,1], name="cmCint") emCint <- mxAlgebra(expression=Em[2,1]/Vm[2,1], name="emCint") afCint <- mxAlgebra(expression=Af[2,1]/Vf[2,1], name="afCint") cfCint <- mxAlgebra(expression=Cf[2,1]/Vf[2,1], name="cfCint") efCint <- mxAlgebra(expression=Ef[2,1]/Vf[2,1], name="efCint") amint_open <- mxAlgebra(expression=Am[3,2]/Vm[3,2], name="amint_open") cmint_open <- mxAlgebra(expression=Cm[3,2]/Vm[3,2], name="cmint_open") emint_open <- mxAlgebra(expression=Em[3,2]/Vm[3,2], name="emint_open") afint_open <- mxAlgebra(expression=Af[3,2]/Vf[3,2], name="afint_open") cfint_open <- mxAlgebra(expression=Cf[3,2]/Vf[3,2], name="cfint_open") efint_open <- mxAlgebra(expression=Ef[3,2]/Vf[3,2], name="efint_open") amcor <- mxAlgebra(expression=solve(sqrt(I*Am)) %*% Am %*% solve(sqrt(I*Am)), name="amcor") cmcor <- mxAlgebra(expression=solve(sqrt(I*Cm)) %*% Cm %*% solve(sqrt(I*Cm)), name="cmcor") emcor <- mxAlgebra(expression=solve(sqrt(I*Em)) %*% Em %*% solve(sqrt(I*Em)), name="emcor") afcor <- mxAlgebra(expression=solve(sqrt(I*Af)) %*% Af %*% solve(sqrt(I*Af)), name="afcor") cfcor <- mxAlgebra(expression=solve(sqrt(I*Cf)) %*% Cf %*% solve(sqrt(I*Cf)), name="cfcor") efcor <- mxAlgebra(expression=solve(sqrt(I*Ef)) %*% Ef %*% solve(sqrt(I*Ef)), name="efcor") # Data objects for Multiple Groups dataMZM <- mxData( observed=mzmData, type="raw" ) dataDZM <- mxData( observed=dzmData, type="raw" ) dataMZF <- mxData( observed=mzfData, type="raw" ) dataDZF <- mxData( observed=dzfData, type="raw" ) dataDOSmf <- mxData( observed=dosmfData, type="raw" ) # Objectives MZMObjective <- mxFIMLObjective(covariance="expCovMZM", means="expMeanMZM", dimnames=selVars) DZMObjective <- mxFIMLObjective(covariance="expCovDZM", means="expMeanDZM", dimnames=selVars) MZFObjective <- mxFIMLObjective(covariance="expCovMZF", means="expMeanMZF", dimnames=selVars) DZFObjective <- mxFIMLObjective(covariance="expCovDZF", means="expMeanDZF", dimnames=selVars) DOSmfObjective <- mxFIMLObjective(covariance="expCovDOSmf", means="expMeanDOSmf", dimnames=selVars) # Combine Groups pars <- list( pathAm, pathCm, pathEm, pathAf, pathCf, pathEf, covAm, covCm, covEm, covAf, covCf, covEf, covVm,covVf, matI,iSDm,iSDf,corPhm,corPhf, grandMeanm,grandMeanf,expMeanMZM,expMeanDZM,expMeanMZF,expMeanDZF,expMeanDOSmf, amCITOsc_, cmCITOsc_, emCITOsc_, afCITOsc_, cfCITOsc_, efCITOsc_, amintellect_, cmintellect_, emintellect_, afintellect_, cfintellect_, efintellect_, amopen_, cmopen_, emopen_, afopen_, cfopen_, efopen_, amCint, cmCint, emCint, afCint, cfCint, efCint, amint_open, cmint_open, emint_open, afint_open, cfint_open, efint_open, amcor, cmcor, emcor, afcor, cfcor, efcor) # Models MZMmodel <- mxModel("MZMmodel", dataMZM, pars, expCovMZM, MZMObjective) DZMmodel <- mxModel("DZMmodel", dataDZM, pars, expCovDZM, DZMObjective) MZFmodel <- mxModel("MZFmodel", dataMZF, pars, expCovMZF, MZFObjective) DZFmodel <- mxModel("DZFmodel", dataDZF, pars, expCovDZF, DZFObjective) DOSmfmodel <- mxModel("DOSmfmodel", dataDOSmf, pars, expCovDOSmf, DOSmfObjective) # Objective min2sumll <- mxAlgebra( expression = MZMmodel.objective + DZMmodel.objective + MZFmodel.objective + DZFmodel.objective + DOSmfmodel.objective , name="min2sumll" ) objective <- mxAlgebraObjective("min2sumll") # Cholesky ACE model CholACEModel <- mxModel("Chol", MZMmodel, DZMmodel, MZFmodel, DZFmodel, DOSmfmodel, min2sumll, objective) # CI added here ## Fitting CholACEFit <- mxRun(CholACEModel, intervals=F) # Generate ACE Output # ----------------------------------------------------------------------- summary(CholACEFit) parameterSpecifications(CholACEFit) expectedMeansCovariances(CholACEFit) tableFitStatistics(CholACEFit) ## ACE estimates estimatedAm <- CholACEFit@submodels$MZMmodel@algebras$Am@result/CholACEFit@submodels$MZMmodel@algebras$Vm@result estimatedCm <- CholACEFit@submodels$MZMmodel@algebras$Cm@result/CholACEFit@submodels$MZMmodel@algebras$Vf@result estimatedEm <- CholACEFit@submodels$MZMmodel@algebras$Em@result/CholACEFit@submodels$MZMmodel@algebras$Vm@result estimatedAf <- CholACEFit@submodels$MZMmodel@algebras$Af@result/CholACEFit@submodels$MZMmodel@algebras$Vf@result estimatedCf <- CholACEFit@submodels$MZMmodel@algebras$Cf@result/CholACEFit@submodels$MZMmodel@algebras$Vf@result estimatedEf <- CholACEFit@submodels$MZMmodel@algebras$Ef@result/CholACEFit@submodels$MZMmodel@algebras$Vf@result print(estimatedAm) print(estimatedCm) print(estimatedEm) print(estimatedAf) print(estimatedCf) print(estimatedEf) ## Correlations correlationsAm <- solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Am@result)) %*% CholACEFit@submodels$MZMmodel@algebras$Am@result %*% solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Am@result)) correlationsCm <- solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Cm@result)) %*% CholACEFit@submodels$MZMmodel@algebras$Cm@result %*% solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Cm@result)) correlationsEm <- solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Em@result)) %*% CholACEFit@submodels$MZMmodel@algebras$Em@result %*% solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Em@result)) correlationsAf <- solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Af@result)) %*% CholACEFit@submodels$MZMmodel@algebras$Af@result %*% solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Af@result)) correlationsCf <- solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Cf@result)) %*% CholACEFit@submodels$MZMmodel@algebras$Cf@result %*% solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Cf@result)) correlationsEf <- solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Ef@result)) %*% CholACEFit@submodels$MZMmodel@algebras$Ef@result %*% solve(sqrt(CholACEFit@submodels$MZMmodel@matrices$I@values*CholACEFit@submodels$MZMmodel@algebras$Ef@result)) print(correlationsAm) print(correlationsCm) print(correlationsEm) print(correlationsAf) print(correlationsCf) print(correlationsEf) #------------------------------------------------------ #Cholesky Ace decomposition with c21m, c22m, c32m = 0 #------------------------------------------------------ CholAceModel <- mxModel(CholACEModel) CholAceModel <- omxSetParameters(CholAEModel, labels=c("c21m","c22m","c32m"), free=F, values = c( 0, 0, 0) ) CholAceModel <- omxSetParameters(CholAEModel, labels=c("c11m","c31m","c33m"), free=T, values =c(1, .8, 1 ) CholAceFit <- mxRun(CholAceModel) # Generate Ace Output # ----------------------------------------------------------------------- summary(CholAceFit) parameterSpecifications(CholAceFit) expectedMeansCovariances(CholAceFit) tableFitStatistics(CholAEFit) ### CholAceFitLLL <- CholAceFit@output$Minus2LogLikelihood CholACEFitLLL <- CholACEFit@output$Minus2LogLikelihood chi2CholAce <- CholAceFitLLL-CholACEFitLLL p_CholAce <- pchisq( chi2CholAce, lower.tail=F, 2) chi2CholAce; p_CholAce ## ACE estimates estimatedAm <- CholAceFit@submodels$MZMmodel@algebras$Am@result/CholACEFit@submodels$MZMmodel@algebras$Vm@result estimatedCm <- CholAceFit@submodels$MZMmodel@algebras$Cm@result/CholACEFit@submodels$MZMmodel@algebras$Vm@result estimatedEm <- CholAceFit@submodels$MZMmodel@algebras$Em@result/CholACEFit@submodels$MZMmodel@algebras$Vm@result print(estimatedAm) print(estimatedCm) print(estimatedEm)