#-------------------------------------------------------------------------------------
# Script: BivThreshLiab.R 
 
# Bivariate Twin model for ordinal data
# Data File:	ADHDiqage.csv
# Phenotypes
# Variable 1: 	IQ categorised in 2 classes: high and low cognitive functioning
# Variable 2: 	ADHD categorised in 3 classes: 0 (unafected), 
#			1 (above TH for a clinical diagnosis), 2 (sub-clinical group) 
# Type of data:	Raw ordinal data
# ZYG: 		1=MZ 2=DZ
# -------------------------------------------------------------------------------------
# Note: in reality you would rather analyse the continuous distributions
 
#rm(list=ls())
#ls()

require(OpenMx)
require(psych)

Bivdata <- read.table ('ADHDiqage.csv', header=T, sep=",", na.strings="NA")

names(Bivdata)
summary(Bivdata)
str(Bivdata)
describe(Bivdata)


#Check the Distribution
#Histograms
par(mfrow=c(1,2))
hist(Bivdata$iq1); hist(Bivdata$iq2)

# categorizing the variables

BivdataOrd		<- data.frame(Bivdata) 
BivdataOrd$Oiq1	<- cut(Bivdata$iq1, breaks=c(-Inf, 115, Inf), right=F, ordered_result=T, labels=c(0,1)) 
BivdataOrd$Oiq2	<- cut(Bivdata$iq2, breaks=c(-Inf, 115, Inf), right=F, ordered_result=T, labels=c(0,1)) 
BivdataOrd$Oadhd1	<- cut(Bivdata$adhd1, breaks=c(-Inf, 4, 6, Inf), right=F, ordered_result=T, labels=c(0,1,2)) 
BivdataOrd$Oadhd2	<- cut(Bivdata$adhd2, breaks=c(-Inf, 4, 6, Inf), right=F, ordered_result=T, labels=c(0,1,2)) 

describe(BivdataOrd)

#Make the integer variables ordered factors (Note:because of the 'cuting' above, this is not necessary here)
BivdataOrd$OFiq1	 	<-mxFactor(BivdataOrd$Oiq1, levels=c(0:1) )
BivdataOrd$OFiq2	 	<-mxFactor(BivdataOrd$Oiq2, levels=c(0:1) )
BivdataOrd$OFadhd1 	<-mxFactor(BivdataOrd$Oadhd1, levels=c(0:2) )
BivdataOrd$OFadhd2 	<-mxFactor(BivdataOrd$Oadhd2, levels=c(0:2) )

describe(BivdataOrd)

vars		<-c('iq','adhd')
selVars	<-c('OFiq1', 'OFadhd1', 'OFiq2', 'OFadhd2' )
useVars	<-c('OFiq1', 'OFadhd1', 'OFiq2', 'OFadhd2', 'age1', 'age2')

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

# To get the CTs for IQ
table(mzData$OFiq1, mzData$OFiq2)
table(dzData$OFiq1, dzData$OFiq2)
# To get the CTs for ADHD
table(mzData$OFadhd1, mzData$OFadhd2)
table(dzData$OFadhd1, dzData$OFadhd2)

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

# 1) Fits a constrained Polychoric correlation model 
# TH same across twins but different across zyg groups
# Age effect is different acros variables, but same across thresholds within variables (if c>2)
# There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric
# ------------------------------------------------------------------------------------------------------------------------------


# CREATE LABELS & START VALUES as objects(to ease specification)
LabThMZ	<-c('Tmz_11','imz_11','Tmz_21','imz_21')			# THs for var 1 and 2 for a twin individual (mz)
LabCorMZ	<-c('r21','rMZ1','MZxtxt','MZxtxt','rMZ2','r21')
LabThDZ	<-c('Tdz_11','idz_11','Tdz_21','idz_21')			# THs for var 1 and 2 for a twin individual (mz)
LabCorDZ	<-c('r21','rDZ1','DZxtxt','DZxtxt','rDZ2','r21')
LabCov	<-c('BageThiq', 'BageThiq', 'BageThadhd', 'BageThadhd') 
ThPat		<-c(T,F,T,T)
StCorMZ	<-c(-.3, .7, -.3,-.3,   .7, -.3)
StCorDZ	<-c(-.3, .4, -.15,-.15, .3, -.3)
StTH		<-c(1.5, 0, -.8, 1.4 )

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

# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )
betaA		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCov, name="BageTH" )

Tmz		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3, labels=LabThMZ, name="ThMZ")
inc		<-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low") 
ThresMZ	<-mxAlgebra( expression= cbind(Low%*%ThMZ + BageTH%x%Age1, Low%*%ThMZ + BageTH%x%Age2), name="expThresMZ")
CorMZ		<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZ, lbound=-.99, ubound=.99, name="expCorMZ") 

Tdz		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3, labels=LabThDZ, name="ThDZ")
ThresDZ	<-mxAlgebra( expression= cbind(Low%*%ThDZ + BageTH%x%Age1, Low%*%ThDZ + BageTH%x%Age2), name="expThresDZ")
CorDZ		<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZ, 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( obsAge1, obsAge2, Mean, betaA, Tmz, inc, ThresMZ, CorMZ, dataMZ, objMZ, name="MZ" )
modelDZ	<- mxModel( obsAge1, obsAge2, Mean, betaA, Tdz, inc, 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]', 'MZ.BageTH[1,1]', 'MZ.BageTH[1,2]') )
SatModel 	<- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf )

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

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

# Generate SatModel Output

round(SatFit@output$estimate,4)

SatFit$MZ.expCorMZ
SatFit$DZ.expCorDZ

(ExpTetraCorMZ	<-mxEval(MZ.expCorMZ, SatFit))
(ExpTetraCorDZ	<-mxEval(DZ.expCorDZ, SatFit))

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

Sub1Model	<- mxModel(SatModel, name="sub1")
Sub1Model	<- omxSetParameters(Sub1Model, labels=LabThDZ, newlabels=LabThMZ, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3 )
Sub1Fit	<- mxRun(Sub1Model, intervals=F)
(Sub1Summ	<- summary(Sub1Fit))
round(Sub1Fit@output$estimate,4)

# Generate Sub1Model Output
(ExpTetraCorMZ	<-mxEval(MZ.expCorMZ, Sub1Fit))
(ExpTetraCorDZ	<-mxEval(DZ.expCorDZ, Sub1Fit))
mxCompare(SatFit,Sub1Fit)


# CREATE LABELS for Cholesky decomposition with a function
laLower	<- function(la,nv) { paste(la,rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") }
LabTh		<-c('T_11','i_11','T_12','i_12')

# 3) Specify ACE Model, with ONE overall set of Thresholds with Age effects on thresholds
# ------------------------------------------------------------------------------------------

# CREATE LABELS for Cholesky decomposition with a function
laLower	<- function(la,nv) { paste(la,rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") }
LabTh		<-c('T_11','i_11','T_12','i_12')

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

# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )
betaA		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCov, name="BageTH" )
Tr		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, labels=LabTh, lbound=c(-3,.001), ubound=3, name="Th")
inc		<-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low") 
Thres		<-mxAlgebra( expression= cbind(Low%*%Th + BageTH%x%Age1, Low%*%Th + BageTH%x%Age2), name="expThres")

# Matrices declared to store a, c, and e Path Coefficients
pathA	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("a",nv), name="a" ) 
pathC	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.3,-.1,.3), label=laLower("c",nv), name="c" )
pathE	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("e",nv), 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(obsAge1, obsAge2, Mean, betaA, Tr, Thres, pathA, pathC, pathE, covA, covC, covE, covP, mUnv, varL, inc )
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=F)
(AceSumm  <- summary(AceFit))
round(AceFit@output$estimate,4)

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

AceFit$MZ.V

mxCompare(SatFit,AceFit)