family link functions v. threshold model on polychoric

raw <- sample[, c(vars1, vars2)]
 
raw$item1_ordinal <- mxFactor(raw$item1_ordinal, levels = c(0,1,2,3))
raw$item2_ordinal <- mxFactor(raw$item2_ordinal, levels = c(0,1,2,3))
raw$item3_ordinal <- mxFactor(raw$item3_ordinal, levels = c(0,1,2,3))
raw$item4_ordinal <- mxFactor(raw$item4_ordinal, levels = c(0,1,2,3))
raw$item5_ordinal <- mxFactor(raw$item5_ordinal, levels = c(0,1,2,3))
raw$item6_ordinal <- mxFactor(raw$item6_ordinal, levels = c(0,1,2,3))
raw$item7_ordinal <- mxFactor(raw$item7_ordinal, levels = c(0,1,2,3))
raw$item8_ordinal <- mxFactor(raw$item8_ordinal, levels = c(0,1,2,3))
 
dataRaw <- mxData(observed = raw, type="raw")
 
# residual variances
resVars <- mxPath(from = c(vars1, vars2), arrows = 2, # how about 8?
                  free = TRUE, values = c(1,1,1,1,1,1,1,1),
                  labels = c("e1", "e2", "e3", "e4", "e5", "e6", "e7", "e8"))
 
# latent variance and covariance
latVars <- mxPath(from = c("F1", "F2"), arrows = 2, connect = "unique.pairs", # is "unique.pairs" correct?
                  free = TRUE, values = c(1, .5, 1), labels = c("varF1", "cov", "varF2"))
 
# factor loadings
facLoadsX <- mxPath(from="F1", to=vars1, arrows=1,
                   free=c(FALSE,TRUE,TRUE,TRUE), values=c(1,1,1,1), labels = c("11","12","13","14"))
facLoadsY <- mxPath(from="F2", to=vars2, arrows=1,
                    free=c(FALSE,TRUE,TRUE,TRUE), values=c(1,1,1,1), labels = c("15","16","17","18"))
# means
means <- mxPath(from = "one", to = c(vars1,vars2,"F1","F2"),
                arrows = 1, 
                free = c(T,T,T,T,T,T,T,T,F,F),
                values = c(1,1,1,1,1,1,1,1,0,0),
                labels = c("meanitem1","meanitem2","meanitem3","meanitem4","meanitem5","meanitem6","meanitem7", "meanitem8", NA, NA))
  # We also need means and intercepts in our model. Exogenous or independent variables have means, while endogenous
  # or dependent variables have intercepts. These can be included by regressing both x and y on a constant, which can
  # be refered to in OpenMx by "one".   
 
# thresholds
mxThreshold(vars = c(vars1, vars2), nThresh = c(3,3,3,3,3,3,3,3), free = TRUE, values = c(0,1,2,0,1,2,0,1,2,0,1,2,0,1,2,0,1,2,0,1,2,0,1,2)) # four ordinal values implies: three thresholds and three values (for each variable) 
  # It is also important to remember that specifying thresholds is not sufficient to get an ordinal data model to run. In
  # addition, the scale of each ordinal variable must be identified just like the scale of a latent variable. The most common
  # method for this involves constraining a ordinal item’s mean to zero and either its total or residual variance to a constant
  # value (i.e., one). For variables with two or more thresholds, ordinal variables may also be identified by constraining
  # two thresholds to fixed values
 
twoFactormodel <- mxmodel("Two Factor Model Path Specification", type="RAM", # is that the type?
                          manifestVars=c(vars1, vars2), latentVars=c("F1", "F2"),
                          dataRaw, resVars, latVars, facLoadsX, facLoadsY, means, thresholds)
  # When using a path specification of the model, the fit function is always RAM which is indicated by using the type
  # argument. We don’t have to specify the fit function explicitly with an mxExpectation() and FitFunction()
  # argument, instead we simply add the following argument to the model.
 
twoFactorResults <- mxRun(twoFactormodel)
 
summary(twoFactorResults)
output(twoFactorResults)