# ------------------------------------------------------------------------------
# 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("miFunctions2.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
lbVas[lower.tri(lbVas)] <- -10     # lower bounds for below diagonal elements
lbVas[upper.tri(lbVas)] <- NA      # lower bounds for above diagonal elements

# 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
cholMZ    <- mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvMZ, name="cholMZ" )
cholDZ    <- mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvDZ, name="cholDZ" )
covMZ     <- mxAlgebra( expression=cholMZ %*% t(cholMZ), name="covMZ" )
covDZ     <- mxAlgebra( expression=cholDZ %*% t(cholDZ), 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, cholMZ, covMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( "DZ", meanDZ, cholDZ, 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)
fitExpc2(fit)

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