#---------------------------------------------------------------------------
# Script: ThreshLiab.R 
 
# Univariate Twin model for ordinal data with two thresholds (3 categories)
# Data File:		CASTage8.dat (TEDS Study)
# Phenotypes:		Childhood Autism Spectrum Test scores, categorized
# Type of data:		Raw ordinal data
# ZYG: 			1=MZ 2=DZ
# ------------------------------------------------------------------------------------------------------------------------

rm(list=ls())

require(OpenMx)
require(psych)

mxOption(NULL, "Default optimizer", "NPSOL")

Castdata <- read.table ('CASTage8.csv', header=T, sep=",", na.strings=".")
summary(Castdata)
str(Castdata)
describe(Castdata)


#Make the integer variables ordered factors
Castdata$Ocast1 	<-mxFactor(Castdata$Ocast1, levels=c(0:2) )
Castdata$Ocast2 	<-mxFactor(Castdata$Ocast2, levels=c(0:2) )
describe(Castdata)

nv 	<- 1		# number of variables per twin
ntv 	<- nv*2	# number of variables per pair
nth	<- 2		# number of max thresholds

selVars	<- c('Ocast1' , 'Ocast2')

# Select Data for Analysis
mzData <- subset(Castdata, zyg==1, selVars)
dzData <- subset(Castdata, zyg==2, selVars)
describe(mzData)
describe(dzData)

# get CT for Ordinal variable
table(mzData$Ocast1, mzData$Ocast2)
table(dzData$Ocast1, dzData$Ocast2)


# 1) Specify Saturated Model (Tetrachoric correlations) 
# ------------------------------------------------------------------------------------------------------------------------------

# Matrix & Algebra for expected means (SND), Thresholds and correlation
MeanL	<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )

Inc	<-mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="L" )
Tmz	<-mxMatrix(type="Full", nrow=nth, ncol=ntv, free=TRUE, values=c(.9, 1, .9, 1), lbound=c(-3, .001, -3, .001 ), ubound=(3),  
	labels=c("Tmz11","imz11", "Tmz12","imz12"), name="ThMZ" )
ThreMZ	<-mxAlgebra( expression= L %*% ThMZ, name="expThreMZ" )

Tdz	<-mxMatrix(type="Full", nrow=nth, ncol=ntv, free=TRUE, values=c(.9, 1, .9, 1), lbound=c(-3, .001, -3, .001 ), ubound=(3),  
	labels=c("Tdz11","idz11", "Tdz12","idz12"), name="ThDZ" )
ThreDZ	<-mxAlgebra( expression= L %*% ThDZ, name="expThreDZ" )

CorMZ	<-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=.8, lbound=-.999, ubound=.999, name="expCorMZ") 
CorDZ	<-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=.4, lbound=-.999, ubound=.999, name="expCorDZ") 

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


# Objective objects for Multiple Groups
objMZ    <- mxExpectationNormal( covariance="expCorMZ", means="M", dimnames=selVars, thresholds="expThreMZ" )
objDZ    <- mxExpectationNormal( covariance="expCorDZ", means="M", dimnames=selVars, thresholds="expThreDZ" )

fitFunction <- mxFitFunctionML()

# Combine Groups
modelMZ	<- mxModel( MeanL, Inc, Tmz, ThreMZ, CorMZ, dataMZ, objMZ, fitFunction, name="MZ" )
modelDZ	<- mxModel( MeanL, Inc, Tdz, ThreDZ, CorDZ, dataDZ, objDZ, fitFunction, name="DZ" )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxFitFunctionAlgebra( "m2LL" )
Conf		<- mxCI (c ('MZ.expCorMZ[2,1]', 'DZ.expCorDZ[2,1]') )
SatModel	<- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf )

# -------------------------------------------------------------------------------------------------------------------------------
# 1) RUN Saturated Model (Tetrachoric correlations)

#SatFit		<- mxRun(SatModel, intervals=T)
# Note: mxTryHardOrdinal replaces 'mxRun'
# wrapper functions to the function mxTryHard, makes multiple attempts to fit an MxModel 
# until the optimizer yields an acceptable solution or the maximum number of attempts is reached.

