# ------------------------------------------------------------------------------ # Program: twinMulAceCon.R # Author: Hermine Maes # Date: 02 25 2012 # # Edited: Katrina Grasby, May 2012, Aug 2015 # # Common Pathway Multivariate Twin Analysis model # Matrix style model - Raw data - Continuous data # ------------------------------------------------------------------------------ setwd("H:/Articles 2016/Twins/Oct 16") rm(list=ls(all=TRUE)) # Load Libraries require(psych) require(OpenMx) #omxGetNPSOL() #mxOption(NULL,"Default optimizer","NPSOL") mxOption(NULL,"Default optimizer","SLSQP") #THIS ONE WORKS!! mxOption(NULL,"Default optimizer","NPSOL") mxOption(NULL, 'Default optimizer', 'CSOLNP') #SLSQP is the one that seems to work setwd("H:/Articles 2016/Twins/Oct 16") omxGetNPSOL() detach('package:OpenMx',unload=TRUE) source('http://openmx.psyc.virginia.edu/getOpenMx.R') # -------------------------------------------------------------------- # PREPARE DATA twinx <- read.csv(file = "twin40.csv", header = TRUE, sep = ",") names(twinx) # Select Variables for Analysis. Use these lists to subset only the data that we want to use. nv <- 3 # number of variables nt <- 2 # number of people ntv <- nv*nt # number of total variables Vars <- c('ed','occ','inc') selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") defVars <- c('sex1', 'sex2', 'age1', 'age2') useVars <- c("zyg", selVars,defVars) #data; twinx$sex1 <- 2-twinx$gender_cde_twin1 twinx$sex2 <- 2-twinx$gender_cde_twin1 twinx$age1 <- twinx$C_age_twin1 twinx$age2 <- twinx$C_age_twin2 twinx$ed1 <- twinx$C_yrs_ed_twin1 twinx$ed2 <- twinx$C_yrs_ed_twin2 twinx$occ1 <- twinx$C_AUSEI_twin1 twinx$occ2 <- twinx$C_AUSEI_twin2 twinx$inc1 <- (twinx$C_inc_twin1)/10000 twinx$inc2 <- (twinx$C_inc_twin2)/10000 #twinx$pared1 <- twinx$yrs_ed_par_twin1 #twinx$pared2 <- twinx$yrs_ed_par_twin2 #twinx$pared1 <- ifelse(!is.na(twinx$yrs_ed_par_twin1), twinx$yrs_ed_par_twin1, twinx$yrs_ed_par_twin2) #twinx$pared2 <- ifelse(!is.na(twinx$yrs_ed_par_twin2), twinx$yrs_ed_par_twin2, twinx$yrs_ed_par_twin1) #help(ifelse) #twinx <- twinx[!(is.na(twinx$ed1)) & !(is.na(twinx$ed2)),] #removes missing values #twinx$pared1[is.na(twinx$pared1)] = mean(twinx$pared1, na.rm=TRUE) #substitutes means for missing #twinx$pared2[is.na(twinx$pared2)] = mean(twinx$pared2, na.rm=TRUE) #substitutes means for missing mean(twinx$inc2, na.rm=T) mean(twinx$ed2) # ------------------------------------------------------------------------------ # PREPARE MEANS MODEL # Regression effects # Matrices to hold beta weights beta1 <- mxMatrix( type='Full', nrow=1, ncol=nv, free=FALSE, values=0, labels=c('bSexEd','bSexOcc','bSexInc'), name='beta1' ) beta2 <- mxMatrix( type='Full', nrow=1, ncol=nv, free=TRUE, values=c(-0.04, 0, -0.02, -0.04, 0, -0.02), labels=c('bAgeEd','bAgeOcc','bAgeInc'), name='beta2' ) # Matrices to hold covariate data obsSex <- mxMatrix( type='Full', nrow=1, ncol=nt, free=FALSE, labels=c('data.sex1','data.sex2'), name='Sex') obsAge <- mxMatrix( type='Full', nrow=1, ncol=nt, free=FALSE, labels=c('data.age1','data.age2'), name='Age') # Matrices to hold intercept/mean intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=c(T, T, T, T, T, T), values=0, labels=c('IntEd', 'IntOcc','IntInc', 'IntEd', 'IntOcc','IntInc'), name="intercept") #intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=6, labels="mean", name="intercept" ) # Regression for mean adjusted for covariates expMean <- mxAlgebra(intercept + Sex%x%beta1 + Age%x%beta2, name="expMean") # List of objects for Means Model inclusions <- list (obsSex, beta1, obsAge, beta2, intercept, expMean) var_ed= var(twinx$ed1, na.rm=T) st=(var_ed/3)**0.5 st_a_ed=((var_ed)*0.44)**0.5 st_c_ed=((var_ed)*0.15)**0.5 st_e_ed=((var_ed)*0.41)**0.5 var_occ= var(twinx$occ1, na.rm=T) st=(var_occ/3)**0.5 st_a_occ=((var_occ)*0.36)**0.5 st_c_occ=((var_occ)*0.02)**0.5 st_e_occ=((var_occ)*0.62)**0.5 var_inc= var(twinx$inc1, na.rm=T) st=(var_inc/3)**0.5 st_a_inc=((var_inc)*0.26)**0.5 st_c_inc=((var_inc)*0.18)**0.5 st_e_inc=((var_inc)*0.56)**0.5 # ------------------------------------------------------------------------------ # DATA # Set up the groups # Data objects for Multiple Groups mzfData <- subset(twinx, zyg==1, useVars) dataMZf <- mxData( observed=mzfData, type="raw" ) mzmData <- subset(twinx, zyg==2, useVars) dataMZm <- mxData( observed=mzmData, type="raw" ) nrow(mzmData) dzfData <- subset(twinx, zyg==3, useVars) dataDZf <- mxData( observed=dzfData, type="raw" ) dzmData <- subset(twinx, zyg==4, useVars) dataDZm <- mxData( observed=dzmData, type="raw" ) nrow(dzmData) dzfmData <- subset(twinx, zyg==5, useVars) dataDZfm <- mxData( observed=dzfmData, type="raw" ) dzmfData <- subset(twinx, zyg==6, useVars) dataDZmf <- mxData( observed=dzmfData, type="raw" ) mean(mzData$ed1, na.rm=T) mean(mzData$ed2, na.rm=T) mean(mzData$occ1, na.rm=T) mean(mzData$occ2, na.rm=T) mean(mzData$inc1, na.rm=T) mean(mzData$inc2, na.rm=T) cor(twinx$occ1, occ2) cor(mzData,use="pairwise") s1 <- subset(twinx, select=c("ed1", "occ1", "inc1")) cov(s1,use="complete") s2 <- subset(twinx, select=c("ed2", "occ2", "inc2")) cov(s2,use="complete", decimals=2) help(cov) # ------------------------------------------------------------------------------ # Cholesky Decomposition ACE Model # ------------------------------------------------------------------------------ # Create Labels for Lower Triangular Matrices afLabs <- paste("af",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") cfLabs <- paste("cf",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") efLabs <- paste("ef",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") amLabs <- paste("am",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") cmLabs <- paste("cm",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") emLabs <- paste("em",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") help(rev) # Matrices declared to store a, c, and e Path Coefficients no ogg elements in lower triangle ((n+1) x n)/2 pathAf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T, T, T, T, T, T), values=c(-1.8, -9.2, -0.5, 5.4, 0.57, 0.91), label=afLabs, name="af" ) pathCf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T, T, T, T, T, T), values=c(0.93, 4.51, 1.29, 4.37, 0.795, -1.27), label=cfLabs, name="cf") pathEf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T, T, T, T, T, T), values=c(1.71, 4.77, 0.48, -15.27, -0.485, 2.92), label=efLabs, name="ef" ) pathAm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T, T, T, T, T, T), values=c(-1.95, -10.83, -0.537, 5.64, 0.786, 0.275), label=amLabs, name="am" ) pathCm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T, T, T, T, T, T), values=c(0.58, -1.29, 1.32, -2.059, -0.757, -2.43), label=cmLabs, name="cm" ) pathEm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T, T, T, T, T, T), values=c(1.79, 5.53, 0.41, 15.77, 0.67, -3.04), label=emLabs, name="em" ) pathRg <- mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=1, label="rg11", name="rg" , lbound=0, ubound=1) # Matrices generated to hold A, C, and E computed Variance Components covAf <- mxAlgebra( af %*% t(af), name="Af" ) covCf <- mxAlgebra( cf %*% t(cf), name="Cf" ) covEf <- mxAlgebra( ef %*% t(ef), name="Ef" ) covAm <- mxAlgebra( am %*% t(am), name="Am" ) covCm <- mxAlgebra( cm %*% t(cm), name="Cm" ) covEm <- mxAlgebra( em %*% t(em), name="Em" ) # Algebra to compute total variances and standard deviations (diagonal only) covPf <- mxAlgebra( Af+Cf+Ef, name="Vf" ) covPm <- mxAlgebra( Am+Cm+Em, name="Vm" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSDf <- mxAlgebra( solve(sqrt(I*Vf)), name="iSDf") invSDm <- mxAlgebra( solve(sqrt(I*Vm)), name="iSDm") # Matrices to hold standardised a, c, and e stpathaf <- mxAlgebra( iSDf %*% af, name="staf") stpathcf <- mxAlgebra( iSDf %*% cf, name="stcf") stpathef <- mxAlgebra( iSDf %*% ef, name="stef") stpatham <- mxAlgebra( iSDm %*% am, name="stam") stpathcm <- mxAlgebra( iSDm %*% cm, name="stcm") stpathem <- mxAlgebra( iSDm %*% em, name="stem") # Matrices to hold standardised A, C, and E stcovAf <- mxAlgebra( Af/Vf, name="stAf") stcovCf <- mxAlgebra( Cf/Vf, name="stCf") stcovEf <- mxAlgebra( Ef/Vf, name="stEf") stcovAm <- mxAlgebra( Am/Vm, name="stAm") stcovCm <- mxAlgebra( Cm/Vm, name="stCm") stcovEm <- mxAlgebra( Em/Vm, name="stEm") # Algebra to compute correlations corAf <- mxAlgebra(expression = solve(sqrt(I*Af))%&%Af, name ="rAf") corCf <- mxAlgebra(expression = solve(sqrt(I*Cf))%&%Cf, name ="rCf") corEf <- mxAlgebra(expression = solve(sqrt(I*Ef))%&%Ef, name ="rEf") corAm <- mxAlgebra(expression = solve(sqrt(I*Am))%&%Am, name ="rAm") corCm <- mxAlgebra(expression = solve(sqrt(I*Cm))%&%Cm, name ="rCm") corEm <- mxAlgebra(expression = solve(sqrt(I*Em))%&%Em, name ="rEm") corPf <- mxAlgebra(expression = iSDf%&%Vf, name ="rPf") corPm <- mxAlgebra(expression = iSDm%&%Vm, name ="rPm") # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins covMZf <- mxAlgebra( rbind( cbind(Af+Cf+Ef , Af+Cf), cbind(Af+Cf , Af+Cf+Ef)), name="expCovMZf" ) covDZf <- mxAlgebra( rbind( cbind(Af+Cf+Ef , 0.5%x%Af+Cf), cbind(0.5%x%Af+Cf , Af+Cf+Ef)), name="expCovDZf" ) covMZm <- mxAlgebra( rbind( cbind(Am+Cm+Em , Am+Cm), cbind(Am+Cm , Am+Cm+Em)), name="expCovMZm" ) covDZm <- mxAlgebra( rbind( cbind(Am+Cm+Em , 0.5%x%Am+Cm), cbind(0.5%x%Am+Cm , Am+Cm+Em)), name="expCovDZm" ) covDZfm <- mxAlgebra( rbind( cbind(Af+Cf+Ef , 0.5%*%rg%x%(af%*%t(am))+cf%*%t(cm)), cbind(0.5%*%rg%x%(am%*%t(af))+cm%*%t(cf) , Am+Cm+Em)), name="expCovDZfm" ) covDZmf <- mxAlgebra( rbind( cbind(Am+Cm+Em , 0.5%*%rg%x%(af%*%t(am))+cf%*%t(cm)), cbind(0.5%*%rg%x%(am%*%t(af))+cm%*%t(cf) , Af+Cf+Ef)), name="expCovDZmf" ) # Objective objects for Multiple Groups fitfun <- mxFitFunctionML() expMZf <- mxExpectationNormal( covariance="expCovMZf", means="expMean", dimnames=selVars ) expDZf <- mxExpectationNormal( covariance="expCovDZf", means="expMean", dimnames=selVars ) expMZm <- mxExpectationNormal( covariance="expCovMZm", means="expMean", dimnames=selVars ) expDZm <- mxExpectationNormal( covariance="expCovDZm", means="expMean", dimnames=selVars ) expDZfm <- mxExpectationNormal( covariance="expCovDZfm", means="expMean", dimnames=selVars ) expDZmf <- mxExpectationNormal( covariance="expCovDZmf", means="expMean", dimnames=selVars ) # Combine Groups parsZf <- list( pathAf, pathCf, pathEf, covAf, covCf, covEf, covPf, matI, invSDf, stpathaf, stpathcf, stpathef, stcovAf, stcovCf, stcovEf, corAf, corCf, corEf, corPf) parsZm <- list( pathAm, pathCm, pathEm, covAm, covCm, covEm, covPm, matI, invSDm, stpatham, stpathcm, stpathem, stcovAm, stcovCm, stcovEm, corAm, corCm, corEm, corPm) modelMZf <- mxModel( parsZf, inclusions, covMZf, dataMZf, expMZf, fitfun, name="MZf" ) modelDZf <- mxModel( parsZf, inclusions, covDZf, dataDZf, expDZf, fitfun, name="DZf" ) modelMZm <- mxModel( parsZm, inclusions, covMZm, dataMZm, expMZm, fitfun, name="MZm" ) modelDZm <- mxModel( parsZm, inclusions, covDZm, dataDZm, expDZm, fitfun, name="DZm" ) modelDZfm <- mxModel( parsZf, pathRg, parsZm, inclusions, covDZfm, dataDZfm, expDZfm, fitfun, name="DZfm" ) modelDZmf <- mxModel( parsZf, pathRg, parsZm, inclusions, covDZmf, dataDZmf, expDZmf, fitfun, name="DZmf" ) multi <- mxFitFunctionMultigroup( c("MZf","MZm","DZf","DZm","DZfm","DZmf")) CIsf <- mxCI( c('stAf','stCf','stEf')) CIsm <- mxCI( c('stAm','stCm','stEm')) modelChol <- mxModel( "CholACE", parsZf, parsZm, modelMZf, modelMZm, modelDZf, modelDZm, modelDZfm, modelDZmf, multi, CIsf, CIsm) # Run model fitCholAce <- mxRun(modelChol, intervals=F) sumCholAce <- summary(fitCholAce) #HERE HERE HERE #Modify Cholesky CholMod <- mxModel(fitCholAce, name="Mod Chol" ) #BEST MODEL FOLLOWS #THis is MINIMAL C THIS ONE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CholMod <- omxSetParameters(CholMod, labels=c( "af_3_3", "cf_2_1", "cf_2_2", "cf_3_2", "cf_3_3", "am_3_3", "cm_2_1", "cm_2_2", "cm_3_2", "cm_3_3") , free=FALSE, values=0) #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #"af_3_1" "am_3_1" #THis is Model 5, THIS IS IT NOV 11 PREFERRED MODEL CholMod <- omxSetParameters(CholMod, labels=c("af_3_3", "af_3_1", "cf_2_1","cf_2_2", "cf_3_2", "cf_3_3", "am_3_3", "am_3_1", "cm_2_1","cm_2_2", "cm_3_2", "cm_3_3") , free=FALSE, values=0) CholMod <- omxSetParameters(CholMod, labels=c("cf_1_1", "cf_2_1", "cf_3_1", "cf_2_2", "cf_3_2","cf_3_3", "cm_1_1", "cm_2_1", "cm_3_1", "cm_2_2", "cm_3_2","cm_3_3"), free=FALSE, values=0) CholMod <- omxSetParameters(CholMod, labels=c("a_1_1", "a_2_1", "a_3_1", "a_2_2", "a_3_2", "a_3_3") ,free=FALSE, values=0) #Test PREFERRED MODEL NOW CholMod <- omxSetParameters(CholMod, labels=c("af_2_2", "af_3_2", "af_3_3", "am_2_2", "am_3_2", "am_3_3", "cf_2_1", "cf_2_2", "cf_3_1", "cf_3_2", "cm_2_1", "cm_2_2", "cm_3_1", "cm_3_2") , free=FALSE, values=0) #CholMod <- omxSetParameters(CholMod, labels=c("a_3_2", "c_1_1", "c_2_1", "c_3_1", "c_2_2", "c_3_2") , free=FALSE, values=0) #THis is Model 5, THIS IS IT, TUES NOV 8 NOV 10 CholMod <- omxSetParameters(CholMod, labels=c("af_3_1", "cf_2_1","cf_2_2", "cf_3_2", "cf_3_3", "am_3_1", "cm_2_1","cm_2_2", "cm_3_2", "cm_3_3") , free=FALSE, values=0) #This Model 6, No Suspect SEs; CholMod <- omxSetParameters(CholMod, labels=c("cf_2_1","cf_2_2", "cf_3_2", "cf_3_3", "af_3_3", "cf_1_1", "cm_2_1","cm_2_2", "cm_3_2", "cm_3_3", "am_2_2", "am_3_3", "am_3_2") , free=FALSE, values=0) #Model in the Paper CholMod <- omxSetParameters(CholMod, labels=c("af_2_2", "af_3_2", "cf_2_1", "cf_2_2", "cf_3_2", "cf_3_1", "am_2_2", "am_3_2", "cm_2_1", "cm_2_2", "cm_3_2", "cm_3_1" ), free=FALSE, values=0) #"af_2_2", "af_3_2", , #"cf_2_1","cf_2_2", "cf_3_2", CholMod <- omxSetParameters(CholMod, labels=c("a_3_3", "c_3_1", "c_2_1", "c_2_2", "c_3_2") , free=FALSE, values=0) #CholMod <- omxSetParameters(CholMod, labels=c("c_2_1", "c_2_2", "c_3_2", "c_3_3") , free=FALSE, values=0) CholMod <- omxSetParameters(CholMod, labels=c("c_1_1", "c_2_1", "c_3_1", "c_2_2", "c_3_2", "a_2_2", "a_3_2", "a_3_3") , free=FALSE, values=0) #CholMod2 <- omxSetParameters(CholMod, labels=c("e_2_1", "e_3_1", "e_3_2") , free=FALSE, values=0) FitCholMod <- mxRun(CholMod, intervals =FALSE) SumCholMod <- summary(FitCholMod, verbose=FALSE) FitCholMod$output$status$code SumCholMod$seSuspect mxCompare(fitCholAce,FitCholMod) print(FitCholMod@output$algebras) FitCholMod$af FitCholMod$Vf FitCholMod$Vm FitCholMod$Af FitCholMod$Am FitCholMod$Em FitCholMod$C FitCholMod$I FitCholMod$stAf FitCholMod$stAm FitCholMod$stCf FitCholMod$stCm FitCholMod$stEf FitCholMod$stEm FitCholMod$staf FitCholMod$stam FitCholMod$stcf FitCholMod$stcm FitCholMod$stef FitCholMod$stem FitCholMod$rAf FitCholMod$rAm FitCholMod$rCf FitCholMod$rE FitCholMod$stc FitCholMod$ste FitCholMod$rP estCovMZ <- mxEval(expCovMZ, FitCholMod$MZ) estcorC <- mxEval(expCovDZ, FitCholMod$C) corC <- mxAlgebra( solve(sqrt(I*C))%&%C, name ="rC")