# -----------------------------------------------------------------------
# Program: UnivariateTwinAnalysis_MatrixRawCon.R  
#  Author: Hermine Maes
#    Date: 08 01 2009 
#
# Univariate Twin Analysis model to estimate causes of variation (ADE) for continuous data
# Matrix style model input - Raw data input
#
# Revision History
#   Hermine Maes -- 02 01 2010 updated & reformatted
#   modified for heterogeneity session by dd
#   male data
# -----------------------------------------------------------------------

#make sure GenEpiHelperFunctions in in the proper working directory
setwd("C:/Documents and Settings/Danielle Dick/My Documents/R/OpenMxClass/heterogeneity")


require(OpenMx)
source("GenEpiHelperFunctions.R")

# Prepare Data
# -----------------------------------------------------------------------
data(twinData)
summary(twinData)
Vars <- 'bmi'
nv <- 1
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") 
ntv <- nv*2
mzData <- subset(twinData, zyg==2, selVars)
dzData <- subset(twinData, zyg==4, selVars)

# Print Descriptive Statistics
# ----------------------------------------------------------
mean(twinData[twinData$zyg==2,'bmi1'], na.rm=TRUE)  #MZm
mean(twinData[twinData$zyg==4,'bmi1'], na.rm=TRUE)  #DZm

cor(twinData[twinData$zyg==2,'bmi1'], twinData[twinData$zyg==2,'bmi2'], use='pairwise.complete.obs')  #MZm
cor(twinData[twinData$zyg==4,'bmi1'], twinData[twinData$zyg==4,'bmi2'], use='pairwise.complete.obs')  #DZm

# Fit ADE Model with RawData and Matrices Input
# -----------------------------------------------------------------------
univADEModel <- mxModel("univADE",
    mxModel("ADE",
    # Matrices a, d, and e to store a, d, and e path coefficients
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ), #X
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="d11", name="d" ), #Y
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" ), #Z
    # Matrices A, D, and E compute variance components
        mxAlgebra( expression=a %*% t(a), name="A" ),
        mxAlgebra( expression=d %*% t(d), name="D" ),
        mxAlgebra( expression=e %*% t(e), name="E" ),
    # Algebra to compute total variances and standard deviations (diagonal only)
        mxAlgebra( expression=A+D+E, name="V" ),
        mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
        mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"),
   ## Note that the rest of the mxModel statements do not change for bivariate/multivariate case
    # Matrix & Algebra for expected means vector
        mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 20, label="mean", name="Mean" ),
        mxAlgebra( expression= cbind(Mean,Mean), name="expMean"),
    # Algebra for expected variance/covariance matrix in MZ
        mxAlgebra( expression= rbind  ( cbind(A+D+E , A+D),
                                        cbind(A+D   , A+D+E)), name="expCovMZ" ),
    # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
        mxAlgebra( expression= rbind  ( cbind(A+D+E     , 0.5%x%A+0.25%x%D),
                                        cbind(0.5%x%A+0.25%x%D , A+D+E)),  name="expCovDZ" ) 
    ),

    mxModel("MZ",
        mxData( observed=mzData, type="raw" ),
        mxFIMLObjective( covariance="ADE.expCovMZ", means="ADE.expMean", dimnames=selVars )
    ),
    mxModel("DZ", 
        mxData( observed=dzData, type="raw" ),
        mxFIMLObjective( covariance="ADE.expCovDZ", means="ADE.expMean", dimnames=selVars ) 
    ),
    mxAlgebra( expression=MZ.objective + DZ.objective, name="-2sumll" ),
    mxAlgebraObjective("-2sumll")
)

univADEFit <- mxRun(univADEModel)
univADESumm <- summary(univADEFit)
univADESumm

# Generate ADE Output
# -----------------------------------------------------------------------
parameterSpecifications(univADEFit)
expectedMeansCovariances(univADEFit)
tableFitStatistics(univADEFit)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesADE <- mxEval(cbind(ADE.a,ADE.d,ADE.e), univADEFit)
    rownames(pathEstimatesADE) <- 'pathEstimates'
    colnames(pathEstimatesADE) <- c('a','d','e')
pathEstimatesADE

varComponentsADE <- mxEval(cbind(ADE.A/ADE.V,ADE.D/ADE.V,ADE.E/ADE.V), univADEFit)
    rownames(varComponentsADE) <- 'varComponents'
    colnames(varComponentsADE) <- c('%a^2','%d^2','%e^2')
varComponentsADE