SatFit	<- mxTryHardOrdinal( SatModel, greenOK = TRUE, checkHess = FALSE, finetuneGradient=FALSE, exhaustive=TRUE, OKstatuscodes=c(0,1,5,6), wtgcsv=c("prev"),intervals=T )

(SatSum	<- summary(SatFit))

# Generate SatModel Output
round(SatFit@output$estimate,4)

SatFit$MZ$expCorMZ
SatFit$DZ$expCorDZ

mxEval(MZ.expThreMZ, SatFit)
mxEval(MZ.expCorMZ, SatFit)
mxEval(DZ.expThreDZ, SatFit)
mxEval(DZ.expCorDZ, SatFit)


# ---------------------------------------------------------------------------------------------------
# 2) Specify and Run Sub1Model:  equating Thresholds across Twin 1 and Twin 2 within MZ group
# Note: omxSetParameters: function to modify the attributes of parameters in a model
# without having to re-specify the model

Sub1Model	<- mxModel(SatModel, name="sub1")
Sub1Model	<- omxSetParameters(Sub1Model, labels=c("Tmz11","imz11", "Tmz12","imz12"), 
		newlabels=c("Tmz11","imz11", "Tmz11","imz11"), free=T, values=c(.9, .8, .9, .8), 
		lbound=c(-3, .001, -3, .001), ubound=(3))
Sub1Fit	<- mxTryHardOrdinal( Sub1Model, greenOK = TRUE, checkHess = FALSE, finetuneGradient=FALSE, exhaustive=TRUE, OKstatuscodes=c(0,1,5,6), wtgcsv=c("prev"), intervals=T )
(Sub1Sum	<- summary(Sub1Fit))

# Generate Sub1Model Output
round(Sub1Fit@output$estimate,4)
Sub1Fit$MZ$expCorMZ
Sub1Fit$DZ$expCorDZ

mxEval(MZ.expThreMZ, Sub1Fit)
mxEval(MZ.expCorMZ, Sub1Fit)
mxEval(DZ.expThreDZ, Sub1Fit)
mxEval(DZ.expCorDZ, Sub1Fit)

# ---------------------------------------------------------------------------------------------------
# 3) Specify and Run Sub2Model:  equating Thresholds across Twin 1 and Twin 2 within DZ group

Sub2Model	<- mxModel(SatModel, name="sub2")
Sub2Model	<- omxSetParameters(Sub2Model, labels=c("Tdz11","idz11","Tdz12","idz12"), 
		newlabels=c("Tdz11","idz11","Tdz11","idz11"), free=T, values=c(.9, .8, .9, .8), 
		lbound=c(-3, .001, -3, .001), ubound=(3))
Sub2Fit	<- mxTryHardOrdinal( Sub2Model, greenOK = TRUE, checkHess = FALSE, finetuneGradient=FALSE, exhaustive=TRUE, OKstatuscodes=c(0,1,5,6), wtgcsv=c("prev"), intervals=T )
(Sub2Sum	<- summary(Sub2Fit))

# Generate Sub2Model Output
round(Sub2Fit@output$estimate,4)
Sub2Fit$MZ$expCorMZ
Sub2Fit$DZ$expCorDZ

mxEval(MZ.expThreMZ, Sub2Fit)
mxEval(MZ.expCorMZ, Sub2Fit)
mxEval(DZ.expThreDZ, Sub2Fit)
mxEval(DZ.expCorDZ, Sub2Fit)

# ---------------------------------------------------------------------------------------------------
# 4) Specify and Run Sub3Model:  equating Thresholds across Twin 1 and Twin 2 and zygosity groups

Sub3Model	<- mxModel(Sub1Model, name="sub3")
Sub3Model	<- omxSetParameters(Sub3Model, labels=c("Tdz11","idz11","Tdz12","idz12"), 
		newlabels=c("Tmz11","imz11","Tmz11","imz11"), free=T, values=c(.9, .8, .9, .8), 
		lbound=c(-3, .001, -3, .001), ubound=(3))
Sub3Fit	<- mxTryHardOrdinal( Sub3Model, greenOK = TRUE, checkHess = FALSE, finetuneGradient=FALSE, exhaustive=TRUE, OKstatuscodes=c(0,1,5,6), wtgcsv=c("prev"), intervals=T )
(Sub3Sum	<- summary(Sub3Fit))

