# ------------------------------------------------------------------------------
# Program: mulSATc.R  
#  Author: Hermine Maes
#    Date: 02 25 2016 
#
# Twin Multivariate Saturated model to estimate means and (co)variances across multiple groups
# Matrix style model - Raw data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|

# Load Libraries & Options
library(OpenMx)
library(psych)
source("miFunctions.R")

# Create Output 
filename    <- "mulSATc2"
sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE)

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

# Load Data
nl <- read.table ("DHBQ_bs.dat", header=T, na=-999)
describe(nl, skew=F)

# Recode Data for Analysis - Rescale variables to have variances around 1.0
nl$family1  <- nl$gff1/5
nl$family2  <- nl$gff2/5
nl$happy1   <- nl$hap1/4
nl$happy2   <- nl$hap2/4
nl$life1    <- nl$sat1/5
nl$life2    <- nl$sat2/5
nl$anxdep1  <- nl$AD1/4
nl$anxdep2  <- nl$AD2/4
nl$somatic1 <- nl$SOMA1/2
nl$somatic2 <- nl$SOMA2/2
nl$social1  <- nl$SOC1/1.5
nl$social2  <- nl$SOC2/1.5

# Select Variables for Analysis
vars      <- c('family','happy','life','anxdep','somatic','social') 
nv        <- 6       # number of variables
ntv       <- nv*2    # number of total variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Random Subset to reduce time to Fit Examples
testData  <- head(nl,n=500)

# Select Data for Analysis
mzData    <- subset(testData, zyg2==1, selVars)
dzData    <- subset(testData, zyg2==2, selVars)

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

# Set Starting Values 
svMe      <- c(7,5,5,1,1,1)        # start value for means
svVa      <- 1                     # start value for variances
svVas     <- diag(svVa,ntv,ntv)    # assign start values to diagonal of matrix
lbVa      <- .0001                 # start value for lower bounds
lbVas     <- diag(lbVa,ntv,ntv)    # assign lower bounds values to diagonal of matrix

# Create Labels
labMeMZ   <- paste("meanMZ",selVars,sep="_")
labMeDZ   <- paste("meanDZ",selVars,sep="_")
labMeZ    <- paste("meanZ",selVars,sep="_")
labCvMZ   <- labLower("covMZ",ntv)
labCvDZ   <- labLower("covDZ",ntv)
labCvZ    <- labLower("covZ",ntv)
labVaMZ   <- labDiag("covMZ",ntv)
labVaDZ   <- labDiag("covDZ",ntv)
labVaZ    <- labDiag("covZ",ntv)

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

# Saturated Model
# Create Algebra for expected Mean Matrices
meanMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZ, name="meanDZ" )

# Create Algebra for expected Variance/Covariance Matrices
covMZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )

# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="covMZ", means="meanMZ", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="covDZ", means="meanDZ", dimnames=selVars )
funML     <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
modelMZ   <- mxModel( "MZ", meanMZ, covMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( "DZ", meanDZ, covDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )

# Create Confidence Interval Objects
ciCov     <- mxCI( c('MZ.covMZ','DZ.covDZ') )
ciMean    <- mxCI( c('MZ.meanMZ','DZ.meanDZ') )

# Build Saturated Model with Confidence Intervals
model     <- mxModel( "mulSATc", modelMZ, modelDZ, multi, ciCov, ciMean )

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

# Run Saturated Model
fit       <- mxRun( model, intervals=F )
sum       <- summary( fit )
#round(fit$output$estimate,4)
fitGofs(fit)
fitExpc(fit)

# ------------------------------------------------------------------------------
# RUN SUBMODELS

# Constrain expected Means to be equal across twin order
modelEMO  <- mxModel( fit, name="eqMeansTwin" )
for (i in 1:nv) {
modelEMO  <- omxSetParameters( modelEMO, label=c(labMeMZ[nv+i],labMeMZ[i]), free=TRUE, values=svMe[i], newlabels=labMeMZ[i] ) 
modelEMO  <- omxSetParameters( modelEMO, label=c(labMeDZ[nv+i],labMeDZ[i]), free=TRUE, values=svMe[i], newlabels=labMeDZ[i] ) }
fitEMO    <- mxRun( modelEMO, intervals=F )
mxCompare( fit, fitEMO )
fitGofs(fitEMO)
fitExpc(fitEMO)

# Constrain expected Means and Variances to be equal across twin order
modelEMVO <- mxModel( fitEMO, name="eqMVarsTwin" )
for (i in 1:nv) {
modelEMVO <- omxSetParameters( modelEMVO, label=c(labVaMZ[nv+i],labVaMZ[i]), free=TRUE, values=svVa, newlabels=labVaMZ[i] ) 
modelEMVO <- omxSetParameters( modelEMVO, label=c(labVaDZ[nv+i],labVaDZ[i]), free=TRUE, values=svVa, newlabels=labVaDZ[i] ) }
fitEMVO   <- mxRun( modelEMVO, intervals=F )
mxCompare( fit, subs <- list(fitEMO, fitEMVO) )
fitGofs(fitEMVO)
fitExpc(fitEMVO)

# Constrain expected Means and Variances to be equal across twin order and zygosity
modelEMVZ <- mxModel( fitEMVO, name="eqMVarsZyg" )
for (i in 1:nv) {
modelEMVZ <- omxSetParameters( modelEMVZ, label=c(labMeMZ[i],labMeDZ[i]), free=TRUE, values=svMe[i], newlabels=labMeZ[i] ) 
modelEMVZ <- omxSetParameters( modelEMVZ, label=c(labVaMZ[i],labVaDZ[i]), free=TRUE, values=svVa, newlabels=labVaZ[i] ) }
fitEMVZ   <- mxRun( modelEMVZ, intervals=F )
mxCompare( fit, subs <- list(fitEMO, fitEMVO, fitEMVZ) )
fitGofs(fitEMVZ)
fitExpc(fitEMVZ)

# ------------------------------------------------------------------------------                                           
sink()
save.image(paste(filename,".Ri",sep=""))