# Program: UnivSat&ACE_Bin.R  

require(psych)
require(OpenMx)

# Reads data from REC ASCI text file (cast90m.dat) with '.' as missing values, space separated
# Variabels: zyg cast90_tw1 cast90_tw2 
# zyg: 1=mz, 2=dz (males only)

Vars	<- 'Cast'
nv 	<- 1		# number of variables per twin
ntv 	<- nv*2		# number of variables per pair
nthresh1<- 1		# number of max thresholds

allVars<- c('zyg', 'cast90_tw1' , 'cast90_tw2')
Castdata <- read.table ('cast90m.dat', header=F, sep="", na.strings=".",col.names=allVars)
Castdata[,c(2,3)] <- mxFactor(Castdata[,c(2,3)], levels=c(0 : nthresh1) )

# MxFactor enforces ordinal data to be sorted (for OpenMx)
# the ‘levels’ parameter is NOT optional, so that missing levels are not skipped, and leaves the list of levels complete


summary(Castdata)
str(Castdata)
Vars <-'cast90'
selVars <- c('cast90_tw1' , 'cast90_tw2')
mzData <- subset(Castdata, zyg==1, selVars)
dzData <- subset(Castdata, zyg==2, selVars)

# Set Starting Values
thVals    <-1      	# start value for thresholds
corValsM  <-.8    	# start value for MZ correlations
corValsD  <-.4     	# start value for DZ correlations


# Print Descriptive Statistics
# ---------------------------------------------------------------------
summary(mzData)
summary(dzData)
table(mzData$cast90_tw1, mzData$cast90_tw2 )
table(dzData$cast90_tw1, dzData$cast90_tw2)

# ___________________________________
# PREPARE SATURATED MODEL
# ___________________________________

# Matrices for expected Means & Thresholds (on liabilities) in MZ & DZ twin pairs
meanG	<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" )
threMZ	<-mxMatrix(type="Full", nrow=1, ncol=ntv, free=TRUE, values=thVals, labels=c("tMZ1","tMZ2"), name="expThreMZ" )
threDZ	<-mxMatrix(type="Full", nrow=1, ncol=ntv, free=TRUE, values=thVals, labels=c("tDZ1","tDZ2"), name="expThreDZ" )
corMZ	<-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=corValsM, lbound=-.99, ubound=.99, labels=c("rMZ"), name="expCorMZ") 
corDZ	<-mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=T, values=corValsD, lbound=-.99, ubound=.99, labels=c("rDZ"), name="expCorDZ") 

# Data objects for Multiple Groups
dataMZ	<-mxData(mzData, type="raw")
dataDZ	<-mxData(dzData, type="raw")

# Objective objects for Multiple Groups
objMZ	<-mxFIMLObjective( covariance="expCorMZ", means="expMean", dimnames=selVars, thresholds="expThreMZ" )
objDZ	<-mxFIMLObjective( covariance="expCorDZ", means="expMean", dimnames=selVars, thresholds="expThreDZ" )
 
# Combine Groups	
groupMZ	<-mxModel("MZ", corMZ, meanG, threMZ, dataMZ, objMZ )
groupDZ	<-mxModel("DZ", corDZ, meanG, threDZ, dataDZ, objDZ )
minus2ll<-mxAlgebra( MZ.objective + DZ.objective, name="minus2sumloglikelihood" )
obj	<-mxAlgebraObjective("minus2sumloglikelihood") 
ciCor	<-mxCI(c('MZ.expCorMZ[2,1]', 'DZ.expCorDZ[2,1]'))
ciThre	<-mxCI( c('MZ.expThreMZ','DZ.expThreDZ' ))
twinSatModel   <-mxModel( "twinSat", minus2ll, obj, groupMZ, groupDZ, ciCor, ciThre ) 
      
