#-------------------------------------------------------------------------------
# Author: Michael D. Hunter
# Date: Mon Mar  4 01:05:36 EST 2024
# Filename: IntroSEMOpenMx.R
# Purpose: Show some basic introductory material for structural equation models
#  and OpenMx
#-------------------------------------------------------------------------------


#-------------------------------------------------------------------------------
# Set working directory


#-------------------------------------------------------------------------------
# Load packages

require(OpenMx)


#-------------------------------------------------------------------------------
# Exercise 1
# Write paths that go with this diagram

# A very literal way of doing this
mxPath(from='one', to='g', arrows=1)
mxPath(from='one', to='x', arrows=1)
mxPath(from='g', to='x', arrows=1)
mxPath(from='x', to='x', arrows=2)

# But arrows=1 is the default
# And you can combine multiple path into a single call
# And you don't need to name all the arguments if they are in the order
#  OpenMx expects
# The below is equivalent to the above
p1 <- mxPath('one', c('g', 'x'))
p2 <- mxPath('g', 'x')
p3 <- mxPath('x', 'x', arrows=2)

mpaths <- list(p1, p2, p3)

# If you want, you can super-duper combine paths
#  and change the values
mikespath <- mxPath(c('one', 'one', 'g', 'x'),
       c('g', 'x', 'x', 'x'),
       arrows=c(1, 1, 1, 2),
       values=c(.5, -4, .8, 2))

mikespath



#-------------------------------------------------------------------------------
# Exercise 2
# Draw a diagram of this model

ex2 <- mxModel('Exercise 2',
               type='RAM', manifestVars='Mike', latentVars='Brad',
               mxPath(from='one', to='Mike', arrows=1, values=1.4, labels='m'),
               mxPath(from='Mike', arrows=2, values=.5, labels='v'),
               mxPath(from='Mike', to='Brad', arrows=1, values=2, labels='r')
)

ex2
ex2$A

require(semPlot)
semPaths(ex2)


#-------------------------------------------------------------------------------
# Exercise 3
# Write a matrix that represents this diagram

v <- letters[24:26]
ex3 <- matrix(c(0, 0, .6, 0, 0, 0, 0, -.3, 0), 3, 3, dimnames=list(v, v))
ex3

# Or do it as an mxMatrix
mex3 <- mxMatrix(nrow=3, ncol=3, values=ex3, name='Coolio')
mex3
# Notice how it kept the dimnames

# Change the matrix so only the non-zero values are free
mex3$free <- (ex3 != 0)
mex3


#-------------------------------------------------------------------------------
# Exercise 4
# Write a matrix that represents this diagram

w <- letters[6:8]
ex4 <- matrix(c(0, 0, 'b', 0, 'a', 0, 'b', 0, 0), 3, 3, dimnames=list(w, w))
ex4

# Do it as an mxMatrix
# First args are type, nrow, and ncol
mex4 <- mxMatrix('Symm', 3, 3, labels=c(NA, NA, 'b', 'a', NA, NA),
                 dimnames=list(w, w))
mex4
# Notice that for the symmetric matrix, you only *have* to write the
#  lower triangular elements.
# But you *can* write all of them if you want.


#-------------------------------------------------------------------------------
# Exercise 5
# Write a matrix that represents this diagram

mex4
# Bit of a trick question.
# It's the same matrix.


#-------------------------------------------------------------------------------
# Exercise 6
# Write a model that represents this diagram

mvars <- c('w', 'x', 'y', 'z')
ex6 <- mxModel('Exercise 6',
               type='RAM',
               latentVars='g',
               manifestVars=mvars,
               mxPath('g', mvars, values=runif(length(mvars))))
# Look at it
ex6
# Poke around
ex6$A
ex6$F

# Use mxEval
mxEval(A, ex6)
mxEval(S, ex6)
mxEval(F, ex6)


#-------------------------------------------------------------------------------
# Exercise 7
# Write a model that represents this diagram

ex6
# Again, bit of a trick question.
# It's the same model.


#-------------------------------------------------------------------------------
# Exercise 8
# Write a matrix for the covariances of A1 & A2.
# Write a matrix for the covariances of C1 & C2.
# Write a matrix for the covariances of E1 & E2.

acov <- mxMatrix('Symm', 2, 2, values=c(1, .5, 1), name='A')
ccov <- mxMatrix('Unit', 2, 2, name='C')
# Could also be
# ccov <- mxMatrix('Symm', 2, 2, values=c(1, 1, 1), name='C')
ecov <- mxMatrix('Iden', 2, 2, name='E')
# Could also be
# ecov <- mxMatrix('Symm', 2, 2, values=c(1, 0, 1), name='E')

acov
ccov
ecov

# Super bonus
# Make this into a model for the covariance
ex8 <- mxModel('Exercise 8 Model',
               acov, ccov, ecov,
               mxMatrix('Full', 1, 3, free=TRUE, values=c(1.2, .5, .2),
                        labels=c('a2', 'c2', 'e2'), name='theta'),
               mxAlgebra(a2 %x% A + c2 %x% C + e2 %x% E, name='Yoshi'),
               mxExpectationNormal(covariance='Yoshi', dimnames=c('p1', 'p2')))
ex8$Yoshi
mxEval(Yoshi, ex8)
mxEval(Yoshi, ex8, compute=TRUE)
1.2 * mxEval(A, ex8)
0.5 * mxEval(C, ex8)
0.2 * mxEval(E, ex8)


#-------------------------------------------------------------------------------
# End