#-------------------------------------------------------------------------------------
# Script: BivOrdCont.R 
 
# Bivariate Twin model for Combined Continuous and Ordered Categorical variables
# Data File:	ADHDiqage.csv
# Phenotypes
# Variable 1: 	IQ continuous variable
# 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: Continuous variable(s) are always selected first, then Ordinal variable(s)
 
#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 ADHD variable

BivdataOrd		<- data.frame(Bivdata) 
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$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('iq1', 'OFadhd1', 'iq2', 'OFadhd2' )
useVars	<-c('iq1', 'OFadhd1', 'iq2', '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 ADHD
table(mzData$OFadhd1, mzData$OFadhd2)
table(dzData$OFadhd1, dzData$OFadhd2)

nv 	<- 2		# number of variables per twin
nvo	<- 1		# number of ordinal variables per twin
nvc	<- nv-nvo	# number of continuous variables per twin
poso	<- nvc+1	# position where ordinal variables start
ntv 	<- nv*2	# number of variables per pair
nth	<- 2		# number of max thresholds
ninc	<- nth-1	# number of max increments

## -----------------------------------------------------------------------------------------
# Part 1:
# Constrained Polyserial correlation model 
# MZ and DZ Thresholds, but constrained across twins
# Age effect is different acros variables, but same across thresholds within ordinal variable(s)if c>2
# There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric
## ------------------------------------------------------------------------------------------------------------------------------

# CREATE LABELS in objects(to ease specification)

LabThMZ	<-c('Tmz_11','imz_11')			# THs for ordinal var for a twin individual (mz)
LabThDZ	<-c('Tdz_11','idz_11')			# THs for ordinal var for a twin individual (dz)

LabCorMZ	<-c('r21','rMZ1','MZxtxt','MZxtxt','rMZ2','r21')
LabCorDZ	<-c('r21','rDZ1','DZxtxt','DZxtxt','rDZ2','r21')

(LabCovM	<-c( paste("Bcv",1:nvc,sep=""), rep(NA, nvo)))
(LabCovTH	<-rep( paste("Bov",1:nvo,sep=""), nth)) 

(LabM		<- c( paste("m",1:nv,sep="")))
(LabSD	<- c( paste("sd",1:nv,sep="")))
(PatSD	<- c( rep(T,nvc), rep(F, nvo)))
(PatM		<- c( rep(T,nvc), rep(F, nvo)))
(PatBetaM	<- c( rep(T,nvc), rep(F, nvo)))

# CREATE START VALUES in objects
#(StM		<-c( colMeans(mzData[,1:nvc],na.rm=TRUE), rep(0,nvo)))
StM		<-c( 100, 0)
(StSD 	<-c( sd(mzData[,1:nvc],na.rm=TRUE), rep(1,nvo)) )
(StBetaM	<-c( rep(.2,nvc), rep(0,nvo)))

StCorMZ	<-c(-.3, .7, -.3,-.3,   .7, -.3)
StCorDZ	<-c(-.3, .4, -.15,-.15, .3, -.3)
StTH		<-c(-.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 observed means, Thresholds, effect of Age on Th and correlations
Mean	<-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatM, values=StM, labels=LabM, name="M" )

betaAm	<-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaM, values=StBetaM, labels=LabCovM, name="BageM" )
betaAth	<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTH, name="BageTH" )

Mu1	<-mxAlgebra( expression= M + t(BageM%*%Age1), name="MeanCtw1" )
Mu2	<-mxAlgebra( expression= M + t(BageM%*%Age2), name="MeanCtw2" )
Mu12	<-mxAlgebra( expression= cbind(MeanCtw1,MeanCtw2), name="ExpMean" )

inc		<-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")
Tmz		<-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThMZ, name="ThMZ")
ThresMZ	<-mxAlgebra( expression= cbind(Low%*%ThMZ + BageTH%x%Age1, Low%*%ThMZ + BageTH%x%Age2), name="ExpThresMZ")

Tdz		<-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThDZ, name="ThDZ")
ThresDZ	<-mxAlgebra( expression= cbind(Low%*%ThDZ + BageTH%x%Age1, Low%*%ThDZ + BageTH%x%Age2), name="ExpThresDZ")

SD	<-mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=c(PatSD,PatSD), values=c(StSD,StSD), labels=c(LabSD,LabSD), name="sd") 
rMZ	<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZ, lbound=-.99, ubound=.99, name="CorMZ") 
rDZ	<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZ, lbound=-.99, ubound=.99, name="CorDZ") 

# Matrices for the Covariance model        
MZCov	<-mxAlgebra( expression=sd %&% CorMZ, name="ExpCovMZ" )
DZCov	<-mxAlgebra( expression=sd %&% CorDZ, name="ExpCovDZ" )

# Data objects 
DataMZ	<-mxData(observed=mzData, type="raw")
DataDZ	<-mxData(observed=dzData, type="raw")

# Objective objects for Multiple Groups
objmz	<-mxFIMLObjective( covariance="ExpCovMZ", means="ExpMean", dimnames=selVars, thresholds="ExpThresMZ", threshnames=c("OFadhd1","OFadhd2"))
objdz	<-mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThresDZ", threshnames=c("OFadhd1","OFadhd2"))