# -----------------------------------------------------------------------
#  RUN SATURATED MODEL (Tetrachoric correlations) 
# -----------------------------------------------------------------------
twinSatFit	<- mxRun(twinSatModel, intervals=F)
twinSatSumm 	<- summary(twinSatFit)
round(twinSatFit@output$estimate,4)
twinSatSumm

# Generate Saturated Model Output
#--------------------------------------------------------------------------
rMZ<- twinSatFit$MZ.expCorMZ@values[2,1]
rDZ<- twinSatFit$DZ.expCorDZ@values[2,1]
tMZ<- twinSatFit$MZ.expThreMZ@values
tDZ<- twinSatFit$DZ.expThreDZ@values

twinSatLLL     <-twinSatFit@output$Minus2LogLikelihood
twinSatAIC     <-twinSatSumm$AIC
twinSatOS      <-twinSatSumm$observedStatistics
twinSatDF      <-twinSatSumm$degreesOfFreedom
twinSatNP      <-length(twinSatSumm$parameters[[1]])

# Use Helper Functions
source("GenEpiHelperFunctions.R")
expectedMeansCovariances(twinSatFit)
tableFitStatistics(twinSatFit)

# -------------------------------------------------------------------------
# RUN SUBMODELS
# -------------------------------------------------------------------------

# SubModel 1: constraining Thresholds across Twin 1 and Twin 2 within zyg group to be equal
# ----------------------------------------------------------------------------------------------------------------------------------------

eqThresholdsTwinModel    <- twinSatFit
eqThresholdsTwinModel    <- omxSetParameters( eqThresholdsTwinModel, label="tMZ1", free=TRUE, values=thVals, newlabels='tMZ' )
eqThresholdsTwinModel    <- omxSetParameters( eqThresholdsTwinModel, label="tMZ2", free=TRUE, values=thVals, newlabels='tMZ' )
eqThresholdsTwinModel    <- omxSetParameters( eqThresholdsTwinModel, label="tDZ1", free=TRUE, values=thVals, newlabels='tDZ' )
eqThresholdsTwinModel    <- omxSetParameters( eqThresholdsTwinModel, label="tDZ2", free=TRUE, values=thVals, newlabels='tDZ' )

eqThresholdsTwinFit     <- mxRun( eqThresholdsTwinModel, intervals=T )
eqThresholdsTwinSumm	<- summary( eqThresholdsTwinFit )
eqThresholdsTwinLLL	<- eqThresholdsTwinFit@output$Minus2LogLikelihood
tableFitStatistics(twinSatFit, eqThresholdsTwinFit)
eqThresholdsTwinSumm

# SubModel 2: constraining Thresholds across Twin 1 and Twin 2 and across zyg group to be equal
# ----------------------------------------------------------------------------------------------------------------------------------------------

eqThresholdsZygModel     <- eqThresholdsTwinModel
eqThresholdsZygModel     <- omxSetParameters( eqThresholdsZygModel, label="tMZ", free=TRUE, values=thVals, newlabels='tZ' )
eqThresholdsZygModel     <- omxSetParameters( eqThresholdsZygModel, label="tDZ", free=TRUE, values=thVals, newlabels='tZ' )

eqThresholdsZygFit      <- mxRun( eqThresholdsZygModel, intervals=F )
eqThresholdsZygSumm	<- summary( eqThresholdsZygFit )
eqThresholdsZygLLL	<- eqThresholdsZygFit@output$Minus2LogLikelihood
tableFitStatistics(eqThresholdsTwinFit, eqThresholdsZygFit)
eqThresholdsZygSumm

# Print Comparative Fit Statistics for Saturated Models
# -----------------------------------------------------------------------------
SatNested <- list(eqThresholdsTwinFit, eqThresholdsZygFit)
tableFitStatistics(twinSatFit, SatNested)


# ___________________________________
# PREPARE GENETIC MODEL
# ___________________________________

# ACE Model with one overall Threshold
# Matrices to store a, c, and e Path Coefficients
pathA    <-mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a11", name="a" ) 
pathC    <-mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c11", name="c" )
pathE    <-mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e11", name="e" )
	
