# Binary

table(data$depvar_1) # 1 Threshold
table(data$depvar_2) # 1 Threshold
data$depvar_1 <- mxFactor(data$depvar_1, levels = c(1,2))
data$depvar_2 <- mxFactor(data$depvar_2, levels = c(1,2))
data <- as.data.frame(data)

nth       <- 1                       
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 00
  # 2 = MZ Low 11
  # 3 = MZ High/Low 01
  # 4 = DZ High 00 
  # 5 = DZ Low 11
  # 6 = DZ High/Low 01

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)


svPa      <- .2                        # start value for path coefficient

labThh   <- c(paste("t",1:nth,"h",sep=""),paste("t",1:nth,"h",sep=""))  # labels for means for FF twins
labThl   <- c(paste("t",1:nth,"l",sep=""),paste("t",1:nth,"l",sep=""))  # labels for means for MM twins
labThhl   <- c(paste("t",1:nth,"h",sep=""),paste("t",1:nth,"l",sep=""))  # labels for means for FM twins



# Matrices for moderator (definition var 
defM      <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.mod1"), name="defM" )

# Main Effects 
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)) 
mean     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values = 0, 
                       label=c("mean_1","mean_2"), name="Mean" )


# expected meanselwert
expMeanh     <- mxAlgebra( expression= Mean + cbind(bh%*%defM,bh%*%defM), name="expMeanh")
expMeanl     <- mxAlgebra( expression= Mean + cbind(bl%*%defM,bl%*%defM), name="expMeanl")
expMeanhl    <- mxAlgebra( expression= Mean + cbind(bh%*%defM,bl%*%defM), name="expMeanhl")

# Thresholds 
thinGh     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0.8, labels=labThh, name="thinGh" )
thinGl     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0.8, labels=labThl, name="thinGl" )
thinGhl     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0.8, labels=labThhl, name="thinGhl" )
 

# unmoderated coefficientienten
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 compn
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" )

# variances + constraints
covPIh      <- mxAlgebra( expression=AIh+CIh+EIh, name="VIh" )
covPIl      <- mxAlgebra( expression=AIl+CIl+EIl, name="VIl" )
var1h     <- mxConstraint(expression = ((ah+ 0%*%aIh) %*% t(ah+ 0%*%aIh))+
                                       ((ch+ 0%*%cIh) %*% t(ch+ 0%*%cIh))+
                                       ((eh+ 0%*%eIh) %*% t(eh+ 0%*%eIh)) == 1, name = "Var1h")
var1l     <- mxConstraint(expression = ((al+ 0%*%aIl) %*% t(al+ 0%*%aIl))+
                                       ((cl+ 0%*%cIl) %*% t(cl+ 0%*%cIl))+
                                       ((el+ 0%*%eIl) %*% t(el+ 0%*%eIl)) == 1, name = "Var1l")

# covariances mzMZ
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="thinGh" )
expMZl     <- mxExpectationNormal( covariance="expCovMZl", means="expMeanl", dimnames=selVars, thresholds="thinGl" )
expMZhl     <- mxExpectationNormal( covariance="expCovMZhl", means="expMeanhl", dimnames=selVars, thresholds="thinGhl" )

expDZh     <- mxExpectationNormal( covariance="expCovDZh", means="expMeanh", dimnames=selVars, thresholds="thinGh" )
expDZl     <- mxExpectationNormal( covariance="expCovDZl", means="expMeanl", dimnames=selVars, thresholds="thinGl" )
expDZhl     <- mxExpectationNormal( covariance="expCovDZhl", means="expMeanhl", dimnames=selVars, thresholds="thinGhl" )

funML     <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
 
parsh <- list(pathBh,mean,thinGh,pathAh, pathCh, pathEh,pathAIh, pathCIh, pathEIh,
              covAIh,covCIh,covEIh,covPIh,var1h)
parsl <- list(pathBl,mean,thinGl,pathAl, pathCl, pathEl,pathAIl, pathCIl, pathEIl,
              covAIl,covCIl,covEIl,covPIl,var1l)
parshl <- list(thinGhl,mean)

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)