coffOrd          <- data.frame(coff) 
summary(coffOrd)

#Make the integer variables ordered factors 
coffOrd$varB_1 <- mxFactor(coffOrd$varB_1, levels=c(0,1))   
coffOrd$varB_2 <- mxFactor(coffOrd$varB_2, levels=c(0,1))
coffOrd$twzyg_v2 <- mxFactor(coffOrd$twzyg_v2, levels=c(1:5))
summary(coffOrd)

vars    <-c('varA_','varB_')
selVars <-c('varA_1', 'varB_1', 'varA_2', 'varB_2')
useVars <-c('varA_1', 'varB_1', 'varA_2', 'varB_2', 'age_cat_v2_1', 'age_cat_v2_2')

# Select Data for Analysis
mzmData <- subset(coffOrd, twzyg_v2==1, useVars)
dzmData <- subset(coffOrd, twzyg_v2==2, useVars)
mzfData <- subset(coffOrd, twzyg_v2==3, useVars)
dzfData <- subset(coffOrd, twzyg_v2==4, useVars)
dosmfData <- subset(coffOrd, twzyg_v2==5, useVars)

summary(mzmData)                     
summary(dzmData)                    
summary(mzfData)                     
summary(dzfData)                     
summary(dosmfData)   

# To get the Cross Tables for ordinal variable
table(mzmData$varB_1, mzmData$varB_2)
table(dzmData$varB_1, dzmData$varB_2)
table(mzfData$varB_1, mzfData$varB_2)
table(dzfData$varB_1, dzfData$varB_2)
table(dosmfData$varB_1, dosmfData$varB_2)

nv      <- 2        # number of variables per twin
nvo     <- 1        # number of ordinal variables per twin
nvc     <- nv-nvo   # number of continuous variables per twin
poso    <- nvc+1    # position where ordinal variables start
ntv     <- nv*2     # number of variables per pair
nth     <- 1        # number of max thresholds
ninc    <- nth-1    # number of max increments
ncovariates <- 1    # number of covariates

## -----------------------------------------------------------------------------------------
# Part 1:
# Constrained Polyserial correlation model 
## ------------------------------------------------------------------------------------------------------------------------------

# CREATE LABELS in objects(to ease specification)
# I constrain across twin order, except for the DOS group because there are mean gender differences

LabThMZM <-c('Tmzm_11','imzm_11')          # THs for ordinal var for a twin individual (mzm)
LabThDZM <-c('Tdzm_11','idzm_11')          # THs for ordinal var for a twin individual (dzm)
LabThMZF <-c('Tmzf_11','imzf_11')          # THs for ordinal var for a twin individual (mzf)
LabThDZF <-c('Tdzf_11','idzf_11')          # THs for ordinal var for a twin individual (dzf)
LabThDOSMF1 <-c('Tdosmf_11','idosmf_11')          # THs for ordinal var for a twin individual (dosmf)
LabThDOSMF2 <-c('Tdosmf_22','idosmf_22')          # THs for ordinal var for a twin individual (dosmf)

LabCorMZM    <-c('r21','rMZM1','MZMxtxt','MZMxtxt','rMZM2','r21')
LabCorDZM    <-c('r21','rDZM1','DZMxtxt','DZMxtxt','rDZM2','r21')
LabCorMZF    <-c('r21','rMZM1','MZMxtxt','MZMxtxt','rMZM2','r21')
LabCorDZF    <-c('r21','rDZM1','DZMxtxt','DZMxtxt','rDZM2','r21')
LabCorDOSMF    <-c('r21','rDZM1','DZMxtxt','DZMxtxt','rDZM2','r21')

(LabCovMM    <-c( paste("BcvM",1:nvc,sep=""), rep(NA, nvo)))
(LabCovTHM   <-rep( paste("BovM",1:nvo,sep=""), nth)) 

(LabCovMF    <-c( paste("BcvF",1:nvc,sep=""), rep(NA, nvo)))
(LabCovTHF   <-rep( paste("BovF",1:nvo,sep=""), nth)) 

(LabSD  <- c( paste("sd",1:nv,sep="")))
(PatSD  <- c( rep(T,nvc), rep(F, nvo)))
(LabMZM       <- c( paste("mmzm_",1:nv,sep="")))
(PatMZM       <- c( rep(T,nvc), rep(F, nvo)))
(PatBetaMZM   <- c( rep(T,nvc), rep(F, nvo)))

