# inst/models/passing/jointFactorModelsTest.R In OpenMx/OpenMx: Extended Structural Equation Modelling

```require(OpenMx)

data("jointdata", package ="OpenMx", verbose= TRUE)
jointData <- jointdata

# specify ordinal columns as ordered factors
jointData[,c(2,4,5)] <- mxFactor(jointData[,c(2,4,5)],
levels=list(c(0,1), c(0, 1, 2, 3), c(0, 1, 2)))

satCov <- mxMatrix("Symm", 5, 5,
free=TRUE, values=diag(5), name="C")
satCov\$free[2,2] <- FALSE
satCov\$free[4,4] <- FALSE
satCov\$free[5,5] <- FALSE

free=TRUE, values=1, name="L", lbound=0)

resid <- mxMatrix("Diag", 5, 5,
free=c(TRUE, FALSE, TRUE, FALSE, FALSE), values=.5, name="U")

means <- mxMatrix("Full", 1, 5,
free=c(TRUE, FALSE, TRUE, FALSE, FALSE), values=0, name="M")

thresh <- mxMatrix("Full", 3, 3, FALSE, 0, name="T")

thresh\$free[,1] <- c(TRUE, FALSE, FALSE)
thresh\$values[,1] <- c(0, NA, NA)
thresh\$labels[,1] <- c("z2t1", NA, NA)

thresh\$free[,2] <- TRUE
thresh\$values[,2] <- c(-1, 0, 1)
thresh\$labels[,2] <- c("z4t1", "z4t2", "z4t3")

thresh\$free[,3] <- c(TRUE, TRUE, FALSE)
thresh\$values[,3] <- c(-1, 1, NA)
thresh\$labels[,3] <- c("z5t1", "z5t2", NA)

omxCheckError(mxExpectationNormal("C", "M", dimnames=names(jointData),
thresholds="T", threshnames=c("z2", "z4", "z2")),
(paste("'threshnames' argument contains 'z2' more than once. \nIf you are having problems with Doppelgangers",
"perhaps you should check the basement for pods :)")))

jointData\$weight <- runif(nrow(jointData), min=.98,max=1.02)
weightedFits <- c()

for (strat in c('auto', 'ordinal', 'continuous')) {
print(paste('***',strat,'***'))
# run factor and saturated models
jointModel1 <- mxModel("ContinuousOrdinalData",
mxData(jointData, "raw"),
mxAlgebra(t(L) %*% L + U, name="C"),
mxFitFunctionML(jointConditionOn=strat),
mxExpectationNormal("C", "M",
dimnames=names(jointData)[1:5],
thresholds="T",
threshnames=c("z2", "z4", "z5"))
)

jointResults1 <- mxRun(jointModel1, suppressWarnings = TRUE)
#summary(jointResults1)
print(jointResults1\$fitfunction\$info)

#cat(deparse(round(coef(jointResults1),3)))
c1 <- c(0.609, 0.579, 0.657, 0.606, 0.165, 0.551, 0.499,
7.978, 2.069, 0.059, -0.387, 0.116, 0.815, -0.633, -0.285)
print(max(abs(coef(jointResults1) - c1)))
omxCheckCloseEnough(coef(jointResults1), c1, 0.005)

omxCheckCloseEnough(jointResults1\$output\$Minus2LogLikelihood, 2683.071, 0.2)

jointModel2 <- mxModel("ContinuousOrdinalData",
mxData(jointData, "raw"),
satCov, means, thresh,
mxFitFunctionML(jointConditionOn=strat),
mxExpectationNormal("C", "M",
dimnames=names(jointData)[1:5],
thresholds="T",
threshnames=c("z2", "z4", "z5"))
)

jointResults2 <- mxRun(jointModel2, suppressWarnings = TRUE)
#summary(jointResults2)

#cat(deparse(round(coef(jointResults2),3)))
c2 <- c(0.923, 0.374, 0.436, 0.367, 0.933, 0.354, 0.48, 0.43,  0.139,
0.266, 0.083, 0.201, 7.978, 2.07, 0.065, -0.413, 0.127,  0.874, -0.872, -0.394)
print(max(abs(coef(jointResults2) - c2)))
omxCheckCloseEnough(coef(jointResults2), c2, 0.005)

omxCheckCloseEnough(jointResults2\$output\$Minus2LogLikelihood, 2674.235, 0.1)

jointModel3 <- mxModel(jointResults2)
jointModel3\$data\$weight <- 'weight'
weightedFit <- mxRun(mxModel(jointModel3, mxComputeOnce('fitfunction','fit')))
weightedFits <- c(weightedFits, weightedFit\$output\$fit)
}

omxCheckCloseEnough(diff(weightedFits), rep(0,2), .04)
```
OpenMx/OpenMx documentation built on Dec. 5, 2019, 4:22 a.m.