# ---------------------------------------------------------------------
# Program: ThreeFactorScale2Test.R
#  Author: Steven M. Boker
#    Date: Sun Mar 14 14:42:16 EDT 2010
#
# This program tests the number of factors in simulated data.
#    using a standard RAM.
#
# ---------------------------------------------------------------------
# Revision History
#    -- Sun Mar 14 14:42:18 EDT 2010
#      Created ThreeFactorScale2Test.R.
#
# ---------------------------------------------------------------------

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

options(width=100)

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

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

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

describe(factorScaleExample2)

numberFactors <- 3
indicators <- names(factorScaleExample2)
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="") 

threeFactorOrthoRaw1 <- mxModel("threeFactorOrthogonal",
    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=indicators, 
           arrows=2, 
           free=TRUE, values=.8, 
           labels=uniqueLabels),
    mxPath(from=latents,
           arrows=2, 
           free=FALSE, values=1, 
           labels=factorVarLabels),
    mxPath(from="one", to=indicators, 
           arrows=1, free=TRUE, values=.1, 
           labels=meanLabels),
    mxData(observed=factorScaleExample2, type="raw")
    )

threeFactorOrthoRaw1Out <- mxRun(threeFactorOrthoRaw1)

summary(threeFactorOrthoRaw1Out)

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

threeFactorObliqueRaw1 <- mxModel(threeFactorOrthoRaw1,
    mxPath(from=latents,to=latents,all=TRUE,
           arrows=2, 
           free=TRUE, values=.3),
    mxPath(from=latents,
           arrows=2, 
           free=FALSE, values=1, 
           labels=factorVarLabels),
    name="threeFactorOblique"
    )

threeFactorObliqueRaw1Out <- mxRun(threeFactorObliqueRaw1)

summary(threeFactorObliqueRaw1Out)

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

threeFactorObliqueRaw2 <- mxModel(threeFactorOrthoRaw1,
    mxPath(from="F2",to="F3",
           arrows=2, 
           free=FALSE, values=1),
    mxPath(from="F1",to="F3",
           arrows=2, 
           free=TRUE, values=.3,labels="C1"),
    mxPath(from="F1",to="F2",
           arrows=2, 
           free=TRUE, values=.3,labels="C1"),
    mxPath(from=latents,
           arrows=2, 
           free=FALSE, values=1, 
           labels=factorVarLabels),
    name="threeFactorOblique2"
    )

threeFactorObliqueRaw2Out <- mxRun(threeFactorObliqueRaw2)

summary(threeFactorObliqueRaw2Out)

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

threeFactorObliqueRaw3 <- mxModel(threeFactorOrthoRaw1,
    mxPath(from=latents, to=latents,
           arrows=2, all=TRUE,
           free=FALSE, values=1, 
           labels=factorVarLabels),
    name="threeFactorAllOne"
    )

threeFactorObliqueRaw3Out <- mxRun(threeFactorObliqueRaw3)

summary(threeFactorObliqueRaw3Out)

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


tableFitStatistics(threeFactorObliqueRaw1Out,
                   c(threeFactorObliqueRaw2Out,threeFactorOrthoRaw1Out,threeFactorObliqueRaw3Out))