(LabDZM       <- c( paste("mdzm_",1:nv,sep="")))
(PatDZM       <- c( rep(T,nvc), rep(F, nvo)))
(PatBetaDZM   <- c( rep(T,nvc), rep(F, nvo)))

(LabMZF       <- c( paste("mmzf_",1:nv,sep="")))
(PatMZF       <- c( rep(T,nvc), rep(F, nvo)))
(PatBetaMZF   <- c( rep(T,nvc), rep(F, nvo)))

(LabDZF       <- c( paste("mdzf_",1:nv,sep="")))
(PatDZF       <- c( rep(T,nvc), rep(F, nvo)))
(PatBetaDZF   <- c( rep(T,nvc), rep(F, nvo)))

(LabDOSMF1       <- c( paste("mdosmf1_",1:nv,sep="")))
(LabDOSMF2       <- c( paste("mdosmf2_",1:nv,sep="")))
(PatDOSMF       <- c( rep(T,nvc), rep(F, nvo)))
(PatBetaDOSMF   <- c( rep(T,nvc), rep(F, nvo)))

# CREATE START VALUES in objects
#(StM       <-c( colMeans(mzmData[,1:nvc],na.rm=TRUE), rep(0,nvo)))
StM     <-c( 100, 0)
(StSD   <-c( sd(dzmData[,1:nvc],na.rm=TRUE), rep(1,nvo)) )
(StBetaM    <-c( rep(.2,nvc), rep(0,nvo)))

StCorMZM <-c(.3, .7, .3, .3,   .7, .3)
StCorDZM <-c(.3, .4, .15,.15,  .3, .3)
StCorMZF <-c(.3, .7, .3, .3,   .7, .3)
StCorDZF <-c(.3, .4, .15, .15, .3, .3)
StCorDOSMF <-c(.3, .4, .15, .15, .3, .3)
StTH    <-c(.8, 0.9)

# Define definition variables to hold the Covariates
obsAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age_cat_v2_1"), name="age_cat_v2_1")
obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age_cat_v2_2"), name="age_cat_v2_2")

# Matrix & Algebra for expected observed means, Thresholds, effect of Age on Th and correlations
MeanMZM    <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatMZM, values=StM, labels=LabMZM, name="MMZM" )
MeanDZM    <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatDZM, values=StM, labels=LabDZM, name="MDZM" )
MeanMZF    <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatMZF, values=StM, labels=LabMZF, name="MMZF" )
MeanDZF    <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatDZF, values=StM, labels=LabDZF, name="MDZF" )
MeanDOSMF1    <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatDOSMF, values=StM, labels=LabDOSMF1, name="MDOSMF1" )
MeanDOSMF2    <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatDOSMF, values=StM, labels=LabDOSMF2, name="MDOSMF2" )

# MALE
# Effect age on ordinal variable
betaAthM <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTHM, name="BageTHM" )
# Effect age on continuous variable
betaAmM  <-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaMZM, values=StBetaM, labels=LabCovMM, name="BageMM" ) 

#FEMALE
# Effect age on ordinal variable
betaAthF <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTHF, name="BageTHF" )
# Effect age on continuous variable
betaAmF  <-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaMZM, values=StBetaM, labels=LabCovMF, name="BageMF" ) 

Mumzm1 <-mxAlgebra( expression= MMZM + t(BageMM%*%age_cat_v2_1), name="MeanmzmCtw1" )
Mumzm2 <-mxAlgebra( expression= MMZM + t(BageMM%*%age_cat_v2_2), name="MeanmzmCtw2" )
Mumzm12    <-mxAlgebra( expression= cbind(MeanmzmCtw1,MeanmzmCtw2), name="ExpMeanMZM" )

Mudzm1 <-mxAlgebra( expression= MDZM + t(BageMM%*%age_cat_v2_1), name="MeandzmCtw1" )
Mudzm2 <-mxAlgebra( expression= MDZM + t(BageMM%*%age_cat_v2_2), name="MeandzmCtw2" )
Mudzm12    <-mxAlgebra( expression= cbind(MeandzmCtw1,MeandzmCtw2), name="ExpMeanDZM" )

