### This is thew basic ordinal indicator latent factor model. 
### it is an adaptation of the script on the openMX website
### The website contains many nice example scripts, this one 
### i adapted to use the same syntax style as you are used to
### here at the boulder workshop.

###############################
### Read in the data        ###
###############################


alldata <- read.table("Boulder 2012 data.dat", header=T,sep="\t",na.string="-1")

## select the variables of intrest:
Ordinal <- alldata[,c(6:8,10,14)]

## Tell R and OpenMx that these data are ordered factors:
Ordinal <- mxFactor(Ordinal,levels=c(0:2))


## create variables indicating the number of variables, Factors and Thresholds:
nVariables<-5
nFactors<-1
nThresholds<-2


###############################
### Matrices and algebra    ###
###############################
loadings  <- mxMatrix(
  type="Full",
  nrow=nVariables,
  ncol=nFactors,
  free=T,
  values=.6,
  lbound=-.99,
  ubound=.99,
  name="facLoadings",
  labels=c("G1L1","G1L2","G1L3","G1L4","G1L5" )
  )

FacVariance <-    mxMatrix(
  type="Symm",
  nrow=nFactors,
  ncol=nFactors,
  free=F,
  values=1,
  name="facVar",
  labels="G1P"
  )


FacMean <-   mxMatrix(
  type="Full",
  nrow=nFactors,
  ncol=nFactors,
  free=F,
  values=0,
  name="facMeans",
  label="G1FacMean"
  )

UnitVector  <-  mxMatrix(
  type="Unit",
  nrow=nVariables,
  ncol=1,
  name="vectorofOnes"
  )

Residuals <-   mxAlgebra(
  expression=vectorofOnes - (diag2vec(facLoadings %*% facVar %*% t(facLoadings))) ,
  name="resVariances"
  )

ExpCovariance <-   mxAlgebra(
  expression=facLoadings %*% facVar %*% t(facLoadings) + vec2diag(resVariances),
  name="expCovariances"
  )

Intercepts <- mxMatrix(
  type="Full",
  nrow=1,
  ncol=nVariables,
  free=F,
  values=0,
  name="varMeans",
  label=c("G1M1","G1M2","G1M3","G1M4","G1M5")         
  )

ThresholdDeviations <- mxMatrix(
  type="Full",
  nrow=nThresholds,
  ncol=nVariables,
  free=TRUE,
  values=.2,
  lbound=rep( c(-Inf,rep(.01,(nThresholds-1))) , nVariables),
  dimnames=list(c(), varNames),
  name="thresholdDeviations",
  labels=letters[1:(nThresholds*nVariables)]
  )

UnitLower <- mxMatrix(
  type="Lower",
  nrow=nThresholds,
  ncol=nThresholds,
  free=FALSE,
  values=1,
  name="unitLower"
  )

ExpThresholds <-    mxAlgebra(
  expression=unitLower %*% thresholdDeviations,
  name="expThresholds"
  )

Data <- mxData(
  observed=MaleOrdinal,
  type='raw'
  )

ExpMeans <- mxAlgebra(expression= varMeans + t(facLoadings %*% facMeans),
                       name="expMeans" 
                       )

Objective <-   mxFIMLObjective(
  covariance="expCovariances",
  means="expMeans",
  dimnames=varNames,
  thresholds="expThresholds"
  )


######################################################
### Build the model from the Matrices and Algebras ###
######################################################

oneFactorModel <- mxModel("Male",loadings,FacVariance,FacMean,UnitVector,
                          Residuals,ExpCovariance, Intercepts,
                          ThresholdDeviations,UnitLower,
                          ExpThresholds,Data,ExpMeans,Objective)


###########
### Run ###
###########

ordinalModelFit <- mxRun(oneFactorModel)