# ---------------------------------------------------------------------
# Program: ThreeLatentMediationTest2.R
#  Author: Steven M. Boker
#    Date: Sun Mar 14 16:28:34 EDT 2010
#
# This program tests variations on a latent mediation model
#    using a standard RAM.
#
# ---------------------------------------------------------------------
# Revision History
#    -- Sun Mar 14 16:28:38 EDT 2010
#      Created ThreeLatentMediationTest2.R.
#
# ---------------------------------------------------------------------

# ----------------------------------
# Read libraries and set options.

options(width=100)

require(psych)
require(OpenMx)
source("GenEpiHelperFunctions.R")

# ----------------------------------
# Read the data and print descriptive statistics.

latentMultipleRegExample2 <- read.csv("latentMultipleRegExample2.csv")

describe(latentMultipleRegExample2)

numberFactors <- 3
indicators <- names(latentMultipleRegExample2)
numberIndicators <- length(indicators)
totalVars <- numberIndicators + numberFactors

# ----------------------------------
# Build an Old-style RAM OpenMx single factor FIML model with fixed variance

latents <- paste("F", 1:numberFactors, sep="")
loadingLabels <- paste("b_F", rep(1:numberFactors, each=numberIndicators), 
                              rep(indicators, numberFactors), sep="") 
loadingLabels

uniqueLabels <- paste("U_", indicators, sep="")
meanLabels <- paste("M_", indicators, sep="")
factorVarLabels <- paste("Var_", latents, sep="")

latents1 <- latents[1]
indicators1 <- indicators[1:4]
loadingLabels1 <- paste("b_F1", indicators[1:4], sep="") 
latents2 <- latents[2]
indicators2 <- indicators[5:8]
loadingLabels2 <- paste("b_F2", indicators[5:8], sep="") 
latents3 <- latents[3]
indicators3 <- indicators[9:12]
loadingLabels3 <- paste("b_F3", indicators[9:12], sep="") 

threeLatentOrthoRaw1 <- mxModel("threeLatentOrthogonal",
    type="RAM",
    manifestVars=indicators,
    latentVars=latents,
    mxPath(from=latents1, to=indicators1, 
           arrows=1, all=TRUE, 
           free=TRUE, values=.2, 
           labels=loadingLabels1),
    mxPath(from=latents2, to=indicators2, 
           arrows=1, all=TRUE, 
           free=TRUE, values=.2, 
           labels=loadingLabels2),
    mxPath(from=latents3, to=indicators3, 
           arrows=1, all=TRUE, 
           free=TRUE, values=.2, 
           labels=loadingLabels3),
    mxPath(from=latents1, to=indicators1[1], 
           arrows=1, 
           free=FALSE, values=1),
    mxPath(from=latents2, to=indicators2[1], 
           arrows=1, 
           free=FALSE, values=1),
    mxPath(from=latents3, to=indicators3[1], 
           arrows=1, 
           free=FALSE, values=1),
    mxPath(from=indicators, 
           arrows=2, 
           free=TRUE, values=.8, 
           labels=uniqueLabels),
    mxPath(from=latents,
           arrows=2, 
           free=TRUE, values=.8, 
           labels=factorVarLabels),
    mxPath(from="one", to=indicators, 
           arrows=1, free=TRUE, values=.1, 
           labels=meanLabels),
    mxData(observed=latentMultipleRegExample2, type="raw")
    )

threeLatentOrthoRaw1Out <- mxRun(threeLatentOrthoRaw1)

summary(threeLatentOrthoRaw1Out)

round(mxEval(A[1:numberIndicators,(numberIndicators+1):totalVars], 
             threeLatentOrthoRaw1Out), 3)

threeLatentMediation1 <- mxModel(threeLatentOrthoRaw1,
    mxPath(from="F1",to="F3",
           arrows=1, 
           free=TRUE, values=.2,
           labels="b13"),
    name="threeLatentMediation1"
    )

threeLatentMediation1Out <- mxRun(threeLatentMediation1)

summary(threeLatentMediation1Out)

round(mxEval(A[1:numberIndicators,(numberIndicators+1):totalVars], threeLatentMediation1Out), 3)

threeLatentMediation2 <- mxModel(threeLatentMediation1,
    mxPath(from="F2",to="F3",
           arrows=1, 
           free=TRUE, values=.2,
           labels="b23"),
    name="threeLatentMediation2"
    )

threeLatentMediation2Out <- mxRun(threeLatentMediation2)

summary(threeLatentMediation2Out)

round(mxEval(A[1:numberIndicators,(numberIndicators+1):totalVars], 
             threeLatentMediation2Out), 3)

threeLatentMediation3 <- mxModel(threeLatentMediation2,
    mxPath(from="F1",to="F2",
           arrows=1, 
           free=TRUE, values=.2,
           labels="b12"),
    name="threeLatentMediation3"
    )

threeLatentMediation3Out <- mxRun(threeLatentMediation3)

summary(threeLatentMediation3Out)

round(mxEval(A[1:numberIndicators,(numberIndicators+1):totalVars], 
             threeLatentMediation3Out), 3)


threeLatentMediation4 <- mxModel(threeLatentMediation3,
    mxPath(from="F1",to="F3",
           arrows=1, 
           free=FALSE, values=0),
    name="threeLatentMediation4"
    )

threeLatentMediation4Out <- mxRun(threeLatentMediation4)

summary(threeLatentMediation4Out)

round(mxEval(A[1:numberIndicators,(numberIndicators+1):totalVars], 
             threeLatentMediation4Out), 3)


tableFitStatistics(threeLatentMediation3Out,
                   c(threeLatentMediation2Out,
                     threeLatentMediation4Out,
                     threeLatentMediation1Out,
                     threeLatentOrthoRaw1Out))