# ------------------------------------------------------------------------------
# Program: twinAceBin.R  
#  Orignal Author: Hermine Maes
#    Date: 02 20 2012 
#  Edited by: Sarah Medland
#    Date: 7 March 2012 
#
# Univariate Twin Analysis model to estimate causes of variation
# Matrix style model - Raw data - Binary data
# -------|---------|---------|---------|---------|---------|---------|---------|

# Load Library
require(OpenMx)
#require(psych)

# ------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
data(twinData)
#describe(twinData, skew=F)

# Select Variables for Analysis
Vars      <- 'bmiO'
nv        <- 1       # number of variables
ntv       <- nv*2    # number of total variables
selVars   <- c('bmiO1','bmiO2')
nth       <- 2       # number of thresholds

# Create binary data and make it an mxFactor

twinData$bmiO1 <-  cut(twinData$bmi1, breaks=c(18,21,23,30), ordered_result=T, labels=c(0,1,2)) 
twinData$bmiO2 <-  cut(twinData$bmi2, breaks=c(18,21,23,30), ordered_result=T, labels=c(0,1,2)) 

twinData$bmiO1 <- mxFactor(twinData$bmiO1, levels=c(0:2) )
twinData$bmiO2 <- mxFactor(twinData$bmiO2, levels=c(0:2) )

# Select Covariates for Analysis
twinData[,'age'] <- twinData[,'age']/100
twinData  <- twinData[-which(is.na(twinData$age)),]
covVars   <- 'age'

# Select Data for Analysis
mzDataBin    <- subset(twinData, zyg==1, c(selVars, covVars))
dzDataBin    <- subset(twinData, zyg==3, c(selVars, covVars))

# Generate Descriptive Statistics
table(mzDataBin$bmiO1,mzDataBin$bmiO2)
table(dzDataBin$bmiO1,dzDataBin$bmiO2)

# ------------------------------------------------------------------------------
# PREPARE MODEL

# ACE Model
# Matrices declared to store a, c, and e Path Coefficients
pathA    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=2, label="a11", name="a" ) 
pathC    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="c11", name="c" )
pathE    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=2, label="e11",lbound=0,  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" )

# Matrix & Algebra for expected means vector and expected thresholds
mean    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1,label="m", name="Mean" )
defAge    	<- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, 
			labels=c("data.age"), name="Age" )
pathB     	<- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, 
			values=.01, label="b11", name="b" )
correctedMean    <- mxAlgebra( expression=  cbind(Mean+(b%*%Age),Mean+(b%*%Age)), name="correctedMean" )

# Create Algebra for expected Threshold Matrices
Inc      <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" )

#if you have more than 2 thresholds - thresholds 3...n are free
#ie free=c(F,F,T,T)
Thre    <- mxMatrix( type="Full", nrow=nth, ncol=nv, free=c(F,F), 
	values=c(0,1),
	labels=c("t1","inc"), 
	name="Thre" )

threshold    <- mxAlgebra( expression= cbind(Inc %*% Thre,Inc %*% Thre), name="threshold" )


# Algebra to compute total variances and standard deviations (diagonal only)
covP     <- mxAlgebra( expression=A+C+E, name="V" )

# Algebras generated 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 Mean and 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=mzDataBin, type="raw" )
dataDZ   <- mxData( observed=dzDataBin, type="raw" )

# Expectation objects for Multiple Groups
expMZ     	<- 	mxExpectationNormal( covariance="expCovMZ", means="correctedMean", dimnames=selVars,
				thresholds="threshold" )
expDZ     	<- 	mxExpectationNormal( covariance="expCovDZ", means="correctedMean", dimnames=selVars,
				thresholds="threshold" )
pars     	<- list( pathA, pathC, pathE, covA, covC, covE, covP,  estVars )
bits 		<-c(mean, threshold, defAge, pathB, Thre, Inc, correctedMean)
funML     	<- 	mxFitFunctionML()
modelMZ  	<- mxModel( pars, bits, covMZ, dataMZ, funML,expMZ, name="MZ" )
modelDZ  	<- mxModel( pars, bits, covDZ, dataDZ, funML,expDZ, name="DZ" )

# Create Confidence Interval Objects
ciCov     	<- mxCI( c('MZ.a','MZ.c','MZ.e'))

# Combine Groups
multi     	<- mxFitFunctionMultigroup( c("MZ","DZ") )
AceModel 	<- mxModel( "ACE", pars,  modelMZ, modelDZ, ciCov, funML, multi  )

# ------------------------------------------------------------------------------
# RUN MODELS

# Run ACE Model
AceFitFixT   <- mxRun(AceModel, intervals=TRUE)
AceSummFixT  <- summary(AceFitFixT)
AceSummFixT
round(AceFitFixT@output$estimate,4)
round(AceFitFixT$Vars@result,4)