# ------------------------------------------------------------------------------
# Program: twinAceBin.R  
#  Author: Hermine Maes
#    Date: 02 20 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
nv        <- 1       # number of variables
ntv       <- nv*2    # number of total variables
selVars   <- c('bmiB1','bmiB2')

# Create binary data and make it an mxFactor
twinData$bmiB1 <- ifelse(twinData$bmi1>22,1.00,0.00)
twinData$bmiB2 <- ifelse(twinData$bmi2>22,1.00,0.00)
twinData$bmiB1 <- mxFactor(twinData$bmiB1, levels=c(0:1) )
twinData$bmiB2 <- mxFactor(twinData$bmiB2, levels=c(0:1) )

# 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$bmiB1,mzDataBin$bmiB2)
table(dzDataBin$bmiB1,dzDataBin$bmiB2)

# ------------------------------------------------------------------------------
# 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=1.5, label="a11", name="a" ) 
pathC    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.5, label="c11", name="c" )
pathE    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="e11", 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="Zero", nrow=1, ncol=ntv, name="Mean" )
threshold 	<- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.8, label="thre", name="threshold" )
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" )

# Create Algebra for expected Threshold Matrices
expThresh   <- mxAlgebra( expression= threshold + cbind(b%*%Age,b%*%Age), 
			name="expThresh" )

# 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="Mean", dimnames=selVars,
				thresholds="expThresh" )
expDZ     	<- 	mxExpectationNormal( covariance="expCovDZ", means="Mean", dimnames=selVars,
				thresholds="expThresh" )


# Combine Groups
pars     <- list( pathA, pathC, pathE, covA, covC, covE, covP, estVars )
bits 		<-c(mean, threshold, defAge, pathB, expThresh)
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','MZ.expThresh'))

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

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

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