R version 3.4.4 (2018-03-15) -- "Someone to Lean On" Copyright (C) 2018 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(psych) > library(OpenMx) > mxOption(NULL, 'Number of Threads', parallel::detectCores()) #now > library(plyr) > > load("caddMJ_Behavioral_withOS.Rdat") > #convert to wide > cadd$FREQ1 <- ifelse(cadd$age1 <= 17 & cadd$age1 >= 13, cadd$FREQ1, NA) > caddT <- subset(cadd, select=c("famid", "twin", "sex", "cohort", "zygosity", "ADHDres", "ADres", "CDres", "age1", "FREQ1", "age2", "FREQ2", "age3", "FREQ3", "pairgender")) > cadd.w <- reshape(caddT, v.names=c("ADHDres", "ADres", "CDres", "age1", "FREQ1", "age2", "FREQ2", "age3", "FREQ3"), timevar="twin", idvar="famid", direction="wide", sep="") > > > ### Create dummy intercept definition variable and dummy 0 variable > cadd.w$intercept <- 1 > cadd.w$zero <- 0 > > #format age and definition variables > cadd.w$age1avg <- (cadd.w$age10 + cadd.w$age11)/2 > cadd.w$age2avg <- (cadd.w$age20 + cadd.w$age21)/2 > cadd.w$age3avg <- (cadd.w$age30 + cadd.w$age31)/2 > > cadd.w$age1twinpair <- ifelse(is.na(cadd.w$age10), cadd.w$age11, cadd.w$age1avg) > cadd.w$age1twinpair <- ifelse(is.na(cadd.w$age11), cadd.w$age10, cadd.w$age1avg) > summary(cadd.w$age10) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 11.58 12.64 14.47 14.84 16.81 19.01 340 > summary(cadd.w$age11) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 11.33 12.87 14.63 14.99 16.83 19.00 66 > summary(cadd.w$age1avg) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 11.58 12.64 14.47 14.84 16.81 19.00 341 > summary(cadd.w$age1twinpair) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 11.58 12.64 14.47 14.84 16.81 19.00 340 > > cadd.w$age2twinpair <- ifelse(is.na(cadd.w$age20), cadd.w$age21, cadd.w$age2avg) > cadd.w$age2twinpair <- ifelse(is.na(cadd.w$age21), cadd.w$age20, cadd.w$age2avg) > summary(cadd.w$age20) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.51 17.29 19.51 20.01 22.22 29.13 447 > summary(cadd.w$age21) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.13 17.73 19.77 20.11 22.14 29.09 224 > summary(cadd.w$age2avg) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.51 17.26 19.44 19.90 22.12 29.06 488 > summary(cadd.w$age2twinpair) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.51 17.29 19.57 20.00 22.24 29.06 447 > > cadd.w$age3twinpair <- ifelse(is.na(cadd.w$age30), cadd.w$age31, cadd.w$age3avg) > cadd.w$age3twinpair <- ifelse(is.na(cadd.w$age31), cadd.w$age30, cadd.w$age3avg) > summary(cadd.w$age30) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 21.15 22.93 24.99 25.34 27.53 34.36 441 > summary(cadd.w$age31) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 21.10 23.30 25.30 25.51 27.55 34.37 200 > summary(cadd.w$age3avg) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 21.14 22.88 24.91 25.27 27.43 34.36 492 > summary(cadd.w$age3twinpair) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 21.14 22.91 24.96 25.33 27.52 34.36 441 > > cadd.w$age1twinpair[is.na(cadd.w$age1twinpair)] <- mean(cadd.w$age1twinpair, na.rm=T) > cadd.w$age2twinpair[is.na(cadd.w$age2twinpair)] <- mean(cadd.w$age2twinpair, na.rm=T) > cadd.w$age3twinpair[is.na(cadd.w$age3twinpair)] <- mean(cadd.w$age3twinpair, na.rm=T) > > cadd.w$FREQ11[cadd.w$age1twinpair > 24] <- NA > cadd.w$FREQ21[cadd.w$age2twinpair > 24] <- NA > cadd.w$FREQ31[cadd.w$age3twinpair > 24] <- NA > > cadd.w$FREQ10[cadd.w$age1twinpair > 24] <- NA > cadd.w$FREQ20[cadd.w$age2twinpair > 24] <- NA > cadd.w$FREQ30[cadd.w$age3twinpair > 24] <- NA > > ###### > cadd.w$age1twinpairC <- cadd.w$age1twinpair - 16.5 > cadd.w$age2twinpairC <- cadd.w$age2twinpair - 16.5 > cadd.w$age3twinpairC <- cadd.w$age3twinpair - 16.5 > > cadd.w$age1twinpairC[cadd.w$age1twinpairC > 7.5] <- 7.5 > cadd.w$age2twinpairC[cadd.w$age2twinpairC > 7.5] <- 7.5 > cadd.w$age3twinpairC[cadd.w$age3twinpairC > 7.5] <- 7.5 > > ########################################################### > #load and format MTFS data > MJWide <- read.csv("mn_twin_mj_freq_wide.csv") > load("formatted_data_mtfs.RData") > MJWide <- subset(MJWide, select=c(2:12)) > MJWide$id <- MJWide$ID > dat.mj.mtfs <- join(mtfs.final, MJWide, type="full", by="id") > mtfs <- subset(dat.mj.mtfs, twin != 2, select=c("id", "famid", "twin", "sex", "cohort", "zygosity", "ADHDres", "ADres", "CDres", "mj_12m_freq_14", "age14", "mj_12m_freq_17", "age17", "mj_12m_freq_20", "age20", "mj_12m_freq_24", "age24")) > > #convert MN data to wide > mtfsT <- subset(mtfs, select=c("famid", "twin", "sex", "cohort", "zygosity", "ADHDres", "ADres", "CDres", "mj_12m_freq_14", "age14", "mj_12m_freq_17", "age17", "mj_12m_freq_20", "age20", "mj_12m_freq_24", "age24")) > mtfs.w <- reshape(mtfsT, v.names=c("ADHDres", "ADres", "CDres", "mj_12m_freq_14", "age14", "mj_12m_freq_17", "age17", "mj_12m_freq_20", "age20", "mj_12m_freq_24", "age24"), timevar="twin", idvar="famid", direction="wide", sep="") > > #format age and definiton variables > mtfs.w$age14avg <- (mtfs.w$age140 + mtfs.w$age141)/2 > mtfs.w$age17avg <- (mtfs.w$age170 + mtfs.w$age171)/2 > mtfs.w$age20avg <- (mtfs.w$age200 + mtfs.w$age201)/2 > mtfs.w$age24avg <- (mtfs.w$age240 + mtfs.w$age241)/2 > > > mtfs.w$age14twinpair <- ifelse(is.na(mtfs.w$age140), mtfs.w$age141, mtfs.w$age14avg) > mtfs.w$age14twinpair <- ifelse(is.na(mtfs.w$age141), mtfs.w$age140, mtfs.w$age14avg) > summary(mtfs.w$age140) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 13.57 14.56 14.90 14.90 15.22 17.02 678 > summary(mtfs.w$age141) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 13.57 14.56 14.90 14.90 15.22 16.72 678 > summary(mtfs.w$age14avg) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 13.57 14.56 14.90 14.90 15.21 16.84 685 > summary(mtfs.w$age14twinpair) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 13.57 14.56 14.91 14.90 15.22 16.84 678 > > mtfs.w$age17twinpair <- ifelse(is.na(mtfs.w$age170), mtfs.w$age171, mtfs.w$age17avg) > mtfs.w$age17twinpair <- ifelse(is.na(mtfs.w$age171), mtfs.w$age170, mtfs.w$age17avg) > summary(mtfs.w$age170) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.55 17.42 17.81 17.85 18.16 20.34 103 > summary(mtfs.w$age171) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.55 17.41 17.81 17.85 18.16 20.12 101 > summary(mtfs.w$age17avg) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.55 17.40 17.81 17.83 18.15 20.11 128 > summary(mtfs.w$age17twinpair) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 16.55 17.42 17.81 17.85 18.17 20.34 103 > > mtfs.w$age20twinpair <- ifelse(is.na(mtfs.w$age200), mtfs.w$age201, mtfs.w$age20avg) > mtfs.w$age20twinpair <- ifelse(is.na(mtfs.w$age201), mtfs.w$age200, mtfs.w$age20avg) > summary(mtfs.w$age200) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 19.22 20.56 20.96 21.07 21.53 24.29 479 > summary(mtfs.w$age201) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 19.45 20.57 20.96 21.08 21.53 24.29 484 > summary(mtfs.w$age20avg) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 19.45 20.57 20.96 21.06 21.49 24.29 520 > summary(mtfs.w$age20twinpair) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 19.22 20.58 20.97 21.07 21.53 24.29 479 > > mtfs.w$age24twinpair <- ifelse(is.na(mtfs.w$age240), mtfs.w$age241, mtfs.w$age24avg) > mtfs.w$age24twinpair <- ifelse(is.na(mtfs.w$age241), mtfs.w$age240, mtfs.w$age24avg) > summary(mtfs.w$age240) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 22.63 24.20 24.84 24.87 25.43 35.29 187 > summary(mtfs.w$age241) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 22.63 24.22 24.85 24.88 25.44 35.29 191 > summary(mtfs.w$age24avg) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 22.63 24.23 24.85 24.86 25.43 35.29 236 > summary(mtfs.w$age24twinpair) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 22.63 24.24 24.86 24.88 25.45 35.29 187 > > > mtfs.w$age14twinpair[is.na(mtfs.w$age14twinpair)] <- mean(mtfs.w$age14twinpair, na.rm=T) > mtfs.w$age17twinpair[is.na(mtfs.w$age17twinpair)] <- mean(mtfs.w$age17twinpair, na.rm=T) > mtfs.w$age20twinpair[is.na(mtfs.w$age20twinpair)] <- mean(mtfs.w$age20twinpair, na.rm=T) > mtfs.w$age24twinpair[is.na(mtfs.w$age24twinpair)] <- mean(mtfs.w$age24twinpair, na.rm=T) > > > mtfs.w$mj_12m_freq_140[mtfs.w$age14twinpair > 24] <- NA > mtfs.w$mj_12m_freq_170[mtfs.w$age17twinpair > 24] <- NA > mtfs.w$mj_12m_freq_200[mtfs.w$age20twinpair > 24] <- NA > mtfs.w$mj_12m_freq_240[mtfs.w$age24twinpair > 24] <- NA > > > mtfs.w$mj_12m_freq_141[mtfs.w$age14twinpair > 24] <- NA > mtfs.w$mj_12m_freq_171[mtfs.w$age17twinpair > 24] <- NA > mtfs.w$mj_12m_freq_201[mtfs.w$age20twinpair > 24] <- NA > mtfs.w$mj_12m_freq_241[mtfs.w$age24twinpair > 24] <- NA > > > mtfs.w$age1twinpairC <- mtfs.w$age14twinpair - 16.5 > mtfs.w$age2twinpairC <- mtfs.w$age17twinpair - 16.5 > mtfs.w$age3twinpairC <- mtfs.w$age20twinpair - 16.5 > mtfs.w$age4twinpairC <- mtfs.w$age24twinpair - 16.5 > > > mtfs.w$age1twinpairC[mtfs.w$age1twinpairC > 7.5] <- 7.5 > mtfs.w$age2twinpairC[mtfs.w$age2twinpairC > 7.5] <- 7.5 > mtfs.w$age3twinpairC[mtfs.w$age3twinpairC > 7.5] <- 7.5 > mtfs.w$age4twinpairC[mtfs.w$age4twinpairC > 7.5] <- 7.5 > > > ### Create dummy intercept definition variable and dummy 0 variable > mtfs.w$intercept <- 1 > mtfs.w$zero <- 0 > > summary(cadd.w$age1twinpairC) Min. 1st Qu. Median Mean 3rd Qu. Max. -4.919 -3.511 -1.662 -1.662 -0.398 2.498 > summary(cadd.w$age2twinpairC) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.01472 1.86584 3.50360 3.43387 4.48665 7.50000 > summary(cadd.w$age3twinpairC) Min. 1st Qu. Median Mean 3rd Qu. Max. 4.639 7.320 7.500 7.111 7.500 7.500 > > hist(cadd.w$age1twinpairC) > hist(cadd.w$age2twinpairC) > hist(cadd.w$age3twinpairC) > > cadd.w.MJ1 <- as.vector(as.matrix(cadd.w[,c("FREQ10", "FREQ11")])) > summary(cadd.w.MJ1) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.0000 0.0000 0.0000 0.2907 0.0000 7.0000 1773 > hist(cadd.w.MJ1) > mean(cadd.w.MJ1, na.rm=T) [1] 0.2907291 > sd(cadd.w.MJ1, na.rm=T) [1] 1.026274 > > cadd.w.MJ2 <- as.vector(as.matrix(cadd.w[,c("FREQ20", "FREQ21")])) > summary(cadd.w.MJ2) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.0000 0.0000 0.0000 0.8097 1.0000 7.0000 814 > hist(cadd.w.MJ2) > mean(cadd.w.MJ2, na.rm=T) [1] 0.8096618 > sd(cadd.w.MJ2, na.rm=T) [1] 1.709865 > > cadd.w.MJ3 <- as.vector(as.matrix(cadd.w[,c("FREQ30", "FREQ31")])) > summary(cadd.w.MJ3) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.000 0.000 0.000 1.071 1.000 7.000 2120 > hist(cadd.w.MJ3) > mean(cadd.w.MJ3, na.rm=T) [1] 1.070681 > sd(cadd.w.MJ3, na.rm=T) [1] 2.03385 > > > summary(mtfs.w$age1twinpairC) Min. 1st Qu. Median Mean 3rd Qu. Max. -2.9274 -1.7312 -1.5977 -1.5977 -1.4812 0.3387 > summary(mtfs.w$age2twinpairC) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04839 0.95699 1.34140 1.34986 1.63710 3.84140 > summary(mtfs.w$age3twinpairC) Min. 1st Qu. Median Mean 3rd Qu. Max. 2.721 4.241 4.573 4.573 4.774 7.500 > summary(mtfs.w$age4twinpairC) Min. 1st Qu. Median Mean 3rd Qu. Max. 6.126 7.500 7.500 7.432 7.500 7.500 > > hist(mtfs.w$age1twinpairC) > hist(mtfs.w$age2twinpairC) > hist(mtfs.w$age3twinpairC) > hist(mtfs.w$age4twinpairC) > > mtfs.w.MJ1 <- as.vector(as.matrix(mtfs.w[,c("mj_12m_freq_140", "mj_12m_freq_141")])) > summary(mtfs.w.MJ1) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.0000 0.0000 0.0000 0.1967 0.0000 6.0000 1398 > hist(mtfs.w.MJ1) > mean(mtfs.w.MJ1, na.rm=T) [1] 0.1966846 > sd(mtfs.w.MJ1, na.rm=T) [1] 0.7895143 > > mtfs.w.MJ2 <- as.vector(as.matrix(mtfs.w[,c("mj_12m_freq_170", "mj_12m_freq_171")])) > summary(mtfs.w.MJ2) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.0000 0.0000 0.0000 0.6432 0.0000 7.0000 267 > hist(mtfs.w.MJ2) > mean(mtfs.w.MJ2, na.rm=T) [1] 0.6431757 > sd(mtfs.w.MJ2, na.rm=T) [1] 1.578391 > > mtfs.w.MJ3 <- as.vector(as.matrix(mtfs.w[,c("mj_12m_freq_200", "mj_12m_freq_201")])) > summary(mtfs.w.MJ3) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.0000 0.0000 0.0000 0.8956 1.0000 7.0000 1014 > hist(mtfs.w.MJ3) > mean(mtfs.w.MJ3, na.rm=T) [1] 0.8956422 > sd(mtfs.w.MJ3, na.rm=T) [1] 1.822371 > > mtfs.w.MJ4 <- as.vector(as.matrix(mtfs.w[,c("mj_12m_freq_240", "mj_12m_freq_241")])) > summary(mtfs.w.MJ4) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 0.000 0.000 0.000 1.035 1.000 7.000 3143 > hist(mtfs.w.MJ4) > mean(mtfs.w.MJ4, na.rm=T) [1] 1.034908 > sd(mtfs.w.MJ4, na.rm=T) [1] 2.054007 > > ############################## > ######################## > ### Biometric Models ### > ######################## > ### > ### Select observed variables > VarsCO <- c("CDres", "ADres", "ADHDres", "FREQ1", "FREQ2", "FREQ3") > nvCO <- length(VarsCO) > ntvCO <- nvCO*2 > selVarsCO <- paste(VarsCO, c(rep(0,nvCO), rep(1,nvCO)), sep="") > > VarsMN <- c("CDres", "ADres", "ADHDres", "mj_12m_freq_14", "mj_12m_freq_17", "mj_12m_freq_20", "mj_12m_freq_24") > nvMN <- length(VarsMN) > ntvMN <- nvMN*2 > selVarsMN <- paste(VarsMN, c(rep(0,nvMN), rep(1,nvMN)), sep="") > > ############ > #definition variables > DefVarsCO <- c("age1twinpairC", "age2twinpairC", "age3twinpairC", "intercept", "sex", "zero") > DefVarsMN <- c("age1twinpairC", "age2twinpairC", "age3twinpairC", "age4twinpairC", "intercept", "sex", "zero") > > > > ### Select Data for Analysis > caddMzData <- subset(cadd.w, zygosity==1, select=c(selVarsCO, DefVarsCO)) > caddDzData <- subset(cadd.w, zygosity==2, select=c(selVarsCO, DefVarsCO)) > > caddDataMZ <- mxData( observed=caddMzData, type="raw" ) > caddDataDZ <- mxData( observed=caddDzData, type="raw" ) > > mtfsMzData <- subset(mtfs.w, zygosity==1, select=c(selVarsMN, DefVarsMN)) > mtfsDzData <- subset(mtfs.w, zygosity==2, select=c(selVarsMN, DefVarsMN)) > > mtfsDataMZ <- mxData( observed=mtfsMzData, type="raw" ) > mtfsDataDZ <- mxData( observed=mtfsDzData, type="raw" ) > > ################################################################### > #base model > ################################################################### > > nl <- 3 # number of latent factors > # Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s) > XMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, 0), + labels=c("x11MN", + "x21MN", "x22MN", + "x31MN", "x32MN", "x33MN"), name="XMN" ) > > YMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("y11MN", + "y21MN", "y22MN", + "y31MN", "y32MN", "y33MN"), name="YMN" ) > > ZMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("z11MN", + "z21MN", "z22MN", + "z31MN", "z32MN", "z33MN"), name="ZMN" ) > > XCO <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, 0), + labels=c("x11CO", + "x21CO", "x22CO", + "x31CO", "x32CO", "x33CO"), name="XCO" ) > > YCO <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("y11CO", + "y21CO", "y22CO", + "y31CO", "y32CO", "y33CO"), name="YCO" ) > > ZCO <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("z11CO", + "z21CO", "z22CO", + "z31CO", "z32CO", "z33CO"), name="ZCO" ) > > > AlMN <- mxAlgebra(XMN %*% t(XMN), name="AlMN") > ClMN <- mxAlgebra(YMN %*% t(YMN), name="ClMN") > ElMN <- mxAlgebra(ZMN %*% t(ZMN), name="ElMN") > > AlCO <- mxAlgebra(XCO %*% t(XCO), name="AlCO") > ClCO <- mxAlgebra(YCO %*% t(YCO), name="ClCO") > ElCO <- mxAlgebra(ZCO %*% t(ZCO), name="ElCO") > > VarLMN <- mxAlgebra(expression = AlMN + ClMN + ElMN, name="VarLMN") > VarLCO <- mxAlgebra(expression = AlCO + ClCO + ElCO, name="VarLCO") > > corLMN <- mxAlgebra(cov2cor(VarLMN), name="corLMN") > corLCO <- mxAlgebra(cov2cor(VarLCO), name="corLCO") > > # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors > pathAsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3, + labels=c("a11MN", "a22MN", "a33MN", "a44MN", "a55MN", "a66MN", "a77MN"), + lbound=.00001, name="asMN" ) > pathCsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3, + labels=c("c11MN","c22MN", "c33MN", "c44MN", "c55MN", "c66MN", "c77MN"), + lbound=.00001, name="csMN" ) > pathEsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3, + labels=c("e11MN", "e22MN", "e33MN", "e44MN", "e55MN", "e66MN", "e77MN"), + lbound=.00001, name="esMN" ) > > pathAsCO <- mxMatrix(type="Diag", nrow=nvCO, ncol=nvCO, free=TRUE, values=.3, + labels=c("a11CO", "a22CO", "a33CO", "a44CO", "a55CO", "a66CO"), + lbound=.00001, name="asCO" ) > pathCsCO <- mxMatrix(type="Diag", nrow=nvCO, ncol=nvCO, free=TRUE, values=.3, + labels=c("c11CO","c22CO", "c33CO", "c44CO", "c55CO", "c66CO"), + lbound=.00001, name="csCO" ) > pathEsCO <- mxMatrix(type="Diag", nrow=nvCO, ncol=nvCO, free=TRUE, values=.3, + labels=c("e11CO", "e22CO", "e33CO", "e44CO", "e55CO", "e66CO"), + lbound=.00001, name="esCO" ) > > > > > pathFlGrowthMN <- mxMatrix(type="Full", nrow=4, ncol=nl, free=FALSE, + labels=c("data.zero", "data.intercept", "data.age1twinpairC", + "data.zero", "data.intercept", "data.age2twinpairC", + "data.zero", "data.intercept", "data.age3twinpairC", + "data.zero", "data.intercept", "data.age4twinpairC"), name="flGrowthMN", byrow=T) > > pathFlGrowthCO <- mxMatrix(type="Full", nrow=3, ncol=nl, free=FALSE, + labels=c("data.zero", "data.intercept", "data.age1twinpairC", + "data.zero", "data.intercept", "data.age2twinpairC", + "data.zero", "data.intercept", "data.age3twinpairC"), name="flGrowthCO", byrow=T) > > > flFree <- c(FALSE, FALSE, FALSE, + TRUE, FALSE, FALSE, + TRUE, FALSE, FALSE) > > flValues <- c(1, 0, 0, + .6, 0, 0, + .6, 0, 0) > > #set first loading to 1 and make it fixed for BD to manifest variables > pathFlBDMN <- mxMatrix(type="Full", nrow=3, ncol=nl, free=flFree, + values=flValues, lbound=.000001, labels=c("fBD1MN", NA, NA, + "fBD2MN", NA, NA, + "fBD3MN", NA, NA), name="flBDMN", byrow=T) > > pathFlBDCO <- mxMatrix(type="Full", nrow=3, ncol=nl, free=flFree, + values=flValues, lbound=.000001, labels=c("fBD1CO", NA, NA, + "fBD2CO", NA, NA, + "fBD3CO", NA, NA), name="flBDCO", byrow=T) > > #combine factor loading matrices > pathflMN <- mxAlgebra(rbind(flBDMN, flGrowthMN), name="flMN") > > pathflCO <- mxAlgebra(rbind(flBDCO, flGrowthCO), name="flCO") > > covAMN <- mxAlgebra( expression=flMN %*% AlMN %*% t(flMN) + asMN %*% t(asMN), name="AMN" ) > covCMN <- mxAlgebra( expression=flMN %*% ClMN %*% t(flMN) + csMN %*% t(csMN), name="CMN" ) > covEMN <- mxAlgebra( expression=flMN %*% ElMN %*% t(flMN) + esMN %*% t(esMN), name="EMN" ) > > covACO <- mxAlgebra( expression=flCO %*% AlCO %*% t(flCO) + asCO %*% t(asCO), name="ACO" ) > covCCO <- mxAlgebra( expression=flCO %*% ClCO %*% t(flCO) + csCO %*% t(csCO), name="CCO" ) > covECO <- mxAlgebra( expression=flCO %*% ElCO %*% t(flCO) + esCO %*% t(esCO), name="ECO" ) > > matIMN <- mxMatrix( type="Iden", nrow=3, ncol=3, name="IMN") > covPMN <- mxAlgebra( expression=AMN[1:3, 1:3]+CMN[1:3, 1:3]+EMN[1:3, 1:3], name="VMN" ) > invSDMN <- mxAlgebra( expression=solve(sqrt(IMN*VMN)), name="iSDMN") > > matICO <- mxMatrix( type="Iden", nrow=3, ncol=3, name="ICO") > covPCO <- mxAlgebra( expression=ACO[1:3, 1:3]+CCO[1:3, 1:3]+ECO[1:3, 1:3], name="VCO" ) > invSDCO <- mxAlgebra( expression=solve(sqrt(ICO*VCO)), name="iSDCO") > > #estimate means > #estimate means for the intercept and slope > FacMeansMN <- mxMatrix("Full",nrow=1,ncol=nl, free=c(TRUE, TRUE, TRUE), values=c(0,.8,.3), labels=c("mean1MN", "mean2MN", "mean3MN"), name="FacMeansMN") > regCoefMN <- mxMatrix(type="Full", nrow=1, ncol=nl, free=TRUE, values=0, + labels=c("beta1MN", "beta2MN", "beta3MN"), name="betaMN") > sexCovMN <- mxMatrix(type="Full", nrow=1, ncol=nl, free=FALSE, labels="data.sex", name="sexMN") > FacMeans2MN <- mxAlgebra(expression= t(FacMeansMN+(betaMN*sexMN)), name="FacMeans2MN") > meansMN <- mxAlgebra(t(flMN %*% FacMeans2MN), name="meansMN") > expMeanMN <- mxAlgebra(cbind(meansMN, meansMN), name="expMeanMN") > > > FacMeansCO <- mxMatrix("Full",nrow=1,ncol=nl, free=c(TRUE, TRUE, TRUE), values=c(0,.8,.3), labels=c("mean1CO", "mean2CO", "mean3CO"), name="FacMeansCO") > regCoefCO <- mxMatrix(type="Full", nrow=1, ncol=nl, free=TRUE, values=0, + labels=c("beta1CO", "beta2CO", "beta3CO"), name="betaCO") > sexCovCO <- mxMatrix(type="Full", nrow=1, ncol=nl, free=FALSE, labels="data.sex", name="sexCO") > FacMeans2CO <- mxAlgebra(expression= t(FacMeansCO+(betaCO*sexCO)), name="FacMeans2CO") > meansCO <- mxAlgebra(t(flCO %*% FacMeans2CO), name="meansCO") > expMeanCO <- mxAlgebra(cbind(meansCO, meansCO), name="expMeanCO") > > > covMZMN <- mxAlgebra( expression= rbind( cbind(AMN+CMN+EMN, AMN+CMN), + cbind( AMN+CMN, AMN+CMN+EMN)), name="expCovMZMN" ) > covDZMN <- mxAlgebra( expression= rbind( cbind(AMN+CMN+EMN, 0.5%x%AMN+CMN), + cbind(0.5%x%AMN+CMN, AMN+CMN+EMN)), name="expCovDZMN" ) > > covMZCO <- mxAlgebra( expression= rbind( cbind(ACO+CCO+ECO, ACO+CCO), + cbind( ACO+CCO, ACO+CCO+ECO)), name="expCovMZCO" ) > covDZCO <- mxAlgebra( expression= rbind( cbind(ACO+CCO+ECO, 0.5%x%ACO+CCO), + cbind(0.5%x%ACO+CCO, ACO+CCO+ECO)), name="expCovDZCO" ) > > > ### Combine Groups > objMZMN <- mxExpectationNormal( covariance="expCovMZMN", means="expMeanMN", dimnames=selVarsMN ) > objDZMN <- mxExpectationNormal( covariance="expCovDZMN", means="expMeanMN", dimnames=selVarsMN ) > > objMZCO <- mxExpectationNormal( covariance="expCovMZCO", means="expMeanCO", dimnames=selVarsCO ) > objDZCO <- mxExpectationNormal( covariance="expCovDZCO", means="expMeanCO", dimnames=selVarsCO ) > > parsMN <- list( XMN, YMN, ZMN, AlMN, ClMN, ElMN, pathAsMN, pathCsMN, pathEsMN, pathflMN, covAMN, covCMN, + covEMN, expMeanMN, regCoefMN, sexCovMN, pathFlGrowthMN, pathFlBDMN, covPMN, matIMN, invSDMN, + FacMeansMN, FacMeans2MN, meansMN, VarLMN, corLMN) > > parsCO <- list( XCO, YCO, ZCO, AlCO, ClCO, ElCO, pathAsCO, pathCsCO, pathEsCO, pathflCO, covACO, covCCO, + covECO, expMeanCO, regCoefCO, sexCovCO, pathFlGrowthCO, pathFlBDCO, covPCO, matICO, invSDCO, + FacMeansCO, FacMeans2CO, meansCO, VarLCO, corLCO) > > funML <- mxFitFunctionML() > > mtfsModelMZ <- mxModel( covMZMN, mtfsDataMZ, objMZMN, funML, name="MZMN", parsMN) > mtfsModelDZ <- mxModel( covDZMN, mtfsDataDZ, objDZMN, funML, name="DZMN", parsMN) > > caddModelMZ <- mxModel( covMZCO, caddDataMZ, objMZCO, funML, name="MZCO", parsCO) > caddModelDZ <- mxModel( covDZCO, caddDataDZ, objDZCO, funML, name="DZCO", parsCO) > #run model > fitML <- mxFitFunctionMultigroup(c("MZMN.fitfunction","DZMN.fitfunction", "MZCO.fitfunction","DZCO.fitfunction")) > LG <- mxModel("LinearGrowthACE", mtfsModelMZ, mtfsModelDZ, caddModelMZ, caddModelDZ, fitML) > mxOption(NULL,"Default optimizer","SLSQP") > mxOption(model= LG, key="Number of Threads", value= (omxDetectCores() - 1)) MxModel 'LinearGrowthACE' type : default $matrices : NULL $algebras : NULL $constraints : NULL $intervals : NULL $latentVars : none $manifestVars : none $data : NULL $submodels : 'MZMN', 'DZMN', 'MZCO', and 'DZCO' $expectation : NULL $fitfunction : MxFitFunctionMultigroup $compute : NULL $independent : FALSE $options : 'Number of Threads' = '23' $output : FALSE > LGFit <- mxRun(LG) ### Run model > RefModLG <- mxRefModels(LGFit, run=TRUE) > summary(LGFit, refModels=RefModLG) Summary of LinearGrowthACE free parameters: name matrix row col Estimate Std.Error A lbound ubound 1 x11MN MZMN.XMN 1 1 2.350231e-01 0.019665428 2 x21MN MZMN.XMN 2 1 4.476988e-01 0.049404162 3 x22MN MZMN.XMN 3 1 7.697653e-02 0.016325681 4 x31MN MZMN.XMN 2 2 2.459550e-01 0.045230706 5 x32MN MZMN.XMN 3 2 1.592254e-01 0.012289353 6 x33MN MZMN.XMN 3 3 -8.703119e-05 0.534074871 ! 7 y11MN MZMN.YMN 1 1 1.189033e-01 0.030692705 8 y21MN MZMN.YMN 2 1 3.907908e-01 0.054240560 9 y22MN MZMN.YMN 3 1 5.934731e-03 0.021086477 10 y31MN MZMN.YMN 2 2 -1.954162e-04 0.229984237 ! 11 y32MN MZMN.YMN 3 2 -1.030907e-05 0.033594477 12 y33MN MZMN.YMN 3 3 3.912847e-06 0.032006974 ! 13 z11MN MZMN.ZMN 1 1 -1.414033e-01 0.013703879 14 z21MN MZMN.ZMN 2 1 -1.622993e-01 0.030548302 15 z22MN MZMN.ZMN 3 1 -5.940402e-03 0.010998580 16 z31MN MZMN.ZMN 2 2 2.164102e-01 0.025476446 17 z32MN MZMN.ZMN 3 2 1.050715e-01 0.009676249 18 z33MN MZMN.ZMN 3 3 -7.730085e-06 0.049689885 ! 19 a11MN MZMN.asMN 1 1 1.000000e-05 0.115148493 ! 0! 20 a22MN MZMN.asMN 2 2 2.276161e-01 0.014195507 1e-05 21 a33MN MZMN.asMN 3 3 3.238890e-01 0.016990215 1e-05 22 a44MN MZMN.asMN 4 4 4.586061e-01 0.034692234 1e-05 23 a55MN MZMN.asMN 5 5 7.666176e-01 0.032189822 1e-05 24 a66MN MZMN.asMN 6 6 7.601724e-01 0.058784399 ! 1e-05 25 a77MN MZMN.asMN 7 7 6.238509e-01 0.214958692 ! 1e-05 26 c11MN MZMN.csMN 1 1 1.536932e-01 0.019499833 1e-05 27 c22MN MZMN.csMN 2 2 1.000000e-05 0.054805955 ! 0! 28 c33MN MZMN.csMN 3 3 1.000000e-05 0.080068093 ! 0! 29 c44MN MZMN.csMN 4 4 1.000000e-05 0.604540141 ! 0! 30 c55MN MZMN.csMN 5 5 4.961608e-05 0.310152093 ! 0! 31 c66MN MZMN.csMN 6 6 1.000000e-05 0.391927291 ! 0! 32 c77MN MZMN.csMN 7 7 1.000000e-05 0.225401233 ! 0! 33 e11MN MZMN.esMN 1 1 3.539216e-01 0.008992914 1e-05 34 e22MN MZMN.esMN 2 2 2.847221e-01 0.009069862 1e-05 35 e33MN MZMN.esMN 3 3 5.032282e-01 0.010178353 1e-05 36 e44MN MZMN.esMN 4 4 4.610736e-01 0.015285214 1e-05 37 e55MN MZMN.esMN 5 5 8.875363e-01 0.020789991 1e-05 38 e66MN MZMN.esMN 6 6 7.631839e-01 0.042206866 1e-05 39 e77MN MZMN.esMN 7 7 6.875564e-01 0.136930450 ! 1e-05 40 beta1MN MZMN.betaMN 1 1 2.182662e-03 0.016516195 41 beta2MN MZMN.betaMN 1 2 -1.575059e-01 0.039707114 42 beta3MN MZMN.betaMN 1 3 -8.573552e-02 0.012788117 43 fBD2MN MZMN.flBDMN 2 1 1.046883e+00 0.060069237 1e-06 44 fBD3MN MZMN.flBDMN 3 1 6.453401e-01 0.051645729 1e-06 45 mean1MN MZMN.FacMeansMN 1 1 1.407761e-03 0.026453943 46 mean2MN MZMN.FacMeansMN 1 2 6.557157e-01 0.063642868 47 mean3MN MZMN.FacMeansMN 1 3 2.454352e-01 0.020404539 48 x11CO MZCO.XCO 1 1 2.163413e-01 0.041563101 49 x21CO MZCO.XCO 2 1 4.189341e-01 0.076498364 50 x22CO MZCO.XCO 3 1 6.148707e-02 0.025123303 51 x31CO MZCO.XCO 2 2 -3.085483e-05 0.304159696 ! 52 x32CO MZCO.XCO 3 2 -8.155192e-06 0.095907159 53 x33CO MZCO.XCO 3 3 2.016162e-08 0.033452399 54 y11CO MZCO.YCO 1 1 2.114612e-01 0.038946958 55 y21CO MZCO.YCO 2 1 3.464137e-01 0.082569035 56 y22CO MZCO.YCO 3 1 3.544332e-02 0.024916935 57 y31CO MZCO.YCO 2 2 2.667176e-01 0.062601049 58 y32CO MZCO.YCO 3 2 8.645020e-02 0.019488048 59 y33CO MZCO.YCO 3 3 -9.569009e-06 0.038104004 ! 60 z11CO MZCO.ZCO 1 1 -1.621558e-01 0.017504139 61 z21CO MZCO.ZCO 2 1 -6.077832e-02 0.039972323 62 z22CO MZCO.ZCO 3 1 1.348051e-02 0.013937947 63 z31CO MZCO.ZCO 2 2 2.049667e-01 0.040222460 64 z32CO MZCO.ZCO 3 2 8.115915e-02 0.015811182 65 z33CO MZCO.ZCO 3 3 1.670152e-06 0.030837893 66 a11CO MZCO.asCO 1 1 1.868966e-01 0.071957080 1e-05 67 a22CO MZCO.asCO 2 2 2.468359e-01 0.055529349 1e-05 68 a33CO MZCO.asCO 3 3 3.347518e-01 0.023771534 1e-05 69 a44CO MZCO.asCO 4 4 5.169549e-01 0.115270819 1e-05 70 a55CO MZCO.asCO 5 5 8.255215e-01 0.068254041 1e-05 71 a66CO MZCO.asCO 6 6 1.000000e-05 0.624370165 ! 0! 72 c11CO MZCO.csCO 1 1 4.479302e-02 0.257786879 ! 0! 73 c22CO MZCO.csCO 2 2 1.378717e-01 0.088433212 1e-05 74 c33CO MZCO.csCO 3 3 1.000000e-05 0.150798502 ! 0! 75 c44CO MZCO.csCO 4 4 5.054821e-01 0.110462613 1e-05 76 c55CO MZCO.csCO 5 5 4.627966e-05 0.257589610 ! 0! 77 c66CO MZCO.csCO 6 6 7.083170e-01 0.164366439 ! 1e-05 78 e11CO MZCO.esCO 1 1 3.356970e-01 0.012746842 1e-05 79 e22CO MZCO.esCO 2 2 3.346053e-01 0.011965271 1e-05 80 e33CO MZCO.esCO 3 3 5.378941e-01 0.014031685 1e-05 81 e44CO MZCO.esCO 4 4 5.466994e-01 0.024825428 1e-05 82 e55CO MZCO.esCO 5 5 1.012656e+00 0.046479030 1e-05 83 e66CO MZCO.esCO 6 6 1.146675e+00 0.083891143 1e-05 84 beta1CO MZCO.betaCO 1 1 -1.051730e-02 0.022342006 85 beta2CO MZCO.betaCO 1 2 -2.182118e-01 0.058558478 86 beta3CO MZCO.betaCO 1 3 -1.092777e-01 0.015895981 87 fBD2CO MZCO.flBDCO 2 1 9.066037e-01 0.075418821 1e-06 88 fBD3CO MZCO.flBDCO 3 1 6.820435e-01 0.059948591 1e-06 89 mean1CO MZCO.FacMeansCO 1 1 1.560580e-02 0.035973597 90 mean2CO MZCO.FacMeansCO 1 2 8.289390e-01 0.094090533 91 mean3CO MZCO.FacMeansCO 1 3 2.690860e-01 0.025747922 Model Statistics: | Parameters | Degrees of Freedom | Fit (-2lnL units) Model: 91 28874 65337.71 Saturated: 418 28547 64480.15 Independence: 104 28861 71648.89 Number of observations/statistics: 3257/28965 chi-square: χ² ( df=327 ) = 857.5611, p = 3.400343e-49 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 7589.714 65519.71 65525.00 BIC: -168211.420 66073.77 65784.63 CFI: 0.9225993 TLI: 0.9256764 (also known as NNFI) RMSEA: 0.02231953 [95% CI (0.02014759, 0.024493)] Prob(RMSEA <= 0.05): 1 timestamp: 2018-08-19 23:59:37 Wall clock time: 1048.378 secs optimizer: SLSQP OpenMx version number: 2.9.6 Need help? See help(mxSummary) > > > > ######################################################################## > #correlation constrained models > ######################################################################## > > nl <- 3 # number of latent factors > # Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s) > XMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, 0), + labels=c("x11MN", + "x21MN", "x22MN", + "x31MN", "x32MN", "x33MN"), name="XMN" ) > > YMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("y11MN", + "y21MN", "y22MN", + "y31MN", "y32MN", "y33MN"), name="YMN" ) > > ZMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("z11MN", + "z21MN", "z22MN", + "z31MN", "z32MN", "z33MN"), name="ZMN" ) > > XCO <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, 0), + labels=c("x11CO", + "x21CO", "x22CO", + "x31CO", "x32CO", "x33CO"), name="XCO" ) > > YCO <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("y11CO", + "y21CO", "y22CO", + "y31CO", "y32CO", "y33CO"), name="YCO" ) > > ZCO <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, + values=c(.6, + .1, .1, + .1, .1, .1), + labels=c("z11CO", + "z21CO", "z22CO", + "z31CO", "z32CO", "z33CO"), name="ZCO" ) > > > AlMN <- mxAlgebra(XMN %*% t(XMN), name="AlMN") > ClMN <- mxAlgebra(YMN %*% t(YMN), name="ClMN") > ElMN <- mxAlgebra(ZMN %*% t(ZMN), name="ElMN") > > AlCO <- mxAlgebra(XCO %*% t(XCO), name="AlCO") > ClCO <- mxAlgebra(YCO %*% t(YCO), name="ClCO") > ElCO <- mxAlgebra(ZCO %*% t(ZCO), name="ElCO") > > VarLMN <- mxAlgebra(expression = AlMN + ClMN + ElMN, name="VarLMN") > VarLCO <- mxAlgebra(expression = AlCO + ClCO + ElCO, name="VarLCO") > > corLMN <- mxAlgebra(cov2cor(VarLMN), name="corLMN") > corLCO <- mxAlgebra(cov2cor(VarLCO), name="corLCO") > > CO_ACons2 <- mxMatrix("Full", 1, 1, free=FALSE, values=0, labels="AlCO[2,2]", name="CO_ACons2") > > # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors > pathAsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3, + labels=c("a11MN", "a22MN", "a33MN", "a44MN", "a55MN", "a66MN", "a77MN"), + lbound=.00001, name="asMN" ) > pathCsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3, + labels=c("c11MN","c22MN", "c33MN", "c44MN", "c55MN", "c66MN", "c77MN"), + lbound=.00001, name="csMN" ) > pathEsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3, + labels=c("e11MN", "e22MN", "e33MN", "e44MN", "e55MN", "e66MN", "e77MN"), + lbound=.00001, name="esMN" ) > > pathAsCO <- mxMatrix(type="Diag", nrow=nvCO, ncol=nvCO, free=TRUE, values=.3, + labels=c("a11CO", "a22CO", "a33CO", "a44CO", "a55CO", "a66CO"), + lbound=.00001, name="asCO" ) > pathCsCO <- mxMatrix(type="Diag", nrow=nvCO, ncol=nvCO, free=TRUE, values=.3, + labels=c("c11CO","c22CO", "c33CO", "c44CO", "c55CO", "c66CO"), + lbound=.00001, name="csCO" ) > pathEsCO <- mxMatrix(type="Diag", nrow=nvCO, ncol=nvCO, free=TRUE, values=.3, + labels=c("e11CO", "e22CO", "e33CO", "e44CO", "e55CO", "e66CO"), + lbound=.00001, name="esCO" ) > > > > > pathFlGrowthMN <- mxMatrix(type="Full", nrow=4, ncol=nl, free=FALSE, + labels=c("data.zero", "data.intercept", "data.age1twinpairC", + "data.zero", "data.intercept", "data.age2twinpairC", + "data.zero", "data.intercept", "data.age3twinpairC", + "data.zero", "data.intercept", "data.age4twinpairC"), name="flGrowthMN", byrow=T) > > pathFlGrowthCO <- mxMatrix(type="Full", nrow=3, ncol=nl, free=FALSE, + labels=c("data.zero", "data.intercept", "data.age1twinpairC", + "data.zero", "data.intercept", "data.age2twinpairC", + "data.zero", "data.intercept", "data.age3twinpairC"), name="flGrowthCO", byrow=T) > > > flFree <- c(FALSE, FALSE, FALSE, + TRUE, FALSE, FALSE, + TRUE, FALSE, FALSE) > > flValues <- c(1, 0, 0, + .6, 0, 0, + .6, 0, 0) > > #set first loading to 1 and make it fixed for BD to manifest variables > pathFlBDMN <- mxMatrix(type="Full", nrow=3, ncol=nl, free=flFree, + values=flValues, lbound=.000001, labels=c("fBD1MN", NA, NA, + "fBD2MN", NA, NA, + "fBD3MN", NA, NA), name="flBDMN", byrow=T) > > pathFlBDCO <- mxMatrix(type="Full", nrow=3, ncol=nl, free=flFree, + values=flValues, lbound=.000001, labels=c("fBD1CO", NA, NA, + "fBD2CO", NA, NA, + "fBD3CO", NA, NA), name="flBDCO", byrow=T) > > #combine factor loading matrices > pathflMN <- mxAlgebra(rbind(flBDMN, flGrowthMN), name="flMN") > > pathflCO <- mxAlgebra(rbind(flBDCO, flGrowthCO), name="flCO") > > covAMN <- mxAlgebra( expression=flMN %*% AlMN %*% t(flMN) + asMN %*% t(asMN), name="AMN" ) > covCMN <- mxAlgebra( expression=flMN %*% ClMN %*% t(flMN) + csMN %*% t(csMN), name="CMN" ) > covEMN <- mxAlgebra( expression=flMN %*% ElMN %*% t(flMN) + esMN %*% t(esMN), name="EMN" ) > > covACO <- mxAlgebra( expression=flCO %*% AlCO %*% t(flCO) + asCO %*% t(asCO), name="ACO" ) > covCCO <- mxAlgebra( expression=flCO %*% ClCO %*% t(flCO) + csCO %*% t(csCO), name="CCO" ) > covECO <- mxAlgebra( expression=flCO %*% ElCO %*% t(flCO) + esCO %*% t(esCO), name="ECO" ) > > matIMN <- mxMatrix( type="Iden", nrow=3, ncol=3, name="IMN") > covPMN <- mxAlgebra( expression=AMN[1:3, 1:3]+CMN[1:3, 1:3]+EMN[1:3, 1:3], name="VMN" ) > invSDMN <- mxAlgebra( expression=solve(sqrt(IMN*VMN)), name="iSDMN") > > matICO <- mxMatrix( type="Iden", nrow=3, ncol=3, name="ICO") > covPCO <- mxAlgebra( expression=ACO[1:3, 1:3]+CCO[1:3, 1:3]+ECO[1:3, 1:3], name="VCO" ) > invSDCO <- mxAlgebra( expression=solve(sqrt(ICO*VCO)), name="iSDCO") > > #estimate means > #estimate means for the intercept and slope > FacMeansMN <- mxMatrix("Full",nrow=1,ncol=nl, free=c(TRUE, TRUE, TRUE), values=c(0,.8,.3), labels=c("mean1MN", "mean2MN", "mean3MN"), name="FacMeansMN") > regCoefMN <- mxMatrix(type="Full", nrow=1, ncol=nl, free=TRUE, values=0, + labels=c("beta1MN", "beta2MN", "beta3MN"), name="betaMN") > sexCovMN <- mxMatrix(type="Full", nrow=1, ncol=nl, free=FALSE, labels="data.sex", name="sexMN") > FacMeans2MN <- mxAlgebra(expression= t(FacMeansMN+(betaMN*sexMN)), name="FacMeans2MN") > meansMN <- mxAlgebra(t(flMN %*% FacMeans2MN), name="meansMN") > expMeanMN <- mxAlgebra(cbind(meansMN, meansMN), name="expMeanMN") > > > FacMeansCO <- mxMatrix("Full",nrow=1,ncol=nl, free=c(TRUE, TRUE, TRUE), values=c(0,.8,.3), labels=c("mean1CO", "mean2CO", "mean3CO"), name="FacMeansCO") > regCoefCO <- mxMatrix(type="Full", nrow=1, ncol=nl, free=TRUE, values=0, + labels=c("beta1CO", "beta2CO", "beta3CO"), name="betaCO") > sexCovCO <- mxMatrix(type="Full", nrow=1, ncol=nl, free=FALSE, labels="data.sex", name="sexCO") > FacMeans2CO <- mxAlgebra(expression= t(FacMeansCO+(betaCO*sexCO)), name="FacMeans2CO") > meansCO <- mxAlgebra(t(flCO %*% FacMeans2CO), name="meansCO") > expMeanCO <- mxAlgebra(cbind(meansCO, meansCO), name="expMeanCO") > > > covMZMN <- mxAlgebra( expression= rbind( cbind(AMN+CMN+EMN, AMN+CMN), + cbind( AMN+CMN, AMN+CMN+EMN)), name="expCovMZMN" ) > covDZMN <- mxAlgebra( expression= rbind( cbind(AMN+CMN+EMN, 0.5%x%AMN+CMN), + cbind(0.5%x%AMN+CMN, AMN+CMN+EMN)), name="expCovDZMN" ) > > covMZCO <- mxAlgebra( expression= rbind( cbind(ACO+CCO+ECO, ACO+CCO), + cbind( ACO+CCO, ACO+CCO+ECO)), name="expCovMZCO" ) > covDZCO <- mxAlgebra( expression= rbind( cbind(ACO+CCO+ECO, 0.5%x%ACO+CCO), + cbind(0.5%x%ACO+CCO, ACO+CCO+ECO)), name="expCovDZCO" ) > > > ### Combine Groups > objMZMN <- mxExpectationNormal( covariance="expCovMZMN", means="expMeanMN", dimnames=selVarsMN ) > objDZMN <- mxExpectationNormal( covariance="expCovDZMN", means="expMeanMN", dimnames=selVarsMN ) > > objMZCO <- mxExpectationNormal( covariance="expCovMZCO", means="expMeanCO", dimnames=selVarsCO ) > objDZCO <- mxExpectationNormal( covariance="expCovDZCO", means="expMeanCO", dimnames=selVarsCO ) > > parsMN_C <- list( XMN, YMN, ZMN, AlMN, ClMN, ElMN, pathAsMN, pathCsMN, pathEsMN, pathflMN, covAMN, covCMN, + covEMN, expMeanMN, regCoefMN, sexCovMN, pathFlGrowthMN, pathFlBDMN, covPMN, matIMN, invSDMN, + FacMeansMN, FacMeans2MN, meansMN, VarLMN, corLMN) > > parsCO_C <- list( XCO, YCO, ZCO, AlCO, ClCO, ElCO, pathAsCO, pathCsCO, pathEsCO, pathflCO, covACO, covCCO, + covECO, expMeanCO, regCoefCO, sexCovCO, pathFlGrowthCO, pathFlBDCO, covPCO, matICO, invSDCO, + FacMeansCO, FacMeans2CO, meansCO, VarLCO, corLCO, CO_ACons2) > > funML <- mxFitFunctionML() > > mtfsModelMZ <- mxModel( covMZMN, mtfsDataMZ, objMZMN, funML, name="MZMN", parsMN_C) > mtfsModelDZ <- mxModel( covDZMN, mtfsDataDZ, objDZMN, funML, name="DZMN", parsMN_C) > > caddModelMZ <- mxModel( covMZCO, caddDataMZ, objMZCO, funML, name="MZCO", parsCO_C) > caddModelDZ <- mxModel( covDZCO, caddDataDZ, objDZCO, funML, name="DZCO", parsCO_C) > #run model > fitML <- mxFitFunctionMultigroup(c("MZMN.fitfunction","DZMN.fitfunction", "MZCO.fitfunction","DZCO.fitfunction")) > LG_Con <- mxModel("LinearGrowthACE", mtfsModelMZ, mtfsModelDZ, caddModelMZ, caddModelDZ, fitML) > mxOption(NULL,"Default optimizer","SLSQP") > mxOption(model= LG_Con, key="Number of Threads", value= (omxDetectCores() - 1)) MxModel 'LinearGrowthACE' type : default $matrices : NULL $algebras : NULL $constraints : NULL $intervals : NULL $latentVars : none $manifestVars : none $data : NULL $submodels : 'MZMN', 'DZMN', 'MZCO', and 'DZCO' $expectation : NULL $fitfunction : MxFitFunctionMultigroup $compute : NULL $independent : FALSE $options : 'Number of Threads' = '23' $output : FALSE > LGFit_Con <- mxRun(LG_Con) ### Run model > RefModLG_Con <- mxRefModels(LGFit_Con, run=TRUE) > summary(LGFit_Con, refModels=RefModLG_Con) Summary of LinearGrowthACE free parameters: name matrix row col Estimate Std.Error A lbound ubound 1 x11MN MZMN.XMN 1 1 2.350210e-01 0.019665682 2 x21MN MZMN.XMN 2 1 4.476925e-01 0.049401871 3 x22MN MZMN.XMN 3 1 7.697643e-02 0.016326450 4 x31MN MZMN.XMN 2 2 2.459597e-01 0.045234969 5 x32MN MZMN.XMN 3 2 1.592253e-01 0.012297308 6 x33MN MZMN.XMN 3 3 -1.887677e-04 0.545491743 ! 7 y11MN MZMN.YMN 1 1 1.189077e-01 0.030689172 8 y21MN MZMN.YMN 2 1 3.907960e-01 0.054236445 9 y22MN MZMN.YMN 3 1 5.935173e-03 0.021087719 10 y31MN MZMN.YMN 2 2 5.907868e-05 0.230514745 ! 11 y32MN MZMN.YMN 3 2 7.731839e-06 0.033595739 ! 12 y33MN MZMN.YMN 3 3 4.024638e-06 0.032006153 ! 13 z11MN MZMN.ZMN 1 1 -1.414029e-01 0.013703807 14 z21MN MZMN.ZMN 2 1 -1.623016e-01 0.030548170 15 z22MN MZMN.ZMN 3 1 -5.941509e-03 0.010997539 16 z31MN MZMN.ZMN 2 2 2.164095e-01 0.025477161 17 z32MN MZMN.ZMN 3 2 1.050712e-01 0.009667683 18 z33MN MZMN.ZMN 3 3 -2.899721e-05 0.049702257 ! 19 a11MN MZMN.asMN 1 1 3.554837e-05 0.115093564 ! 0! 20 a22MN MZMN.asMN 2 2 2.276158e-01 0.014194924 1e-05 21 a33MN MZMN.asMN 3 3 3.238884e-01 0.016990249 1e-05 22 a44MN MZMN.asMN 4 4 4.586061e-01 0.034692385 ! 1e-05 23 a55MN MZMN.asMN 5 5 7.666174e-01 0.032189767 1e-05 24 a66MN MZMN.asMN 6 6 7.601754e-01 0.058770245 1e-05 25 a77MN MZMN.asMN 7 7 6.238699e-01 0.214088736 ! 1e-05 26 c11MN MZMN.csMN 1 1 1.536918e-01 0.019499370 1e-05 27 c22MN MZMN.csMN 2 2 1.000000e-05 0.054808445 ! 0! 28 c33MN MZMN.csMN 3 3 1.000000e-05 0.080051039 ! 0! 29 c44MN MZMN.csMN 4 4 1.000000e-05 0.597298805 ! 0! 30 c55MN MZMN.csMN 5 5 1.000000e-05 0.310799579 ! 0! 31 c66MN MZMN.csMN 6 6 1.000000e-05 0.390330992 ! 0! 32 c77MN MZMN.csMN 7 7 1.000000e-05 0.225028938 ! 0! 33 e11MN MZMN.esMN 1 1 3.539219e-01 0.008992756 1e-05 34 e22MN MZMN.esMN 2 2 2.847225e-01 0.009069645 1e-05 35 e33MN MZMN.esMN 3 3 5.032281e-01 0.010178326 1e-05 36 e44MN MZMN.esMN 4 4 4.610728e-01 0.015284901 1e-05 37 e55MN MZMN.esMN 5 5 8.875360e-01 0.020790054 1e-05 38 e66MN MZMN.esMN 6 6 7.631799e-01 0.042182759 ! 1e-05 39 e77MN MZMN.esMN 7 7 6.875365e-01 0.136465965 ! 1e-05 40 beta1MN MZMN.betaMN 1 1 2.181751e-03 0.016513574 41 beta2MN MZMN.betaMN 1 2 -1.575067e-01 0.039693841 42 beta3MN MZMN.betaMN 1 3 -8.573576e-02 0.012787403 43 fBD2MN MZMN.flBDMN 2 1 1.046881e+00 0.060054278 1e-06 44 fBD3MN MZMN.flBDMN 3 1 6.453401e-01 0.051644595 1e-06 45 mean1MN MZMN.FacMeansMN 1 1 1.408951e-03 0.026449748 46 mean2MN MZMN.FacMeansMN 1 2 6.557163e-01 0.063621512 47 mean3MN MZMN.FacMeansMN 1 3 2.454356e-01 0.020403399 48 x11CO MZCO.XCO 1 1 2.163309e-01 0.041570791 49 x21CO MZCO.XCO 2 1 4.189580e-01 0.076480119 50 x22CO MZCO.XCO 3 1 6.149654e-02 0.025123985 51 x31CO MZCO.XCO 2 2 4.223360e-06 0.304349924 ! 52 x32CO MZCO.XCO 3 2 3.520017e-06 0.095964242 53 x33CO MZCO.XCO 3 3 3.538853e-06 0.033451296 ! 54 y11CO MZCO.YCO 1 1 2.114679e-01 0.038944234 55 y21CO MZCO.YCO 2 1 3.463907e-01 0.082550059 56 y22CO MZCO.YCO 3 1 3.543392e-02 0.024913916 57 y31CO MZCO.YCO 2 2 2.667154e-01 0.062601804 58 y32CO MZCO.YCO 3 2 8.644977e-02 0.019495414 59 y33CO MZCO.YCO 3 3 3.694467e-06 0.038105596 ! 60 z11CO MZCO.ZCO 1 1 -1.621564e-01 0.017502642 61 z21CO MZCO.ZCO 2 1 -6.077706e-02 0.039972634 62 z22CO MZCO.ZCO 3 1 1.348076e-02 0.013937975 63 z31CO MZCO.ZCO 2 2 2.049606e-01 0.040229671 64 z32CO MZCO.ZCO 3 2 8.115663e-02 0.015816398 65 z33CO MZCO.ZCO 3 3 3.673791e-06 0.030837436 ! 66 a11CO MZCO.asCO 1 1 1.869567e-01 0.072174752 1e-05 67 a22CO MZCO.asCO 2 2 2.468433e-01 0.055520377 1e-05 68 a33CO MZCO.asCO 3 3 3.347534e-01 0.023771021 1e-05 69 a44CO MZCO.asCO 4 4 5.169486e-01 0.115110314 1e-05 70 a55CO MZCO.asCO 5 5 8.255102e-01 0.068277240 1e-05 71 a66CO MZCO.asCO 6 6 1.974372e-05 0.624448577 ! 0! 72 c11CO MZCO.csCO 1 1 4.456158e-02 0.260124472 ! 0! 73 c22CO MZCO.csCO 2 2 1.378579e-01 0.088452315 1e-05 74 c33CO MZCO.csCO 3 3 1.116386e-05 0.150691905 ! 0! 75 c44CO MZCO.csCO 4 4 5.054873e-01 0.110290655 1e-05 76 c55CO MZCO.csCO 5 5 1.000000e-05 0.257526262 ! 0! 77 c66CO MZCO.csCO 6 6 7.083022e-01 0.164494493 1e-05 78 e11CO MZCO.esCO 1 1 3.356941e-01 0.012749793 1e-05 79 e22CO MZCO.esCO 2 2 3.346045e-01 0.011965760 1e-05 80 e33CO MZCO.esCO 3 3 5.378944e-01 0.014031373 1e-05 81 e44CO MZCO.esCO 4 4 5.466991e-01 0.024824319 1e-05 82 e55CO MZCO.esCO 5 5 1.012660e+00 0.046502080 1e-05 83 e66CO MZCO.esCO 6 6 1.146683e+00 0.083950104 1e-05 84 beta1CO MZCO.betaCO 1 1 -1.052049e-02 0.022349846 85 beta2CO MZCO.betaCO 1 2 -2.182137e-01 0.058629040 86 beta3CO MZCO.betaCO 1 3 -1.092777e-01 0.015896896 87 fBD2CO MZCO.flBDCO 2 1 9.066088e-01 0.075440682 1e-06 88 fBD3CO MZCO.flBDCO 3 1 6.820440e-01 0.059952083 1e-06 89 mean1CO MZCO.FacMeansCO 1 1 1.561032e-02 0.035986050 90 mean2CO MZCO.FacMeansCO 1 2 8.289420e-01 0.094207417 91 mean3CO MZCO.FacMeansCO 1 3 2.690863e-01 0.025749418 Model Statistics: | Parameters | Degrees of Freedom | Fit (-2lnL units) Model: 91 28874 65337.71 Saturated: 418 28547 64480.15 Independence: 104 28861 71648.89 Number of observations/statistics: 3257/28965 chi-square: χ² ( df=327 ) = 857.5611, p = 3.400344e-49 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 7589.714 65519.71 65525.00 BIC: -168211.420 66073.77 65784.63 CFI: 0.9225993 TLI: 0.9256764 (also known as NNFI) RMSEA: 0.02231953 [95% CI (0.02014759, 0.024493)] Prob(RMSEA <= 0.05): 1 timestamp: 2018-08-20 00:28:48 Wall clock time: 1078.656 secs optimizer: SLSQP OpenMx version number: 2.9.6 Need help? See help(mxSummary) > > mxCompare(LGFit, LGFit_Con) base comparison ep minus2LL df AIC diffLL 1 LinearGrowthACE 91 65337.71 28874 7589.714 NA 2 LinearGrowthACE LinearGrowthACE 91 65337.71 28874 7589.714 -1.298889e-06 diffdf p 1 NA NA 2 0 NA > > > LGFit$output$algebras$MZCO.AlCO [,1] [,2] [,3] [1,] 0.04680354 0.09063273 0.01330219 [2,] 0.09063273 0.17550576 0.02575903 [3,] 0.01330219 0.02575903 0.00378066 > LGFit_Con$output$algebras$MZCO.AlCO [,1] [,2] [,3] [1,] 0.04679905 0.09063356 0.013303601 [2,] 0.09063356 0.17552583 0.025764470 [3,] 0.01330360 0.02576447 0.003781825 > > save.image() >