On dealing with non-normality: transform the raw variables or the residuals?

# Load Library
require(OpenMx)
require(psych)

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

# Select Variables for Analysis
Vars <- c('yrsTraining','Age','Sex','MTNcdW95BCnew')
nvInit <- 4
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(Vars,c(rep(1,nvInit),rep(2,nvInit)),sep="")
#=============================
# Select Data for Analysis

mzData <- subset(twinData2015, ZygosityGrp==1, selVars)
dzData <- subset(twinData2015, ZygosityGrp==2, selVars)
colnames(mzData)=colnames(dzData)=c('yrsTraining1','Age1','Sex1','t1','yrsTraining2','Age2','Sex2','t2')

# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")

mzData$t1[is.na(mzData$yrsTraining1)] <- NA
mzData$yrsTraining1[is.na(mzData$yrsTraining1)] <- -999
mzData$t2[is.na(mzData$yrsTraining2)] <- NA
mzData$yrsTraining2[is.na(mzData$yrsTraining2)] <- -999

dzData$t1[is.na(dzData$yrsTraining1)] <- NA
dzData$yrsTraining1[is.na(dzData$yrsTraining1)] <- -999
dzData$t2[is.na(dzData$yrsTraining2)] <- NA
dzData$yrsTraining2[is.na(dzData$yrsTraining2)] <- -999

mzData = mxData(subset(mzData, (!is.na(mzData$t1)|!is.na(mzData$t2))), "raw")
dzData = mxData(subset(dzData, (!is.na(dzData$t1)|!is.na(dzData$t2))), "raw")

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

# ACE Model

pathA <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE,
values=.6, label="a11", name="a" )
pathC <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE,
values=.6, label="c11", name="c" )
pathE <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE,
values=.6, label="e11", name="e" )

covA <- mxAlgebra( a %*% t(a), name="A" )
covC <- mxAlgebra( c %*% t(c), name="C" )
covE <- mxAlgebra( e %*% t(e), name="E" )
covP <- mxAlgebra( A+C+E, name="V" )

# Matrix for expected Mean
intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values= 5, label="mean", name="Mean" )

#***************************************************************

# Matrix for moderating/interacting variable
defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
labels=c("data.Sex1","data.Sex2"), name="Sex")

# Matrices declared to store linear Coefficients for covariate
B_Sex <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values= .01, label="betaSex", name="bSex" )

meanSex <- mxAlgebra( bSex%*%Sex, name="SexR")

#age
defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
labels=c("data.Age1","data.Age2"), name="Age")

# Matrices declared to store linear Coefficients for covariate
B_Age <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values= .01, label="betaAge", name="bAge" )

meanAge <- mxAlgebra( bAge%*%Age, name="AgeR")

#yrsTraining
defYTrain <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
labels=c("data.yrsTraining1","data.yrsTraining2"), name="YTrain")

# Matrices declared to store linear Coefficients for covariate
B_YTrain <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values= .01, label="betaYTrain", name="bYTrain" )

meanYTrain <- mxAlgebra( bYTrain%*%YTrain, name="YTrainR")
#***************************************************************

expMean <- mxAlgebra( Mean + YTrainR + AgeR + SexR , name="expMean")

defs <- list( intercept, defYTrain, B_YTrain, meanYTrain, defSex, B_Sex, meanSex, defAge, B_Age, meanAge)

covMZ <- mxAlgebra(
rbind( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" )

covDZ <- mxAlgebra(
rbind( cbind(A+C+E, 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )

# Objective objects for Multiple Groups

objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=c('t1','t2') )
objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=c('t1','t2') )

funML <- mxFitFunctionML()

rowVars <- rep('vars',nv)
colVars <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estVars <- mxAlgebra(cbind(A,C,E,A/V,C/V,E/V), name="Vars",
dimnames=list(rowVars,colVars))

parameters <- list( pathA, pathC, pathE, covA, covC, covE, covP, estVars)

modelMZ <- mxModel( parameters, defs, expMean, covMZ, mzData, objMZ, funML,
name="MZ" )
modelDZ <- mxModel( parameters, defs, expMean, covDZ, dzData, objDZ, funML,
name="DZ" )

fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction"))

AceModel <- mxModel( "ACE", parameters, modelMZ, modelDZ, fitML, mxCI(c("Vars[1,4:6]")))

# ---------------------------------------------------------------
# RUN MODEL

# Run ACE model
AceFit_cov <- mxRun(AceModel,intervals = T)
AceFitsumm <- summary(AceFit_cov)
AceFitsumm

getMean <- function(i){mxEval(AceFit_cov$MZ$expMean, AceFit_cov$MZ, compute=T, defvar.row=i)}
MZresiduals <- AceFit_cov$MZ$data$observed - t(sapply(1:dim(AceFit_cov$MZ$data$observed)[1], getMean))
getMean <- function(i){mxEval(AceFit_cov$DZ$expMean, AceFit_cov$DZ, compute=T, defvar.row=i)}
DZresiduals <- AceFit_cov$DZ$data$observed - t(sapply(1:dim(AceFit_cov$DZ$data$observed)[1], getMean))