require(OpenMx)
require(psych)

# Read in an example dataset supplied with the OpenMx package
data(myFADataRaw)
# Grab three variables
df <- myFADataRaw[, c("x1","z1", "z3")]

# rename them to be more memorable
manifests = c("cont1", "ord1","ord2")
names(df) <-manifests 

# For OpenMx to work with ordinal data, they must be ordered factors!
df$ord1 <- mxFactor(df$ord1, levels=0:1)
df$ord2 <- mxFactor(df$ord2, levels=0:2)

# Let's see what we're working with now
str(df)
'data.frame':	500 obs. of  3 variables:
 $ cont1: num  1 1.63 1.95 2.95 3.66 ...
 $ ord1 : Ord.factor w/ 2 levels "0"<"1": 2 2 2 2 2 2 2 2 2 1 ...
 $ ord2 : Ord.factor w/ 3 levels "0"<"1"<"2": 1 2 2 1 2 2 1 2 1 1 ...

hist(df$cont1) # nice and normal
plot(~df$ord1) # binary - 4x as many 1s as 0s
plot(~df$ord2) # three-values: low endorsement of the highest value

# Just to help keep things simple in the model, let's list up the continuous 
# and the ordinal manifests separately
ordinalVars = manifests[2:3]
continuousVars = manifests[1]

m1 <- mxModel("joint CFA", type="RAM",
	manifestVars = manifests,
	latentVars = "F1",
	# Factor loadings
	mxPath(from="F1", to = manifests, arrows=1, free=T, values=1, labels=c("l1","l2","l3")),

	# Manifest continuous mean & residual variance (both free)
	mxPath(from="one", to = continuousVars, arrows=1, free=T, values=0, labels="meancont1"),
	mxPath(from=continuousVars            , arrows=2, free=T, values=1, labels="e1" ),

	# Ordinal manifests mean and variance fixed at 0 (the threshold matrix will esimate these)
	mxPath(from="one", to = ordinalVars, arrows=1, free=F, values=0, labels=c("meanord1","meanord2")),
	mxPath(from=ordinalVars            , arrows=2, free=F, values=1, labels=c("e2","e3")),
	
	# fix latent variable variance @1 & mean @0 
	mxPath(from="F1"          , arrows=2, free=F, values=1, labels ="varF1"),
	mxPath(from="one", to="F1", arrows=1, free=F, values=0, labels="meanF1"),

	# Make a matrix as wide as the number of ordinal variables (2), and as many rows as the max # of thresholds (3) - 1
	# In this case that is 2 * 2
    # leave it free where there is a real threshold to be set, and fill the rest of each column fixed @NA
	# set the column names to the name of the ordinal variables: that is how OpenMx maps into this threshold matrix
	mxMatrix("Full", nrow=2, ncol=2, byrow=T, name="thresh",
		dimnames = list(c(), ordinalVars),
		free  = c(T, T ,
			      F, T), 
		values = c(0, -1,
			      NA,  1)
	),
	# Tell the mxRAMObjective that we have not only A,S, F and M matrices, but also some ordinal vars with thresholds
	 mxRAMObjective(A="A", S="S", F="F", M="M", thresholds="thresh"),
 	# And last but not least, the data
 	mxData(df, type="raw")
)

# run!
m1 = mxRun(m1); summary(m1)

# free parameters:
#        name matrix   row   col   Estimate  Std.Error lbound ubound
# 1        l1      A cont1    F1 -0.1953311 0.04115263              
# 2        l2      A  ord1    F1 -1.0135426        NaN              
# 3        l3      A  ord2    F1 -4.3049347        NaN              
# 4        e1      S cont1 cont1  0.9571266 0.06018277              
# 5 meancont1      M     1 cont1  2.9876221 0.04461568              
# 6      <NA> thresh     1  ord1 -1.3038918        NaN              
# 7      <NA> thresh     1  ord2 -2.1373824        NaN              
# 8      <NA> thresh     2  ord2  5.7570834        NaN              
# 
# observed statistics:  1500 
# estimated parameters:  8 
# degrees of freedom:  1492 
# -2 log likelihood:  2678.115 
# saturated -2 log likelihood:  NA 
# number of observations:  500 
# chi-square:  NA 
# p:  NA 
# Information Criteria: 
#     df Penalty Parameters Penalty Sample-Size Adjusted
# AIC  -305.8846           2694.115                   NA
# BIC -6594.0798           2727.832              2702.44
# CFI: NA 
# TLI: NA 
# RMSEA:  NA