Mumzf1 <-mxAlgebra( expression= MMZF + t(BageMF%*%age_cat_v2_1), name="MeanmzfCtw1" )
Mumzf2 <-mxAlgebra( expression= MMZF + t(BageMF%*%age_cat_v2_2), name="MeanmzfCtw2" )
Mumzf12    <-mxAlgebra( expression= cbind(MeanmzfCtw1,MeanmzfCtw2), name="ExpMeanMZF" )

Mudzf1 <-mxAlgebra( expression= MDZF + t(BageMF%*%age_cat_v2_1), name="MeandzfCtw1" )
Mudzf2 <-mxAlgebra( expression= MDZF + t(BageMF%*%age_cat_v2_2), name="MeandzfCtw2" )
Mudzf12    <-mxAlgebra( expression= cbind(MeandzfCtw1,MeandzfCtw2), name="ExpMeanDZF" )

Mudosmf1 <-mxAlgebra( expression= MDOSMF1 + t(BageMM%*%age_cat_v2_1), name="MeandosmfCtw1" )
Mudosmf2 <-mxAlgebra( expression= MDOSMF2 + t(BageMF%*%age_cat_v2_2), name="MeandosmfCtw2" )
Mudosmf12    <-mxAlgebra( expression= cbind(MeandosmfCtw1,MeandosmfCtw2), name="ExpMeanDOSMF" )

# thresholds

inc       <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")

Tmzm      <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThMZM, name="ThMZM")
ThresMZM  <-mxAlgebra( expression= cbind(Low%*%ThMZM + BageTHM%x%age_cat_v2_1, Low%*%ThMZM + BageTHM%x%age_cat_v2_2), name="ExpThresMZM")

Tdzm      <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThDZM, name="ThDZM")
ThresDZM  <-mxAlgebra( expression= cbind(Low%*%ThDZM + BageTHM%x%age_cat_v2_1, Low%*%ThDZM + BageTHM%x%age_cat_v2_2), name="ExpThresDZM")

Tmzf      <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThMZF, name="ThMZF")
ThresMZF  <-mxAlgebra( expression= cbind(Low%*%ThMZF + BageTHF%x%age_cat_v2_1, Low%*%ThMZF + BageTHF%x%age_cat_v2_2), name="ExpThresMZF")

Tdzf      <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThDZF, name="ThDZF")
ThresDZF  <-mxAlgebra( expression= cbind(Low%*%ThDZF + BageTHF%x%age_cat_v2_1, Low%*%ThDZF + BageTHF%x%age_cat_v2_2), name="ExpThresDZF")

Tdosmf1    <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThDOSMF1, name="ThDOSMF1")
Tdosmf2    <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThDOSMF2, name="ThDOSMF2")
ThresDOSMF <-mxAlgebra( expression= cbind(Low%*%ThDOSMF1 + BageTHM%x%age_cat_v2_1, Low%*%ThDOSMF2 + BageTHF%x%age_cat_v2_2), name="ExpThresDOSMF")

SD  <-mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=c(PatSD,PatSD), values=c(StSD,StSD), labels=c(LabSD,LabSD), name="sd") 
rMZM <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZM, labels=LabCorMZM, lbound=-.99, ubound=.99, name="CorMZM") 
rDZM <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZM, labels=LabCorDZM, lbound=-.99, ubound=.99, name="CorDZM") 
rMZF <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZF, labels=LabCorMZF, lbound=-.99, ubound=.99, name="CorMZF") 
rDZF <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZF, labels=LabCorDZF, lbound=-.99, ubound=.99, name="CorDZF") 
rDOSMF <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDOSMF, labels=LabCorDOSMF, lbound=-.99, ubound=.99, name="CorDOSMF") 

# Matrices for the Covariance model        
MZMCov   <-mxAlgebra( expression=sd %&% CorMZM, name="ExpCovMZM" )
DZMCov   <-mxAlgebra( expression=sd %&% CorDZM, name="ExpCovDZM" )
MZFCov   <-mxAlgebra( expression=sd %&% CorMZF, name="ExpCovMZF" )
DZFCov   <-mxAlgebra( expression=sd %&% CorDZF, name="ExpCovDZF" )
DOSMFCov   <-mxAlgebra( expression=sd %&% CorDOSMF, name="ExpCovDOSMF" )

# Data objects 
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")

