getwd() 
setwd("C:/")

asdtrainingfile <- read.table('asdtrainingfile.dat', header=FALSE, sep="", na.strings="-999.00")

nv 	<- 2		# number of variables per twin
ntv 	<- nv*2	# number of variables per pair
maxth <- 5	

allVars<- c("zygos", "ASD1", "ASD2", "TotalSSP1", "TotalSSP2", "TotalPRCADS1", "TotalPRCADS2", "SSP1", "SSP2", "TotalPRCADSthresh1", "TotalPRCADSthresh2")
asdtrainingfile  <- read.table ('asdtrainingfile.dat', header=F, sep="", na.strings="-999.00",col.names=allVars)
asdtrainingfile [,c(2,3,8,9,10,11)] <- mxFactor(asdtrainingfile [,c(2,3,8,9,10,11)], levels=c(0 : maxth) )

# MxFactor enforces ordinal data to be sorted (for OpenMx)
# the ‘levels’ parameter is NOT optional, so that missing levels are not skipped, and leaves the list of levels complete

summary(asdtrainingfile)
str(asdtrainingfile )

selVars <- c('ASD1' , 'SSP1' , 'ASD2', 'SSP2')
thres	<- c("Th1", "Th2", "Th3", "Th4", "Th5")
mzData <- subset(asdtrainingfile, zygos==1, selVars)
dzData <- subset(asdtrainingfile,  zygos==2, selVars)


# Print Descriptive Statistics
# ---------------------------------------------------------------------
summary(mzData)
summary(dzData)

table(mzData$ASD1, mzData$ASD2)
table(dzData$ASD1, dzData$ASD2)

table(mzData$SSP1, mzData$SSP2)
table(dzData$SSP1, dzData$SSP2) 

# ___________________________________
# PREPARE SATURATED MODEL
# ___________________________________

# Matrices for expected Means & Thresholds (on liabilities) in MZ & DZ twin pairs
meanG	<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" )

threASDmz	<-mxMatrix(type="Full", nrow=5, ncol=1, free=c(F,F,F,F,F), values=c(2.3,0,0,0,0), lbound=c(-4,.001,.001,.001,.001), labels=c('tv1', 'i11', 'i21', 'i31', 'i41'), name="V1thmz" )
threSSPmz	<-mxMatrix(type="Full", nrow=5, ncol=1, free=TRUE, values=c(-2,1,1,1,1), lbound=c(-4,.001,.001,.001,.001), labels=c('tv2', 'i12', 'i22', 'i32', 'i42'), name="V2thmz" )
threASDdz	<-mxMatrix(type="Full", nrow=5, ncol=1, free=c(F,F,F,F,F), values=c(2.3,0,0,0,0), lbound=c(-4,.001,.001,.001,.001), labels=c('tv1', 'i11', 'i21', 'i31', 'i41'), name="V1thdz" )
threSSPdz	<-mxMatrix(type="Full", nrow=5, ncol=1, free=TRUE, values=c(-2,1,1,1,1), lbound=c(-4,.001,.001,.001,.001), labels=c('tv2', 'i12', 'i22', 'i32', 'i42'), name="V2thdz" )

Low1		<-mxMatrix(type="Full", nrow=5, ncol=5, free=FALSE, values=1, name="L1")

threshTmz	<-mxAlgebra(expression = cbind (    (L1%*%V1thmz) , (L1%*%V2thmz),  (L1%*%V1thmz) , (L1%*%V2thmz)    ) , name="expThremz") 
threshTdz	<-mxAlgebra(expression = cbind (    (L1%*%V1thdz) , (L1%*%V2thdz),  (L1%*%V1thdz) , (L1%*%V2thdz)    ) , name="expThredz") 

corPhmz	<-mxMatrix(type="Stand", nrow=nv, ncol=nv, free=T, values=.4, lbound=-.99, ubound=.99, labels="r12", name="expCorPhmz") 
corPhdz	<-mxMatrix(type="Stand", nrow=nv, ncol=nv, free=T, values=.4, lbound=-.99, ubound=.99, labels="r12", name="expCorPhdz") 
corBl2	<-mxMatrix(type="Symm", nrow=nv, ncol=nv, free=c(F,T, T, T),  values= c(.8, .3, .3, .6), lbound=-.99, ubound=.99, name="Block2") 
corBl3	<-mxMatrix(type="Symm", nrow=nv, ncol=nv, free=c(F,T, T, T),  values= c(.4, .3, .3, .6), lbound=-.99, ubound=.99, name="Block3") 

corMZ		<-mxAlgebra(expression = rbind ( cbind ( expCorPhmz, Block2) , cbind ( Block2, expCorPhmz) ) , name="expCorMZ") 
corDZ		<-mxAlgebra(expression = rbind ( cbind ( expCorPhdz, Block3) , cbind ( Block3, expCorPhdz) ) , name="expCorDZ") 

# Data objects for Multiple Groups
dataMZ	<-mxData(mzData, type="raw")
dataDZ	<-mxData(dzData, type="raw")

# Objective objects for Multiple Groups
objMZ	<-mxFIMLObjective( covariance="expCorMZ", means="expMean", dimnames=selVars, thresholds="expThremz" )
objDZ	<-mxFIMLObjective( covariance="expCorDZ", means="expMean", dimnames=selVars, thresholds="expThredz" )
 
# Combine Groups	
groupMZ	<-mxModel("MZ", corMZ, corPhmz, corBl2, meanG, Low1, threASDmz, threSSPmz, threshTmz, dataMZ, objMZ )
groupDZ	<-mxModel("DZ", corDZ, corPhdz, corBl3, meanG, Low1, threASDdz, threSSPdz, threshTdz, dataDZ, objDZ )

minus2ll  <-mxAlgebra( MZ.objective + DZ.objective, name="minus2sumloglikelihood" )
obj	<-mxAlgebraObjective("minus2sumloglikelihood") 
#ciCorPh	<-mxCI (  'MZ.expCorMZ[2,1]'  )
#ciCorMZ	<-mxCI (c  (    'MZ.expCorMZ[4,1]',   'MZ.expCorMZ[4,2]'    )  )
#ciCorDZ	<-mxCI (c  (    'DZ.expCorMZ[4,1]',   'DZ.expCorMZ[4,2]'  )  )

twinSatModel   <-mxModel( "twinSat", minus2ll, obj, groupMZ, groupDZ ) 
      
# -----------------------------------------------------------------------
#  RUN SATURATED MODEL (Tetrachoric correlations) 
# -----------------------------------------------------------------------
twinSatFit	<- mxRun(twinSatModel, intervals=F)
twinSatSumm 	<- summary(twinSatFit)
round(twinSatFit@output$estimate,4)
twinSatSumm


&&&&&&&&&&&&&&&&