# Open libraries
require(OpenMx)
require(psych)

data	<- read.csv ('dutchdataNopartner.csv', header=T)

data$ebent1 	<-mxFactor(data$ebent1, levels=c(0:1) )
data$ebent2 	<-mxFactor(data$ebent2, levels=c(0:1) )
data$Sex1 		<-mxFactor(data$Sex1, levels=c(0:1) )
data$Sex2 		<-mxFactor(data$Sex2, levels=c(0:1) )

data$PHU1 	<- data$PHU1 -1
data$PHU2 	<- data$PHU2 -1
data$PHU1 	<-mxFactor(data$PHU1, levels=c(0:2) )
data$PHU2 	<-mxFactor(data$PHU2, levels=c(0:2) )



selVars	<- c('ebent1','ebent2')
useVars	<- c('ebent1','ebent2', 'Sex1', 'Sex2','Age1', 'Age2', 'PHU1', 'PHU2')


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

# To get the Cont Tables per zygosity group
table(dzData$ebent1, dzData$PHU1)
table(dzData$ebent2, dzData$PHU2)
table(dzData$ebent1, dzData$ebent2)

nv 	<- 1		# number of variables per twin
ntv 	<- nv*2	# number of variables per pair
nth	<- 1		# number of max thresholds (I changed this to 1 since I only have 2 categories (0,1))
Mage	<- mean(data$Age1)		# Mean age of the sample, to get mean-adjusted covariate effecst for age

# 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")
obsSex	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, 
                    labels=c("data.Sex1","data.Sex2"), name="Sex")
obsPHU	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, 
                    labels=c("data.PHU1","data.PHU2"), name="PHU")
#obspartner	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, 
#                    labels=c("data.partner1","data.partner2"), name="partner")


# Matrix & Algebra for expected means (SND), Cor, Thresholds, Age effect on Th 

Un		<-mxMatrix( type="Unit", nrow=1, ncol=ntv, name="Un2" )

MeanL		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )

betaA		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=c(.00), 
                   labels=c('BaTH'), name="BageTH" )
betaS		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=c(.00), 
                   labels=c('BsTH'), name="BsexTH" )
betaPHU		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=c(0), 
                     labels=c('BpTH'), name="BphuTH" )
#betapart		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=F, values=c(0), 
#                   labels=c('BparTH'), name="BpartTH" )


Tmz		<-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=c(2, 2), 
                 lbound=-4, ubound=3,
                 labels=c("Tmz11","Tmz12"), name="ThMZ")
ThresMZ	<-mxAlgebra( expression= ThMZ + BageTH%x%(Age-Un2%x%Mage) + BsexTH%x%Sex + BphuTH%x%PHU , 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(2, 2), 
                 lbound=-4, ubound=3,
                 labels=c("Tdz11","Tdz12"), name="ThDZ")
ThresDZ	<-mxAlgebra( expression= ThDZ + BageTH%x%(Age-Un2%x%Mage) + BsexTH%x%Sex + BphuTH%x%PHU , 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" )
fitFunction <- mxFitFunctionML()

# Combine Groups
modelMZ	<- mxModel( obsAge, obsSex, obsPHU, MeanL, betaA, betaS, betaPHU,  Tmz, ThresMZ, CorMZ, dataMZ, objMZ, Un, fitFunction, name="MZ" )
modelDZ	<- mxModel( obsAge, obsSex, obsPHU,  MeanL, betaA, betaS, betaPHU,  Tdz, ThresDZ, CorDZ, dataDZ, objDZ, Un, fitFunction, 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]','MZ.BsexTH[1,1]', 'MZ.BphuTH[1,1]', 'MZ.BpartTH[1,1]')  )
Conf		<- mxCI (c ('MZ.expCorMZ[2,1]', 'DZ.expCorDZ[2,1]', 'MZ.BageTH[1,1]','MZ.BsexTH[1,1]' , 'MZ.BphuTH[1,1]' ))
SatModel 	<- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf)

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

SatFit	<- mxRun(SatModel, intervals=F) # This was set to TRUE originally
(SatSumm	<- summary(SatFit))

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

# Generate output using OpenMx functions
mxEval(MZ.expCorMZ, SatFit)
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","Tmz12"), 
                              newlabels=c("Tmz11","Tmz11"), 
                              free=TRUE, values=c(2, 2), lbound=-4, ubound=3)
Sub1Model	<- omxSetParameters(Sub1Model, labels=c("Tdz11","Tdz12"), 
                              newlabels=c("Tdz11","Tdz11"), 
                              free=TRUE, values=c(2, 2), lbound=-4, ubound=3)
Sub1Fit	<- mxRun(Sub1Model, intervals=F)
(Sub1Summ	<- summary(Sub1Fit))

# Generate Sub1Model Output
round(Sub1Fit@output$estimate,4)
mxEval(MZ.expCorMZ, Sub1Fit)
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","Tmz12"), 
                              newlabels=c("Tmz11","Tmz11"), 
                              free=TRUE, values=c(2, 2), lbound=-4, ubound=3)
Sub2Model	<- omxSetParameters(Sub2Model, labels=c("Tdz11","Tdz12"), 
                              newlabels=c("Tmz11","Tmz11"), 
                              free=TRUE, values=c(2, 2), lbound=-4, ubound=3)
Sub2Fit	<- mxRun(Sub2Model, intervals=T)
(Sub2Summ	<- summary(Sub2Fit))

# Generate Sub2Model Output
round(Sub2Fit@output$estimate,4)
mxEval(MZ.expCorMZ, Sub2Fit)
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 Sex Partner PHU 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")
obsSex	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, 
                    labels=c("data.Sex1","data.Sex2"), name="Sex")
obsPHU	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, 
                    labels=c("data.PHU1","data.PHU2"), name="PHU")
#obspartner	<- mxMatrix( type="Full", nrow=1, ncol=2, free=F, 
#                        labels=c("data.partner1","data.partner2"), name="partner")

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

MeanL		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )

betaA		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=c(.00), 
                   labels=c('BaTH'), name="BageTH" )
betaS		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=c(.00), 
                   labels=c('BsTH'), name="BsexTH" )
betaPHU		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=c(0), 
                     labels=c('BpTH'), name="BphuTH" )
#betapart		<-mxMatrix( type="Full", nrow=nth, ncol=1, free=F, values=c(0), 
#                   labels=c('BparTH'), name="BpartTH" )

TH		<-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=c(2, 2), 
                lbound=-4, ubound=3,
                labels=c("Th11","Th11"), name="Th")
Thres		<-mxAlgebra( expression= Th + BageTH%x%(Age-Un2%x%Mage) + BsexTH%x%Sex + BphuTH%x%PHU , 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" )
fitFunction <- mxFitFunctionML()

# Combine Groups
pars		<- list(MeanL, betaA, betaS, betaPHU, TH, Thres, pathA, pathC, pathE, covA, covC, covE, covP, mUnv, varL, Un )
modelMZ	<- mxModel( obsAge, obsSex, obsPHU, pars, covMZ, dataMZ, objMZ, fitFunction, name="MZ" )
modelDZ	<- mxModel( obsAge, obsSex, obsPHU, pars, covDZ, dataDZ, objDZ, fitFunction, 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=T)
(AceSumm  <- summary(AceFit))
round(AceFit@output$estimate,4)

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

mxCompare (SatFit,AceFit)