# Objective objects for Multiple Groups
objmzm   <-mxFIMLObjective(covariance="ExpCovMZM", means="ExpMeanMZM", dimnames=selVars, thresholds="ExpThresMZM", threshnames=c("varB_1", "varB_2"))
objdzm   <-mxFIMLObjective(covariance="ExpCovDZM", means="ExpMeanDZM", dimnames=selVars, thresholds="ExpThresDZM", threshnames=c("varB_1", "varB_2"))
objmzf   <-mxFIMLObjective(covariance="ExpCovMZF", means="ExpMeanMZF", dimnames=selVars, thresholds="ExpThresMZF", threshnames=c("varB_1", "varB_2"))
objdzf   <-mxFIMLObjective(covariance="ExpCovDZF", means="ExpMeanDZF", dimnames=selVars, thresholds="ExpThresDZF", threshnames=c("varB_1", "varB_2"))
objdosmf   <-mxFIMLObjective(covariance="ExpCovDOSMF", means="ExpMeanDOSMF", dimnames=selVars, thresholds="ExpThresDOSMF", threshnames=c("varB_1", "varB_2"))

# Combine Groups
pars       <-list (obsAge1, obsAge2, betaAthM, betaAmM, betaAthF, betaAmF, inc, SD)
modelMZM <-mxModel(pars, Tmzm, ThresMZM, MeanMZM, Mumzm1, Mumzm2, Mumzm12, rMZM, MZMCov, DataMZM, objmzm, name="MZM" )
modelDZM <-mxModel(pars, Tdzm, ThresDZM, MeanDZM, Mudzm1, Mudzm2, Mudzm12, rDZM, DZMCov, DataDZM, objdzm, name="DZM" )
modelMZF <-mxModel(pars, Tmzf, ThresMZF, MeanMZF, Mumzf1, Mumzf2, Mumzf12, rMZF, MZFCov, DataMZF, objmzf, name="MZF" )
modelDZF <-mxModel(pars, Tdzf, ThresDZF, MeanDZF, Mudzf1, Mudzf2, Mudzf12, rDZF, DZFCov, DataDZF, objdzf, name="DZF" )
modelDOSMF <-mxModel(pars, Tdosmf1, Tdosmf2, ThresDOSMF, MeanDOSMF1, MeanDOSMF2, Mudosmf1, Mudosmf2, Mudosmf12, rDOSMF, DOSMFCov, DataDOSMF, objdosmf, name="DOSMF" )

minus2ll    <-mxAlgebra( expression=MZM.objective + DZM.objective + MZF.objective + DZF.objective + DOSMF.objective, name="m2LL" )
obj         <-mxAlgebraObjective("m2LL" )
Conf1       <-mxCI (c  ('MZM.rMZM[2,1]' )  )                      # cof-smk phenotypic correlations
Conf2       <-mxCI (c  ('MZM.rMZM[3,1]', 'MZM.rMZM[4,2]','MZM.rMZM[3,2]')  )      # MZ twin cor var1, var 2 and xtr-xtwin
Conf3       <-mxCI (c  ('DZM.rDZM[3,1]', 'DZM.rDZM[4,2]','DZM.rDZM[3,2]')  )      # DZ twin cor var1, var 2 and xtr-xtwin
Conf4       <-mxCI (c  ('MZF.rMZF[3,1]', 'MZF.rMZF[4,2]','MZF.rMZF[3,2]')  )      # MZ twin cor var1, var 2 and xtr-xtwin
Conf5       <-mxCI (c  ('DZF.rDZF[3,1]', 'DZF.rDZF[4,2]','DZF.rDZF[3,2]')  )      # DZ twin cor var1, var 2 and xtr-xtwin
Conf6       <-mxCI (c  ('DOSMF.rDOSMF[3,1]', 'DOSMF.rDOSMF[4,2]','DOSMF.rDOSMF[3,2]')  )      # DZ twin cor var1, var 2 and xtr-xtwin
CorModel    <-mxModel( "Cor", modelMZM, modelDZM, modelMZF, modelDZF, modelDOSMF, minus2ll, obj, Conf1, Conf2, Conf3, Conf4, Conf5, Conf6) 

# RUN Correlation MODEL 

CorFit  <-mxRun(CorModel, intervals=F)
CorFit2  <-mxRun(CorFit, intervals=F)
(CorSumm    <-summary(CorFit))

round(CorFit@output$estimate,4)