#---------------------------------------------------------------------------
# Script: UnivM&SDord.R 
 
# Univariate Twin model for ordinal data, Estimate Means & SD of categorical variables
# Only possible when n categories is at least 3 
# 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()

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 Polychoric correlation Model with ONE overall Mean with Age effects 
# Separate SD are estimated for twin1 and 2 in MZ and DZ group 
# Note TH1 fixed to 0, TH2 fixed to 1 in both MZ and DZ group
# However, with c>3, it is possible to estimate the other thresholds
# ------------------------------------------------------------------------------------------------------------------------------

# 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 M
Mean		<-mxMatrix( type="Full", nrow=1, ncol=nv, free=T, values=.8, labels=c("M1"), name="M" )
betaAm	<-mxMatrix( type="Full", nrow=nv, ncol=1, free=T, values=.2, labels=c("BaM"), name="BageM" )

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=ntv, free=F, values=c(0,1), name="ThMZ")
ThresMZ	<-mxAlgebra( expression= Low%*%ThMZ, name="ExpThresMZ")

Tdz		<-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=F, values=c(0,1), name="ThDZ")
ThresDZ	<-mxAlgebra( expression= Low%*%ThDZ, name="ExpThresDZ")

SDmz	<-mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=T, values=1, labels=c("Smz1","Smz2"), name="MZsd") 
rMZ	<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=.7, labels=c("rmz"), lbound=-.99, ubound=.99, name="CorMZ") 
MZCov	<-mxAlgebra( expression=MZsd %&% CorMZ, name="ExpCovMZ" )

SDdz	<-mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=T, values=1, labels=c("Sdz1","Sdz2"), name="DZsd") 
rDZ	<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=.7, labels=c("rdz"), lbound=-.99, ubound=.99, name="CorDZ") 
DZCov	<-mxAlgebra( expression=DZsd %&% CorDZ, name="ExpCovDZ" )

# 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="ExpMean", dimnames=selVars, thresholds="ExpThresMZ" )
objDZ    	<- mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThresDZ" )

# Combine Groups
modelMZ	<- mxModel( obsAge1, obsAge2, Mean, betaAm, Mu1, Mu2, Mu12, inc, Tmz, ThresMZ, SDmz, rMZ, MZCov, dataMZ, objMZ, name="MZ" )
modelDZ	<- mxModel( obsAge1, obsAge2, Mean, betaAm, Mu1, Mu2, Mu12, inc, Tdz, ThresDZ, SDdz, rDZ, DZCov, dataDZ, objDZ, name="DZ" )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxAlgebraObjective( "m2LL" )
Conf		<- mxCI (c ('MZ.CorMZ[2,1]', 'DZ.CorDZ[2,1]', 'MZ.BageM[1,1]') )
SatModel 	<- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf )

# -------------------------------------------------------------------------------------------------------------------------------
# 1) RUN Correlation Model

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

round(SatFit@output$estimate,4)
SatFit$MZ.ThMZ
SatFit$MZ.CorMZ
SatFit$DZ.CorDZ


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


Sub1Model	<- mxModel(SatModel, name="sub1")
Sub1Model	<- omxSetParameters(Sub1Model, labels=c("Smz1","Smz2"), newlabels=c("Smz1","Smz1"), free=T, values=1)
Sub1Model	<- omxSetParameters(Sub1Model, labels=c("Sdz1","Sdz2"), newlabels=c("Sdz1","Sdz1"), free=T, values=1)
Sub1Fit	<- mxRun(Sub1Model, intervals=TRUE)
(Sub1Summ	<- summary(Sub1Fit))
round(Sub1Fit@output$estimate,4)

# Generate Sub1Model Output
mxCompare(SatFit,Sub1Fit)

# ---------------------------------------------------------------------------------------------------
# 3) Specify and Run Sub2Model:  equating SD(Variance) across Twin 1 and Twin 2 and zygosity groups

Sub2Model	<- mxModel(SatModel, name="sub2")
Sub2Model	<- omxSetParameters(Sub2Model, labels=c("Smz1","Smz2"), newlabels=c("Smz1","Smz1"), free=T, values=1)
Sub2Model	<- omxSetParameters(Sub2Model, labels=c("Sdz1","Sdz2"), newlabels=c("Smz1","Smz1"), free=T, values=1)
Sub2Fit	<- mxRun(Sub2Model, intervals=TRUE)
(Sub2Summ	<- summary(Sub2Fit))
round(Sub2Fit@output$estimate,4)

# Generate Sub2Model Output
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 Mean with Age effects 
# ---------------------------------------------------------------------
# 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 M
Mean		<-mxMatrix( type="Full", nrow=1, ncol=nv, free=T, values=.8, labels=c("M1"), name="M" )
betaAm	<-mxMatrix( type="Full", nrow=nv, ncol=1, free=T, values=.2, labels=c("BaM"), name="BageM" )

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=ntv, free=F, values=c(0,1), name="ThMZ")
ThresMZ	<-mxAlgebra( expression= Low%*%ThMZ, name="ExpThresMZ")

Tdz		<-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=F, values=c(0,1), name="ThDZ")
ThresDZ	<-mxAlgebra( expression= Low%*%ThDZ, name="ExpThresDZ")

# Matrices declared to store a, c, and e Path Coefficients
pathA	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.8, label="a11", name="a" ) 
pathC	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.3, label="c11", name="c" )
pathE	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, 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 to compute standardized variance components
StA	<- mxAlgebra( expression=A/V, name="h2")
StC	<- mxAlgebra( expression=C/V, name="c2")
StE	<- mxAlgebra( expression=E/V, name="e2")

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

# 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="ExpMean", dimnames=selVars, thresholds="ExpThresMZ" )
objDZ    	<- mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThresDZ" )

# Combine Groups
pars1		<- list(obsAge1, obsAge2, Mean, betaAm, Mu1, Mu2, Mu12, inc)
pars2		<- list(pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE)
modelMZ	<- mxModel( pars1, pars2, Tmz, ThresMZ, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ	<- mxModel( pars1, pars2, Tdz, ThresDZ, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxAlgebraObjective( "m2LL" )
Conf		<- mxCI (c ('h2[1,1]', 'c2[1,1]', 'e2[1,1]') )
AceModel 	<- mxModel( "ACE", pars2, modelMZ, modelDZ, minus2ll, obj, Conf)

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

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

AceFit$h2
AceFit$c2
AceFit$e2
AceFit$V

mxCompare(SatFit,AceFit)