# -----------------------------------------------------------------------------------------------------------
# Script: 2ACE_Chol.R  

# Fits a MV Cholesky ACE Decomposition
# Data File:		MVasb.dat (E-Risk Study)
# Phenotypes:		mothers-, teachers-, examiners-, and  children's self-report on antisocial behavior
# Type of data:		Raw continuous data
# ZYG: 			1=MZ 2=DZ
# Assumptions: 		means and variances can be equated across birth order & zygosity groups
# ------------------------------------------------------------------------------------------------------------

require(OpenMx)
require(psych)

# Read in data file 
W5Data <- read.csv("W5OpenMXF.csv") #header=TRUE, na.strings='-99')
W5S <- subset(W5Data, randomTwin == 1) #tells R to only select cases in which the variable = 1
source("GenEpiHelperFunctions.R")

summary(W5S)
describe(W5S)

nv  <- 3    		# number of variables
ntv	<- 2*nv			# number of variables*max family size

Vars	<- c('dfswanhim','dfswanin','LogSCT')
# This yields a list of the 4 variables for twin 1 and twin 2: 
# Mother1   Teacher1   Examiner1   Self1       Mother2   Teacher2   Examiner2   Self2
selVars	<- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")   

mzData <- subset(W5S, zyg==1, selVars)
dzData <- subset(W5S, zyg==2, selVars)


# To create start values 
Stmean <-colMeans(mzData[,1:3],na.rm=TRUE)

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

# To create Labels for Column and Diagonal Matrices
mLabs	<- paste("m",1:nv,sep="")

# -------------------------------------------------------------------------------------------------
# (model 2) Specify Cholesky ACE Decomposition 
# -------------------------------------------------------------------------------------------------

# Matrices declared to store a, c, and e Path Coefficients
pathA  <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,T,T,T,T,T), values=c(.4,.4,.4,.4,.4,.4), labels=c("a11", "a21", "a31", "a22", "a32", "a33"), name="a" )
pathC	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,T,T,T,T,T), values=c(.1,.1,.1,.1,.1,.1), labels=c("c11", "c21", "c31", "c22", "c32", "c33"), name="c" )
pathE	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,T,T,T,T,T), values=c(.2,.2,.2,.2,.2,.2), labels=c("e11", "e21", "e31", "e22", "e32", "e33"), 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" )
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 Correlations
matI	<- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
Rph	<- mxAlgebra( expression=solve(sqrt(I*V)) %&% V , name="Phcor")
Rg	<- mxAlgebra( expression=solve(sqrt(I*A)) %&% A , name="Acor")
Rc	<- mxAlgebra( expression=solve(sqrt(I*C)) %&% C , name="Ccor")
Re	<- mxAlgebra( expression=solve(sqrt(I*E)) %&% E , name="Ecor")

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins

Mean	<- mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=Stmean, labels=c(mLabs,mLabs), name="ExpMean" )

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	<- mxExpectationNormal( covariance="ExpCovMZ", means="ExpMean", dimnames=selVars )
objDZ	<- mxExpectationNormal( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars )

fitFunction <- mxFitFunctionML()

# Combine Groups
pars	<- list( pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, Rph, Rg, Rc, Re) 
modelMZ	<- mxModel( pars, covMZ, Mean, dataMZ, objMZ, fitFunction, name="MZ" )
modelDZ	<- mxModel( pars, covDZ, Mean, dataDZ, objDZ, fitFunction, name="DZ" )
minus2ll<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj	<- mxFitFunctionAlgebra( "m2LL" )
conf1	<- mxCI (c ('MZ.h2[1,1]', 'MZ.h2[2,2]', 'MZ.h2[3,3]', 'MZ.h2[4,4]') )		# h2
conf2	<- mxCI (c ('MZ.c2[1,1]', 'MZ.c2[2,2]', 'MZ.c2[3,3]', 'MZ.c2[4,4]') )		# c2
conf3	<- mxCI (c ('MZ.e2[1,1]', 'MZ.e2[2,2]', 'MZ.e2[3,3]', 'MZ.e2[4,4]') )		# e2
conf4	<- mxCI (c ('MZ.Acor[2,1]','MZ.Acor[3,1]','MZ.Acor[4,1]','MZ.Acor[3,2]','MZ.Acor[4,2]','MZ.Acor[4,3]'))#Rg
conf5	<- mxCI (c ('MZ.Ccor[2,1]','MZ.Ccor[3,1]','MZ.Ccor[4,1]','MZ.Ccor[3,2]','MZ.Ccor[4,2]','MZ.Ccor[4,3]'))#Rc
conf6	<- mxCI (c ('MZ.Ecor[2,1]','MZ.Ecor[3,1]','MZ.Ecor[4,1]','MZ.Ecor[3,2]','MZ.Ecor[4,2]','MZ.Ecor[4,3]'))#Re
CholAceModel<- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj, conf1, conf2, conf3, conf4, conf5, conf6)

#---------------------------------------------------------------------------------------------------------------------
# (model 2) RUN Cholesky Decomposition ACE Model
CholAceFit	<- mxRun(CholAceModel, intervals=F)
(CholAceSum	<- summary(CholAceFit))

#round(CholAceFit@output$estimate,4)

# Print the h2, c2, e2, rG, rC, rE matrices (note: in the summary output you will find them with 95% CI)
mxEval(CholACE.h2, CholAceFit)
mxEval(CholACE.c2, CholAceFit)
mxEval(CholACE.e2, CholAceFit)

mxEval(CholACE.Acor, CholAceFit)
mxEval(CholACE.Ccor, CholAceFit)
mxEval(CholACE.Ecor, CholAceFit)
mxEval(CholACE.Phcor, CholAceFit)