# Algebra to generate Matrices 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" )

# Matrices for expected Means & Thresholds (on liabilities) 
meanG    <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" )
threT    <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=thVals, label="thre", name="expThre" )

# Algebra to compute Total Variance
covP     <-mxAlgebra( expression=A+C+E, name="V" )

# Algebra to generate Matrices to hold Parameter Estimates and Derived Variance Components
rowVars  <-rep('vars',nv)
colVars  <-rep(c('A','C','E','SA','SC','SE'),each=nv)
estVars  <-mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="Vars", dimnames=list(rowVars,colVars))

# 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 variance of the liability of Binary variables (assumed to have a SND) 
matUnv	<-mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv1" )
var1	<-mxConstraint( expression=diag2vec(V)==Unv1, name="Var1" )

# 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="expThre" )
objDZ	<-mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="expThre" )

# Combine Groups
pars	<-list( pathA, pathC, pathE, covA, covC, covE, covP, matUnv, estVars )
modelMZ	<-mxModel( pars, meanG, threT, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ	<-mxModel( pars, meanG, threT, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll<-mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj	<-mxAlgebraObjective( "m2LL" )
ciVC	<-mxCI(c('A', 'C', 'E'))
AceModel <-mxModel( "ACE", pars, var1, modelMZ, modelDZ, minus2ll, obj, ciVC)


# -------------------------------------------------
# RUN GENETIC MODEL
# -------------------------------------------------

# Run ACE Model
AceFit  <- mxRun(AceModel, intervals=T )
AceSumm <- summary(AceFit)
AceSumm
round(AceFit@output$estimate,4)
round(AceFit$Vars@result,4)

# Generate ACE Model Output (using helper function GenEpiHelperFunctions.R)
parameterSpecifications(AceFit)
expectedMeansCovariances(AceFit)
tableFitStatistics(AceFit)

ACEcovMatrices	<- c("A","C","E","V","A/V","C/V","E/V")
ACEcovLabels	<- c("covs_A","covs_C","covs_E","Var","prop_A","prop_C","prop_E")
formatOutputMatrices(AceFit,ACEcovMatrices,ACEcovLabels,Vars,4)

# -----------------------------------------------------------
# RUN SUBMODELS
# -----------------------------------------------------------

# Fit AE model
# -----------------------------------------------------------------
AeModel	<- mxModel( AceFit, name="AE")
AeModel	<- omxSetParameters( AeModel, labels="c11", free=FALSE, values=0 )
AeFit	<- mxRun(AeModel, intervals=T)
AeSumm  <- summary(AeFit)
AeSumm
round(AeFit@output$estimate,4)
round(AeFit$Vars@result,4)


# Run CE model
# ------------------------------------------------------------------
CeModel	<- mxModel( AceFit, name="CE")
CeModel	<- omxSetParameters( CeModel, labels="a11", free=FALSE, values=0 )
CeFit	<- mxRun(CeModel)
round(CeFit@output$estimate,4)
round(CeFit$Vars@result,4)

# ------------------------------------------------------------------------------

# Run E model
# ------------------------------------------------------------------
eModel	<- mxModel( AeFit, name="E")
eModel	<- omxSetParameters( eModel, labels="a11", free=FALSE, values=0 )
eFit	<- mxRun(eModel)
round(eFit@output$estimate,4)
round(eFit$Vars@result,4)

# ------------------------------------------------------------------------------

# Print Comparative Fit Statistics
AceNested	<- list(AeFit, CeFit, eFit)
tableFitStatistics(AceFit,AceNested)

round(rbind(AceFit@output$estimate,AeFit@output$estimate,CeFit@output$estimate,eFit@output$estimate),4)
round(rbind(AceFit$Vars@result,AeFit$Vars@result,CeFit$Vars@result,eFit$Vars@result),4)