# By Rob Kirkpatrick & Hermine Maes
# Adapted from a script originally written by Timothy C. Bates

require(umx)
umx_set_optimizer("NPSOL")

# =========================================================
# Setup data
# =========================================================

data(GFF)
# gff  = General Family Functioning
# hap  = Happiness
# sat  = Satisfaction with life
# AD   = Anxiety/Depression
# fc   = Family conflict
# qol  = Quality of life

mzData <- subset(GFF, zyg_2grp == "MZ")
dzData <- subset(GFF, zyg_2grp == "DZ")

selDVs <- c("gff", "fc", "qol", "hap", "sat", "AD")

#------------------------------------------------------------------------------
# Correlated Factors ACE Model [CF]
#------------------------------------------------------------------------------

# This is just an "atheoretical" ACE model:
CF <- umxACEv( selDVs=selDVs, sep="_T", dzData=dzData, mzData=mzData, addCI=F, tryHard="yes")
## Unfortunately, the path diagram is kind of messy.  
## Laying out a complicated path diagram in a way that "looks pretty" is hard to automate...

parameterSpecifications(CF)

#------------------------------------------------------------------------------
# 3-Factor Common-Pathway Model [CP3]
#------------------------------------------------------------------------------

# Let's make and run an MxModel object with `umxCP()`:
CP3 <- umxCP( "CP3", selDVs=selDVs, nFac=3, sep="_T", dzData=dzData, mzData=mzData, addCI=F, tryHard="yes", autoRun=T)
## Notice that variable "hap" has no unique variance!  
## This rather degenerate scenario is known as a "Haywood case".

# Run the MxModel
CP3 <- mxTryHard(CP3) 
## Sometimes optimization gets 'stuck'; check if TryHard reaches a better solution

# Plot a path diagram
umxPlotCP(CP3)

# Show which parameters are estimated (marked in square brackets)
parameterSpecifications(CP3)

#------------------------------------------------------------------------------
# 3-Factor Bifactor Common-Pathway Model [biCP3]
#------------------------------------------------------------------------------

# Let's change the CP3 model into a bifactor model:
# Copy the MxModel object to `biCP3`
biCP3 <- mxRename( CP3, "bifCP3")  

# We need to fix certain factor loadings to zero to specify a bifactor model:
# the first common factor will be a group factor for the first 3 phenotypes;
biCP3$top$cp_loadings$free[4:6,1] <- FALSE 
biCP3$top$cp_loadings$values[4:6,1] <- 0 

# the second common factor will be the general factor;
# the third common factor will be a group factor for the second 3 phenotypes;
biCP3$top$cp_loadings$free[1:3,3] <- FALSE
biCP3$top$cp_loadings$values[1:3,3] <- 0

# View the matrix of loadings
biCP3$top$cp_loadings  

# Run the MxModel
biCP3 <- mxTryHard(biCP3) 

# Plot a path diagram
umxPlotCP(biCP3)

#------------------------------------------------------------------------------
# Independent-Pathway Model
#------------------------------------------------------------------------------

# We will reuse CP3 as starting place for IP3, but give it a new "MxName":
IP3 <- mxRename(CP3,"IP3")

# We need to fix some common factors' loadings (onto the higher-order A, C, and E factors) to zero to specify an IP model:
# the first factor will represent common additive-genetic variance in the phenotypes 
# View the common factors' matrix of loadings onto the higher-order 'A' factors, before we modify it
IP3$top$a_cp 
# Modify it
IP3$top$a_cp$values <- diag(c(1,0,0))
IP3$top$a_cp$free <- diag(c(T,F,F))
# View after modification
IP3$top$a_cp

# the second factor will represent common shared-environmental variance;
IP3$top$c_cp$values <- diag(c(0,1,0))
IP3$top$c_cp$free <- diag(c(F,T,F))
# the third factor will represent common nonshared-environmental variance:
IP3$top$e_cp$values <- diag(c(0,0,1))
IP3$top$e_cp$free <- diag(c(F,F,T))

# Run the MxModel
IP3 <- mxTryHard(IP3)

# Plot a path diagram
umxPlotCP(IP3) 

# ============================
# Model comparison
# ============================

# The unstructured model is clearly the AIC-preferred model:
umxCompare(CF, c(CP3, biCP3, IP3))
omxAkaikeWeights(list(CF,CP3,biCP3,IP3))