table(data$depvar_1) # 3 Thresholds table(data$depvar_2) # 3 Thresholds data$depvar_1 <- mxFactor(data$depvar_1, levels = c(1,2,3,4)) data$depvar_2 <- mxFactor(data$depvar_2, levels = c(1,2,3,4)) data <- as.data.frame(data) # some short-cuts nth <- 3 vars <- 'depvar_' nv <- 1 ntv <- nv*2 selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") modVars <- "mod1" # 6 groups: # 1 = MZ High 22 # 2 = MZ Low 11 # 3 = MZ High/Low 21 # 4 = DZ High 22 # 5 = DZ Low 11 # 6 = DZ High/Low 21 data <- data %>% mutate(zyg = ifelse(zyg == 1 & mod2_1 == 2 & mod2_2 == 2,1, ifelse(zyg == 1 & mod2_1 == 1 & mod2_2 == 1,2, ifelse(zyg == 1 & (mod2_1 == 2 & mod2_2 == 1),3, ifelse(zyg == 2 & mod2_1 == 2 & mod2_2 == 2,4, ifelse(zyg == 2 & mod2_1 == 1 & mod2_2 == 1,5, ifelse(zyg == 2 & (mod2_1 == 2 & mod2_2 == 1),6, NA))))))) %>% select(depvar_1, depvar_2, mod1, zyg, mod2_1, mod2_2) #View(data) table(data$zyg) mzhData <- subset(data, zyg==1, c(selVars, modVars)) mzlData <- subset(data, zyg==2, c(selVars, modVars)) mzhlData <- subset(data, zyg==3, c(selVars, modVars)) dzhData <- subset(data, zyg==4, c(selVars, modVars)) dzlData <- subset(data, zyg==5, c(selVars, modVars)) dzhlData <- subset(data, zyg==6, c(selVars, modVars)) dim(mzhData) dim(mzlData) dim(mzhlData) dim(dzhData) dim(dzlData) dim(dzhlData) svITh <- 1 # start value for increments frTh <- matrix(rep(c(F,F,(rep(T,nth-2)))),nrow=nth,ncol=nv) # free thresholds svTh <- matrix(rep(c(0,1,(rep(svITh,nth-2)))),nrow=nth,ncol=nv) # start values thresholds lbTh <- matrix(rep(rep(0.001,nth),nv),nrow=nth,ncol=nv) # lower bounds for thresholds svPa <- .2 # start value for path coefficient labMeh <- c("mh","mh") # labels for means for h twins labMel <- c("ml","ml") # labels for means for l twins labMehl <- c("mh","ml") # labels for means for hl twins labThh <- c(paste("t",1:nth,"h",sep=""),paste("t",1:nth,"h",sep="")) # labels for means for H twins labThl <- c(paste("t",1:nth,"l",sep=""),paste("t",1:nth,"l",sep="")) # labels for means for L twins labThhl <- c(paste("t",1:nth,"h",sep=""),paste("t",1:nth,"l",sep="")) # labels for means for HL twins # Matrices for moderator (definition var) defM <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.mod1"), name="defM" ) # Main Effect pathBh <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= c(0.5), label=c("b_modh"), name="bh" ) pathBl <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= c(0.5), label=c("b_modl"), name="bl" ) # conditional mean (Mu|M=0) meanh <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values = 0, label=labMeh, name="Meanh" ) meanl <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values = 0, label=labMel, name="Meanl" ) meanhl <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values = 0, label=labMehl, name="Meanhl" ) # expected means expMeanh <- mxAlgebra( expression= Meanh + cbind(bh%*%defM,bh%*%defM), name="expMeanh") expMeanl <- mxAlgebra( expression= Meanl + cbind(bl%*%defM,bl%*%defM), name="expMeanl") expMeanhl <- mxAlgebra( expression= Meanhl + cbind(bh%*%defM,bl%*%defM), name="expMeanhl") # Thresholds thinGh <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=frTh, values=svTh, lbound=lbTh, labels=labThh, name="thinGh" ) thinGl <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=frTh, values=svTh, lbound=lbTh, labels=labThl, name="thinGl" ) thinGhl <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=frTh, values=svTh, lbound=lbTh, labels=labThhl, name="thinGhl" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) thruGh <- mxAlgebra( expression= inc %*% thinGh, name="thruGh" ) # multiplication -> Th3>Th2>Th1 thruGl <- mxAlgebra( expression= inc %*% thinGl, name="thruGl" ) thruGhl <- mxAlgebra( expression= inc %*% thinGhl, name="thruGhl" ) # unmoderated coefficients pathAh <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="ah11", name="ah" ) pathCh <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="ch11", name="ch" ) pathEh <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="eh11", name="eh" ) pathAl <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="al11", name="al" ) pathCl <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="cl11", name="cl" ) pathEl <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="el11", name="el" ) # moderated coefficients pathAIh <- mxMatrix( type='Full', nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="aIh11", name="aIh" ) pathCIh <- mxMatrix( type='Full', nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="cIh11", name="cIh" ) pathEIh <- mxMatrix( type='Full', nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="eIh11", name="eIh" ) pathAIl <- mxMatrix( type='Full', nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="aIl11", name="aIl" ) pathCIl <- mxMatrix( type='Full', nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="cIl11", name="cIl" ) pathEIl <- mxMatrix( type='Full', nrow=nv, ncol=nv, free=TRUE, values = 0.5, label="eIl11", name="eIl" ) # moderated var comps covAIh <- mxAlgebra( expression=(ah+ defM%*%aIh) %*% t(ah+ defM%*%aIh), name="AIh" ) covCIh <- mxAlgebra( expression=(ch+ defM%*%cIh) %*% t(ch+ defM%*%cIh), name="CIh" ) covEIh <- mxAlgebra( expression=(eh+ defM%*%eIh) %*% t(eh+ defM%*%eIh), name="EIh" ) covAIl <- mxAlgebra( expression=(al+ defM%*%aIl) %*% t(al+ defM%*%aIl), name="AIl" ) covCIl <- mxAlgebra( expression=(cl+ defM%*%cIl) %*% t(cl+ defM%*%cIl), name="CIl" ) covEIl <- mxAlgebra( expression=(el+ defM%*%eIl) %*% t(el+ defM%*%eIl), name="EIl" ) # variance covPIh <- mxAlgebra( expression=AIh+CIh+EIh, name="VIh" ) covPIl <- mxAlgebra( expression=AIl+CIl+EIl, name="VIl" ) # covariances mz covMZh <- mxAlgebra( expression= AIh+CIh, name="cMZh" ) covMZl <- mxAlgebra( expression= AIl+CIl, name="cMZl" ) covMZhl <- mxAlgebra( expression= ((ah+ defM%*%aIh) %*% t(al+ defM%*%aIl))+ ((ch+ defM%*%cIh) %*% t(cl+ defM%*%cIl)), name="covMZhl" ) covMZlh <- mxAlgebra( expression= ((al+ defM%*%aIl) %*% t(ah+ defM%*%aIh))+ ((cl+ defM%*%cIl) %*% t(ch+ defM%*%cIh)), name="covMZlh" ) # covariances dz covDZh <- mxAlgebra( expression= 0.5%x%AIh+ CIh, name="cDZh" ) covDZl <- mxAlgebra( expression= 0.5%x%AIl+ CIl, name="cDZl" ) covDZhl <- mxAlgebra( expression= 0.5%x%((ah+ defM%*%aIh) %*% t(al+ defM%*%aIl))+ ((ch+ defM%*%cIh) %*% t(cl+ defM%*%cIl)), name="covDZhl" ) covDZlh <- mxAlgebra( expression= 0.5%x%((al+ defM%*%aIl) %*% t(ah+ defM%*%aIh))+ ((cl+ defM%*%cIl) %*% t(ch+ defM%*%cIh)), name="covDZlh" ) # expected covariance matrices expCovMZh <- mxAlgebra( expression= rbind( cbind(VIh, cMZh), cbind(t(cMZh), VIh)), name="expCovMZh" ) expCovMZl <- mxAlgebra( expression= rbind( cbind(VIl, cMZl), cbind(t(cMZl), VIl)), name="expCovMZl" ) expCovMZhl <- mxAlgebra( expression= rbind( cbind(VIh, covMZhl), cbind(covMZlh, VIl)), name="expCovMZhl" ) expCovDZh <- mxAlgebra( expression= rbind( cbind(VIh, cDZh), cbind(t(cDZh), VIh)), name="expCovDZh" ) expCovDZl <- mxAlgebra( expression= rbind( cbind(VIl, cDZl), cbind(t(cDZl), VIl)), name="expCovDZl" ) expCovDZhl <- mxAlgebra( expression= rbind( cbind(VIh, covDZhl), cbind(covDZlh, VIl)), name="expCovDZhl" ) # data dataMZh <- mxData( observed=mzhData, type="raw" ) dataMZl <- mxData( observed=mzlData, type="raw" ) dataMZhl <- mxData( observed=mzhlData, type="raw" ) dataDZh <- mxData( observed=dzhData, type="raw" ) dataDZl <- mxData( observed=dzlData, type="raw" ) dataDZhl <- mxData( observed=dzhlData, type="raw" ) # expectation objects and fit function expMZh <- mxExpectationNormal( covariance="expCovMZh", means="expMeanh", dimnames=selVars, thresholds="thruGh" ) expMZl <- mxExpectationNormal( covariance="expCovMZl", means="expMeanl", dimnames=selVars, thresholds="thruGl" ) expMZhl <- mxExpectationNormal( covariance="expCovMZhl", means="expMeanhl", dimnames=selVars, thresholds="thruGhl" ) expDZh <- mxExpectationNormal( covariance="expCovDZh", means="expMeanh", dimnames=selVars, thresholds="thruGh" ) expDZl <- mxExpectationNormal( covariance="expCovDZl", means="expMeanl", dimnames=selVars, thresholds="thruGl" ) expDZhl <- mxExpectationNormal( covariance="expCovDZhl", means="expMeanhl", dimnames=selVars, thresholds="thruGhl" ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups parsh <- list(pathBh,meanh,thinGh,inc,thruGh,pathAh, pathCh, pathEh,pathAIh, pathCIh, pathEIh, covAIh,covCIh,covEIh,covPIh) parsl <- list(pathBl,meanl,thinGl,inc,thruGl,pathAl, pathCl, pathEl,pathAIl, pathCIl, pathEIl, covAIl,covCIl,covEIl,covPIl) parshl <- list(thinGhl, inc,thruGhl,meanhl) modelMZh <- mxModel( defM,parsh, covMZh, expMeanh, expCovMZh, dataMZh, expMZh, funML, name="MZh" ) modelMZl <- mxModel( defM,parsl, covMZl, expMeanl, expCovMZl, dataMZl, expMZl, funML, name="MZl" ) modelMZhl <- mxModel( defM,parsh, parsl, parshl, covMZhl,covMZlh, expMeanhl, expCovMZhl, dataMZhl, expMZhl, funML, name="MZhl" ) modelDZh <- mxModel( defM,parsh, covDZh, expMeanh, expCovDZh, dataDZh, expDZh, funML, name="DZh" ) modelDZl <- mxModel( defM,parsl, covDZl, expMeanl, expCovDZl, dataDZl, expDZl, funML, name="DZl" ) modelDZhl <- mxModel( defM,parsh,parsl,parshl, covDZhl,covDZlh, expMeanhl, expCovDZhl, dataDZhl, expDZhl, funML, name="DZhl" ) multi <- mxFitFunctionMultigroup( c("MZh","MZl","MZhl","DZh","DZl","DZhl") ) modelACE <- mxModel( "modACE1", defM,parsh, parsl, parshl, modelMZh,modelMZl,modelMZhl, modelDZh,modelDZl,modelDZhl, dataMZh, dataMZl, dataMZhl, dataDZh,dataDZl,dataDZhl, multi) set.seed(1) mxOption(NULL,"Default optimizer","SLSQP") mod1_fit <- mxTryHardOrdinal( modelACE, extraTries = 10, intervals = TRUE, bestInitsOutput = TRUE, silent=FALSE, exhaustive = T)