#variables
nv        <- 1            
ntv       <- nv*2           
nth	    <- 3		

selVars   <- c('var1', 'var2')
defVars <- c('age1','age2')
useVars  <- c(selVars,defVars)

nndata$var1 	<-mxFactor(nndata$var1, levels=c(1:4) )
nndata$var2 	<-mxFactor(nndata$var2, levels=c(1:4) )

#==================================================================================================
# Saturated Model
#==================================================================================================

#starting values for thresholes
stThresh <-c(-1,1,1.5,-1,1,1.5)

# Age corrections
# Regression effects
beta1        <- mxMatrix( type="Full", nrow=nth, ncol=1, free=TRUE,
                          values=.1, labels='betaAge', name="beta1" )

# Independent variables
obsAge       <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE,
                          labels=c("data.age1","data.age2" ), name="Age")

# threshs - fixed
mean         <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="mean" )

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

inclusions   <- list(beta1, obsAge, mean, add)

#means Monozygotic female
mzfdata     <- subset(nndata, zyg==1, useVars)
dataMZf     <- mxData( observed=mzfdata, type="raw" )

#Expected Threshold Monozygotic female
threshfMZ   <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=TRUE, values=stThresh, labels=c("thresh1a","thresh1b",
                "thresh1c","thresh1a1","thresh1b1","thresh1c1"), lbound=c(-3,.001,.001), ubound=3, name="intercept_MZf" )
expThMZf  	<- mxAlgebra( Low%*%intercept_MZf + beta1 %x% Age, name="expThreshfMZ")

#Covariance  Monozygotic female
covMZf     	<- mxMatrix( type='Stand', nrow=ntv, ncol=ntv, free=TRUE, values=0.3, labels=c('cor1'), name='expCovMZf' )

#Expectation Objects
objfMZ    	<- mxExpectationNormal( covariance="expCovMZf", means="mean", dimnames=selVars, thresholds="expThreshfMZ" )
modelMZf    <- mxModel( "MZf", inclusions, threshfMZ, expThMZf, mean, covMZf, dataMZf, objfMZ, funML )


#and so on the similar for other zygosity groups......

#Putting it all together
multi   <- mxFitFunctionMultigroup(c("MZf", "MZm", "DZf", "DZm", "DZfm", "DZmf"))
ciVC	  <-mxCI(c('cor1', 'cor2', "cor3","cor4","cor5","cor6"))
baseModel <- mxModel("baseModel" , multi, modelMZf,modelMZm,modelDZf,modelDZm,modelDZfm,modelDZmf, ciVC)

# RUN MODEL
BaseFit     <- mxTryHardOrdinal(baseModel)
sumBase      <- summary(BaseFit )
sumBase