# Combine Groups
pars		<-list (obsAge1,obsAge2,Mean,betaAm,betaAth,Mu1,Mu2,Mu12,inc,SD)
modelMZ	<-mxModel(pars, Tmz, ThresMZ, rMZ, MZCov, DataMZ, objmz, name="MZ" )
modelDZ	<-mxModel(pars, Tdz, ThresDZ, rDZ, DZCov, DataDZ, objdz, name="DZ" )
minus2ll	<-mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<-mxAlgebraObjective("m2LL" )
Conf1		<-mxCI (c  ('MZ.rMZ[2,1]' )  )						# IQ-ADHD phenotypic correlations
Conf2		<-mxCI (c  ('MZ.rMZ[3,1]', 'MZ.rMZ[4,2]','MZ.rMZ[3,2]')  )		# MZ twin cor var1, var 2 and xtr-xtwin
Conf3		<-mxCI (c  ('DZ.rDZ[3,1]', 'DZ.rDZ[4,2]','DZ.rDZ[3,2]')  )		# DZ twin cor var1, var 2 and xtr-xtwin
CorModel	<-mxModel( "Cor", modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2, Conf3) 

# RUN Correlation MODEL 

CorFit	<-mxRun(CorModel, intervals=F)
(CorSumm	<-summary(CorFit))

round(CorFit@output$estimate,4)
CorFit$Cor$MZ.Rph
CorFit$Cor$MZ.Rbmz
CorFit$Cor$DZ.Rbdz


# ********************************************************************************************************

# -------------------------------------------------------------------------------------------------------------------
# (2) Specify Bivariate ACE Model using Cholesky Decomposition, 
# One set of THRESHOLDS (same across twins and zyg group)
# -------------------------------------------------------------------------------------------------------------------

LabTh		<-c('T_11','i_11')			# THs and inc for ordinal var for a twin individual 

# CREATE LABELS FOR Lower Triangular Matrices
(aLabs 	<- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep=""))
(cLabs 	<- paste("c", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep=""))
(eLabs 	<- paste("e", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep=""))

# CREATE START VALUES FOR Lower Triangular Matrices
(SDcon 	<-(sqrt(vech(cov(mzData[,1:nvc],mzData[,1:nvc],use="complete")))))/3
Stpaths	<- (5, -.2, 3)

# 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 observed means, Thresholds, effect of Age on Th and correlations
Mean		<-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatM, values=StM, labels=LabM, name="M" )

betaAm	<-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaM, values=StBetaM, labels=LabCovM, name="BageM" )
betaAth	<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTH, name="BageTH" )

Mu1		<-mxAlgebra( expression= M + t(BageM%*%Age1), name="MeanCtw1" )
Mu2		<-mxAlgebra( expression= M + t(BageM%*%Age2), name="MeanCtw2" )
Mu12		<-mxAlgebra( expression= cbind(MeanCtw1,MeanCtw2), name="ExpMean" )

inc		<-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")
T		<-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabTh, name="Th")
Thres		<-mxAlgebra( expression= cbind(Low%*%Th + BageTH%x%Age1, Low%*%Th + BageTH%x%Age2), name="ExpThres")

# MATRICES FOR THE ACE COV MODEL        
# Matrices to store a, c, and e Path Coefficients
pathA	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=aLabs, name="a" ) 
pathC	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=cLabs, name="c" )
pathE	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=eLabs, 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" )

# Algebra to compute standardized variance components
covP	<- mxAlgebra( expression=A+C+E, name="V" )
StA	<- mxAlgebra( expression=A/V, name="h2")
StC	<- mxAlgebra( expression=C/V, name="c2")
StE	<- mxAlgebra( expression=E/V, name="e2")

# Algebra to compute Phenotypic, A, C & E correlations
matI	<- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
rph	<- mxAlgebra( expression= solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)), name="Rph")
rA	<- mxAlgebra( expression= solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)), name="Ra" )
rC	<- mxAlgebra( expression= solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)), name="Rc" )
rE	<- mxAlgebra( expression= solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)), name="Re" )

# 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" )

# Unity constraint on total variance of Ordinal variable(s)
mUnv	<- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv" )
varL	<- mxConstraint( expression=diag2vec(V[(poso:nv),(poso:nv)])==Unv, name="VarL" )

# Data objects 
DataMZ	<-mxData(observed=mzData, type="raw")
DataDZ	<-mxData(observed=dzData, type="raw")

# Objective objects for Multiple Groups
objmz	<-mxFIMLObjective( covariance="ExpCovMZ", means="ExpMean", dimnames=selVars, thresholds="ExpThres", threshnames=c("OFadhd1","OFadhd2"))
objdz	<-mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThres", threshnames=c("OFadhd1","OFadhd2"))

# Combine Groups
pars1		<-list(obsAge1,obsAge2,Mean,betaAm,betaAth,Mu1,Mu2,Mu12,inc,T,Thres,mUnv,varL)
pars2		<-list(pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, rph, rA, rC, rE )
modelMZ	<- mxModel( pars1, pars2, covMZ, DataMZ, objmz, name="MZ" )
modelDZ	<- mxModel( pars1, pars2, covDZ, DataDZ, objdz, name="DZ" )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxAlgebraObjective( "m2LL" )
Conf1		<- mxCI (c ('h2[1,1]', 'h2[2,2]', 'c2[1,1]', 'c2[2,2]', 'e2[1,1]', 'e2[2,2]') )
Conf2		<- mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Rc[2,1]', 'Re[2,1]') )
AceModel	<- mxModel( "ACE", pars2, modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2)


# ------------------------------------------------------------------------------
# RUN Constrained Bivariate ACE Model

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

AceFit$ACE.h2
AceFit$ACE.c2
AceFit$ACE.e2
AceFit$ACE.Rph
AceFit$ACE.Ra
AceFit$ACE.Rc
AceFit$ACE.Re

# Generate AceModel Output

# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
mxCompare(CorFit,AceFit)