# Generate Sub3Model Output
round(Sub3Fit@output$estimate,4)
Sub3Fit$MZ$expCorMZ
Sub3Fit$DZ$expCorDZ

mxEval(MZ.expThreMZ, Sub3Fit)
mxEval(MZ.expCorMZ, Sub3Fit)
mxEval(DZ.expThreDZ, Sub3Fit)
mxEval(DZ.expCorDZ, Sub3Fit)




# *********************************************************************
# Print Comparative Fit Statistics for Saturated and all sub Models
# ---------------------------------------------------------------------

(SatNested.fit <- 	rbind(	
			mxCompare(SatFit, Sub1Fit),
			mxCompare(SatFit, Sub2Fit)[2,],
	  		mxCompare(SatFit, Sub3Fit)[2,]    ) )




# ---------------------------------------------------------------------
# 5) Specify ACE Model, with ONE overall set of Thresholds
# Note: There is a constraint on the total variances to be unity
# ---------------------------------------------------------------------
# Matrices declared to store a, c, and e Path Coefficients
pathA	<- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a11", name="a" ) 
pathC	<- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c11", name="c" )
pathE	<- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e11", name="e" )
	
# Matrices generated to hold A, C, and E computed Variance Components
covA	<- mxAlgebra( expression=a %*% t(a), name="A" )
covC	<- mxAlgebra( expression=c %*% t(c), name="C" ) 
covE	<- mxAlgebra( expression=e %*% t(e), name="E" )
covP	<- mxAlgebra( expression=A+C+E, name="V" )

# Matrix & Algebra for expected means vector and expected thresholds
meanL	<- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )
Inc	<- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="L" )
T	<-mxMatrix(type="Full", nrow=nth, ncol=ntv, free=TRUE, values=c(.9, 1, .9, 1), 
	lbound=c(-3, .001, -3, .001 ), ubound=(3),  
	labels=c("T11","i11", "T11","i11"), name="Th" )
Thre	<-mxAlgebra( expression= L %*% Th, name="expThre" )


# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covMZ	<- mxAlgebra( expression= rbind( cbind(A+C+E , A+C),
                                          cbind(A+C   , A+C+E)), name="expCovMZ" )
covDZ	<- mxAlgebra( expression= rbind( cbind(A+C+E     , 0.5%x%A+C),
                                          cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )

# Constraint on total variance of Ordinal variables (A+C+E=1)
matUnv	<- mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv" )
varL		<- mxConstraint( expression=diag2vec(V)==Unv, name="VarL" )

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

# Objective objects for Multiple Groups
objMZ	<- mxExpectationNormal( covariance="expCovMZ", means="M", dimnames=selVars, thresholds="expThre" )
objDZ	<- mxExpectationNormal( covariance="expCovDZ", means="M", dimnames=selVars, thresholds="expThre" )

fitFunction <- mxFitFunctionML()

# Combine Groups
pars		<- list( pathA, pathC, pathE, covA, covC, covE, covP, matUnv )
modelMZ	<- mxModel( pars, meanL, Inc, T, Thre, covMZ, dataMZ, objMZ, fitFunction, name="MZ" )
modelDZ	<- mxModel( pars, meanL, Inc, T, Thre, covDZ, dataDZ, objDZ, fitFunction, name="DZ" )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxFitFunctionAlgebra( "m2LL" )
Conf		<- mxCI (c ('ACE.A[1,1]', 'ACE.C[1,1]', 'ACE.E[1,1]') )
AceModel 	<- mxModel( "ACE", pars, varL, modelMZ, modelDZ, minus2ll, obj, Conf)
# ------------------------------------------------------------------------------
# 5) RUN AceModel

AceFit	<- mxTryHardOrdinal( AceModel, greenOK = TRUE, checkHess = FALSE, finetuneGradient=FALSE, exhaustive=TRUE, OKstatuscodes=c(0,1,5,6), wtgcsv=c("prev"), intervals=TRUE )

(AceSum  	<- summary(AceFit))

round(AceFit@output$estimate,4)

mxEval(MZ.A, AceFit)
mxEval(MZ.C, AceFit)
mxEval(MZ.E, AceFit)