#-----------------------------------------------------------------
# Script: ThresholdLiab.R 
 
# Univariate Twin model for ordinal data
# Data File:		CASTage8.csv (TEDS Study)
# Phenotypes:		Childhood Aspergers Symptom Test scores, 
#				categorized in 3 ordered classes
# Type of data:		Raw ordinal data
# ZYGOS: 			1=MZ 2=DZ
# -----------------------------------------------------------------

#rm(list=ls())
#ls()

# Open libraries
require(OpenMx)
require(psych)

# Reads data from a csv file (CASTage8.dat) 
Castdata <- read.table ('CASTage8.csv', header=T, sep=",", na.strings=".")
names(Castdata)
summary(Castdata)
str(Castdata)
describe(Castdata)

#Make the integer variables ordered factors
Castdata$OFcast1 	<-mxFactor(Castdata$Ccast1, levels=c(0:2) )
Castdata$OFcast2 	<-mxFactor(Castdata$Ccast2, levels=c(0:2) )

vars		<-('cast')
selVars	<-c('OFcast1' , 'OFcast2')
useVars	<-c('OFcast1' , 'OFcast2', 'age1', 'age2')

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

# To get the Cont Tables per zygosity group
table(mzData$OFcast1, mzData$OFcast2)
table(dzData$OFcast1, dzData$OFcast2)

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

# 1) Specify Saturated Model (Polychoric correlations) with Age effects on the thresholds 
# ---------------------------------------------------------------------------------------

# Define definition variables to hold the Covariates
obsAge	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, 
					labels=c("data.age1","data.age2"), name="Age")

# Matrix & Algebra for expected means (SND), Cor, Thresholds, Age effect on Th
Mean		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )

betaA		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=c(-.2), 
					labels=c('BaTH'), name="BageTH" )
add		<-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low") 

Tmz		<-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=c(.8, 1), 
					lbound=c(-3,.001), ubound=3,
					labels=c("Tmz11","imz11","Tmz12","imz12"), name="ThMZ")
ThresMZ	<-mxAlgebra( expression= Low%*%ThMZ + BageTH%x%Age, name="expThresMZ")

CorMZ		<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=.6, 
					lbound=-.99, ubound=.99, name="expCorMZ") 

Tdz		<-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=c(.8, 1), 
					lbound=c(-3,.001), ubound=3,
					labels=c("Tdz11","idz11","Tdz12","idz12"), name="ThDZ")
ThresDZ	<-mxAlgebra( expression= Low%*%ThDZ + BageTH%x%Age, name="expThresDZ")

CorDZ		<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=.3, 
					lbound=-.99, ubound=.99, 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    	<- mxFIMLObjective( covariance="expCorMZ", means="M", dimnames=selVars, thresholds="expThresMZ" )
objDZ    	<- mxFIMLObjective( covariance="expCorDZ", means="M", dimnames=selVars, thresholds="expThresDZ" )

# Combine Groups
modelMZ	<- mxModel( obsAge, Mean, betaA, add, Tmz, ThresMZ, CorMZ, dataMZ, objMZ, name="MZ" )
modelDZ	<- mxModel( obsAge, Mean, betaA, add, Tdz, ThresDZ, CorDZ, dataDZ, objDZ, name="DZ" )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxAlgebraObjective( "m2LL" )
Conf		<- mxCI (c ('MZ.expCorMZ[2,1]', 'DZ.expCorDZ[2,1]', 'MZ.BageTH[1,1]') )
SatModel 	<- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf )

# -----------------------------------------------------------------------------------------------
# 1) RUN Saturated Model

SatFit	<- mxRun(SatModel, intervals=TRUE)
(SatSumm	<- summary(SatFit))

# Different ways of getting output:
round(SatFit@output$estimate,4)
SatFit$MZ.expCorMZ
SatFit$DZ.expCorDZ

# Generate output using OpenMx functions
(ExpTetraCorMZ	<-mxEval(MZ.expCorMZ, SatFit))
(ExpTetraCorDZ	<-mxEval(DZ.expCorDZ, SatFit))


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


Sub1Model	<- mxModel(SatModel, name="sub1")
Sub1Model	<- omxSetParameters(Sub1Model, labels=c("Tmz11","imz11","Tmz12","imz12"), 
							 newlabels=c("Tmz11","imz11","Tmz11","imz11"), 
							 free=TRUE, values=c(.8, 1), lbound=c(-3), ubound=3)
Sub1Model	<- omxSetParameters(Sub1Model, labels=c("Tdz11","idz11","Tdz12","idz12"), 
							 newlabels=c("Tdz11","idz11","Tdz11","idz11"), 
							 free=TRUE, values=c(.8, 1), lbound=c(-3), ubound=3)
