# Program: twins.R  
# Author: Hermine Maes
# Date: 03 03 2014 
# AND
# Program: twinMulAceCon.R  
# Author: Hermine Maes
# Date: 03 04 2014 
# ADAPTED BY CHARLOTTE HUPPERTZ, May 2015

### SIMPLEX SCRIPT, for MZM and DZM twins only ###

rm( list=ls(all=TRUE) ) #clear the working space
setwd( "xxx/" ) #set working directory -> change to yours!

#sink( "Charlotte_Simplex_OUTPUT" ) #combined with the commmand and the end of the script, this will create an output file
#First, run the simplex-part without the sink-command - then, run it with the sink command!

# Load libraries
require( OpenMx )
require( psych )
source( "myFunctions.R" ) #used for labeling in Cholesky-script, e.g. labLower("cm",nv) -> you don't need to specify the labels of the matrix yourself

# ------------------------------------------------------------------------------
# 1) PREPARE DATA
# ------------------------------------------------------------------------------
# Load data
bmidata    <- read.table ( "bmidata.dat", header=T, na=-99,dec="." ) #tell openMx which data to use
head( bmidata, n=10 ) #have a look at the data! Anything strange?
dim( bmidata ) 
describe( bmidata )

Vars       <- c( 'bmi16','bmi17','bmi18' ) #we will use three measurement moments
nv         <- length( Vars )
ntv        <- nv*2
selVars    <- paste( Vars,c(rep(1,nv),rep(2,nv)),sep="" )   
selVars

# Select data for analysis
mzmData    <- subset( bmidata, zyg==1, selVars ) #what are we selecting here?
dzmData    <- subset( bmidata, zyg==2, selVars )

dataMZM    <- mxData( observed=mzmData, type="raw" ) #changes the data to openMx-format
dataDZM    <- mxData( observed=dzmData, type="raw" )

# ------------------------------------------------------------------------------
# 2) SIMPLEX MODEL FOR A, C, AND E COMPONENT
# ------------------------------------------------------------------------------

# Expected Mean Matrices
meanMM     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(20,21,21,20,21,21), 
                        label=c( "meanP1M","meanP2M","meanP3M","meanP1M","meanP2M","meanP3M" ), name="expMeanMM" )

# Create matrices to store a, c, & e path coefficients (t=transmission, i=innovation, m=males)
valuesEm   <- c( 0.8,0.5,0.4,0.7,0.2,0.9 )

pathAtm    <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c( FALSE,TRUE,FALSE,FALSE,TRUE,FALSE ), values=c( 0,1.0,0,0,1.0,0 ), 
                        label=c("atm11","atm21","atm31","atm22","atm32","atm33"), name="atm" )
pathAim    <- mxMatrix( type="Diag",  nrow=nv, ncol=nv, free=TRUE, values=c( 2.0,-0.8,0.8 ), label=c("aim11","aim22","aim33"), name="aim" )
pathEm     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=0.8, label=labLower("em",nv), name="em" )

# Matrices A, C, & E to compute variance components
matI       <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I" )

covAm      <- mxAlgebra( expression=solve(I-atm) %&% (aim %*% t(aim)), name="Am" ) #matrix WITHIN twin
covEm      <- mxAlgebra( expression=em %*% t(em), name="Em" )

# Algebra for expected variance/covariance matrices, WITHIN- and CROSS-TWIN
varM       <- mxAlgebra( expression=Am+Em, name="Vm" ) #WITHIN twin

CovMZM     <- mxAlgebra( expression=rbind( cbind(Vm, Am),                             
                                           cbind(Am, Vm)),       name="expCovMZM" )
CovDZM     <- mxAlgebra( expression=rbind( cbind(Vm, 0.5%x%Am), 
                                           cbind(0.5%x%Am, Vm)), name="expCovDZM" )

# Objective objects for multiple groups
objMZM     <- mxExpectationNormal( covariance="expCovMZM", means="expMeanMM", dimnames=selVars )
objDZM     <- mxExpectationNormal( covariance="expCovDZM", means="expMeanMM", dimnames=selVars )
fitFunction <- mxFitFunctionML( vector=F )

# Combine groups
pars       <- list( pathAtm, pathAim, pathEm, covAm, covEm, varM, matI, fitFunction ) #all objects that are not specific to a zygosity

modelMZM   <- mxModel( pars, meanMM, CovMZM, dataMZM, objMZM, name="MZM" ) #combine all objects into one model for each zygosity
modelDZM   <- mxModel( pars, meanMM, CovDZM, dataDZM, objDZM, name="DZM" )

# Objective  
minus2ll   <- mxAlgebra( expression=MZM.objective + DZM.objective, name="m2LL" )
obj        <- mxFitFunctionAlgebra( "m2LL" )

relativeAm <- mxAlgebra( expression=Am/Vm,name="relAm" ) #the relative contribution of A -> heritability
relativeEm <- mxAlgebra( expression=Em/Vm,name="relEm" )

SimplexAModel   <- mxModel( "SimplexA", pars, modelMZM, modelDZM, minus2ll, obj, relativeAm, relativeEm )

# ------------------------------------------------------------------------------
# RUN SIMPLEX AE MODEL
SimplexAFit    <- mxTryHard( SimplexAModel, extraTries=5 )
summary( SimplexAFit )