Sub1Fit	<- mxRun(Sub1Model, intervals=TRUE)
(Sub1Summ	<- summary(Sub1Fit))

# Generate Sub1Model Output
round(Sub1Fit@output$estimate,4)
(ExpTetraCorMZ	<-mxEval(MZ.expCorMZ, Sub1Fit))
(ExpTetraCorDZ	<-mxEval(DZ.expCorDZ, Sub1Fit))
mxCompare(SatFit,Sub1Fit)


# -------------------------------------------------------------------------------
# 3) Specify and Run Sub2Model:  equating Thresh across Twins and zygosity groups

Sub2Model	<- mxModel(SatModel, name="sub2")
Sub2Model	<- omxSetParameters(Sub2Model, labels=c("Tmz11","imz11","Tmz12","imz12"), 
							 newlabels=c("Tmz11","imz11","Tmz11","imz11"), 
							 free=TRUE, values=c(.8, 1), lbound=c(-3, .001), ubound=3)
Sub2Model	<- omxSetParameters(Sub2Model, labels=c("Tdz11","idz11","Tdz12","idz12"), 
							 newlabels=c("Tmz11","imz11","Tmz11","imz11"), 
							 free=TRUE, values=c(.8, 1), lbound=c(-3, .001), ubound=3)
Sub2Fit	<- mxRun(Sub2Model, intervals=TRUE)
(Sub2Summ	<- summary(Sub2Fit))

# Generate Sub2Model Output
round(Sub2Fit@output$estimate,4)
(ExpTetraCorMZ	<-mxEval(MZ.expCorMZ, Sub2Fit))
(ExpTetraCorDZ	<-mxEval(DZ.expCorDZ, Sub2Fit))
mxCompare(SatFit,Sub2Fit)


# Print Comparative Fit Statistics for Saturated and all sub Models
# ---------------------------------------------------------------------
SatNested <- list(Sub1Fit, Sub2Fit)
mxCompare(SatFit,SatNested)


#******************************************************************************************
# 4) Specify ACE Model, with ONE overall set of Thresholds with Age effects on thresholds
# -----------------------------------------------------------------------------------------
# Define definition variables to hold the Covariates
obsAge	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.age1","data.age2"), name="Age")

# Matrix & Algebra for expected means (SND), Thresholds and correlation
Mean		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )
betaA		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=-.2, 
				labels=c('BaTH'), name="BageTH" )
add		<-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low") 
TH		<-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=c(.8,1), 
				lbound=c(-3,.001), ubound=3, labels=c("T1","inc1"), name="Th")
Thres		<-mxAlgebra( expression= Low%*%Th + BageTH%x%Age, name="expThres")

# Matrices declared to store a, c, and e Path Coefficients
pathA	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=.8, label="a11", name="a" ) 
pathC	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=.3, label="c11", name="c" )
pathE	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=.6, label="e11", name="e" )
	
# Matrices generated to hold A, C, and E components and total Variance
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" )

# 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)
mUnv	<- 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	<- mxFIMLObjective( covariance="expCovMZ", means="M", dimnames=selVars, thresholds="expThres" )
objDZ	<- mxFIMLObjective( covariance="expCovDZ", means="M", dimnames=selVars, thresholds="expThres" )

# Combine Groups
pars		<- list(obsAge, Mean, betaA, TH, Thres, add, pathA, pathC, pathE, covA, covC, covE, covP, mUnv, varL )
modelMZ	<- mxModel( pars, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ	<- mxModel( pars, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxAlgebraObjective( "m2LL" )
Conf		<- mxCI (c ('MZ.A[1,1]', 'MZ.C[1,1]', 'MZ.E[1,1]') )
AceModel 	<- mxModel( "ACE", modelMZ, modelDZ, minus2ll, obj, Conf)

# ------------------------------------------------------------------------------
# 4) RUN AceModel

AceFit   <- mxRun(AceModel, intervals=TRUE)
(AceSumm  <- summary(AceFit))
round(AceFit@output$estimate,4)

(h2	<-mxEval(MZ.A, AceFit))
(c2	<-mxEval(MZ.C, AceFit))
(e2	<-mxEval(MZ.E, AceFit))

mxCompare (SatFit,AceFit)

AEmod <- omxSetParameters(model=AceFit,labels="c11",free=F,values=0,name="AE")
AEfit <- mxRun(AEmod,intervals=T)
(AEsumm <- summary(AEfit))
mxCompare(base=AceFit,comparison=AEfit)
pchisq(AEfit@output$minimum - AceFit@output$minimum,df=1,lower.tail=F)