require(OpenMx)

### we read in the data.

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

## Split the data by sex:
MaleOrdinal <- alldata[alldata[,4]==1,c(6:8,10,14)]
FemaleOrdinal <- alldata[alldata[,4]==2,c(6:8,10,14)]

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

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

## create variable names for use in the OpenMx Objective function:

varNames <- paste("cbcl",1:nVariables,sep="")
names(MaleOrdinal) <- varNames
names(FemaleOrdinal) <-varNames


###################################
### Male Matrices and algebra   ###
###################################
Mloadings  <- mxMatrix(
        type="Full",
        nrow=nVariables,
        ncol=nFactors,
        free=T,
        values=.6,
        lbound=-.99,
        ubound=.99,
        name="facLoadings",
        labels=c("G1L1","G1L2","G1L3","G1L4","G1L5" )
   		 )
   		 
MFacVariance <-    mxMatrix(
    	   type="Symm",
        nrow=nFactors,
        ncol=nFactors,
        free=F,
        values=1,
        name="facVar",
        labels="G1P"
    	)
    	
    	
 MFacMean <-   mxMatrix(
    	type="Full",
    	nrow=nFactors,
    	ncol=nFactors,
    	free=F,
    	values=0,
    	name="facMeans",
    	label="G1FacMean"
    	)
    	
MUnitVector  <-  mxMatrix(
        type="Unit",
        nrow=nVariables,
        ncol=1,
        name="vectorofOnes"
    )
    
 MResiduals <-   mxAlgebra(
        expression=vectorofOnes - (diag2vec(facLoadings %*% facVar %*% t(facLoadings))) ,
        name="resVariances"
    )
    
 MExpCovariance <-   mxAlgebra(
        expression=facLoadings %*% facVar %*% t(facLoadings) + vec2diag(resVariances),
        name="expCovariances"
    )
    
 MIntercepts <-     mxMatrix(
        type="Full",
        nrow=1,
        ncol=nVariables,
        free=F,
		values=0,
        name="varMeans",
        label=c("G1M1","G1M2","G1M3","G1M4","G1M5")         
    )
    
MThresholdDeviations <- 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)]
    )
    
MUnitLower <- mxMatrix(
        type="Lower",
        nrow=nThresholds,
        ncol=nThresholds,
        free=FALSE,
        values=1,
        name="unitLower"
    )
    
MExpThresholds <-    mxAlgebra(
        expression=unitLower %*% thresholdDeviations,
        name="expThresholds"
    )
    
MData <- mxData(
        observed=MaleOrdinal,
        type='raw'
    )
    
MExpMeans <- mxAlgebra(expression= varMeans + t(facLoadings %*% facMeans),
        name="expMeans" 
        )

MObjective <-   mxFIMLObjective(
        covariance="expCovariances",
        means="expMeans",
        dimnames=varNames,
        thresholds="expThresholds"
        )
        
        
###################################
### Female Matrices and algebra ###
###################################


FLoadings <- mxMatrix(
        type="Full",
        nrow=nVariables,
        ncol=nFactors,
        free=T,
        values=.6,
        lbound=-.99,
        ubound=.99,
        name="facLoadings",
        labels=c("G2L1","G2L2","G2L3","G2L4","G2L5" )
    )
FFacVariance <-    mxMatrix(
    	   type="Symm",
        nrow=nFactors,
        ncol=nFactors,
        free=F,
        values=1,
        name="facVar",
        labels="G2P"
    )
FFacMean <-    mxMatrix(
    	type="Full",
    	nrow=nFactors,
    	ncol=nFactors,
    	free=F,
    	values=0,
    	name="facMeans",
    	label="G2FacMean"
    	)
FUnitVector <-    mxMatrix(
        type="Unit",
        nrow=nVariables,
        ncol=1,
        name="vectorofOnes"
    )
FResiduals <-   mxAlgebra(
        expression=vectorofOnes - (diag2vec(facLoadings %*% facVar %*% t(facLoadings))) ,
        name="resVariances"
    )
FExpCovariance <-    mxAlgebra(
        expression=facLoadings %*% facVar %*% t(facLoadings) + vec2diag(resVariances),
        name="expCovariances"
    )
FIntercepts <-   mxMatrix(
        type="Full",
        nrow=1,
        ncol=nVariables,
        free=c(T,T,T,T,T),
		values=0,
        name="varMeans",
        label=c("G2M1","G2M2","G2M3","G2M4","G2M5")               
    )
FThesholdDeviations <-   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)]
    )
FUnitLower <-    mxMatrix(
        type="Lower",
        nrow=nThresholds,
        ncol=nThresholds,
        free=FALSE,
        values=1,
        name="unitLower"
    )
FExpThresholds <-   mxAlgebra(
        expression=unitLower %*% thresholdDeviations,
        name="expThresholds"
    )
FData <-    mxData(
        observed=FemaleOrdinal,
        type='raw'
    )
FExpMeans <-     mxAlgebra(expression= varMeans + t(facLoadings %*% facMeans),
        name="expMeans" )

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


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

maleModel <- mxModel("Male",Mloadings,MFacVariance,MFacMean,MUnitVector,MResiduals,MExpCovariance, MIntercepts,MThresholdDeviations,MUnitLower,MExpThresholds,MData,MExpMeans,MObjective)

femaleModel <- mxModel("Female",FLoadings,FFacVariance,FFacMean,FUnitVector,FResiduals,FExpCovariance, FIntercepts,FThesholdDeviations,FUnitLower,FExpThresholds,FData,FExpMeans,FObjective)

M2LL <-  mxAlgebra(
        Male.objective + Female.objective,
        name="minus2loglikelihood"
    )
Alg <-   mxAlgebraObjective("minus2loglikelihood")

#######################################################
#### Build and run the configural invariance model: ###
#######################################################

combinedModel <- mxModel("configural invariance" ,maleModel,femaleModel,M2LL,Alg)
ConfiguralFit <- mxRun(combinedModel)

#######################################################
#### Build and run the metric invariance model:     ###
#######################################################


### remember the metric invariance model requires that the factor loadings 
### are equal over the 2 groups. We will build a Loadings matrix that can be 
### put in both the Male and the Female model. these loadings will be invariant
### over sex. it will be called ILoadings:
ILoadings  <- mxMatrix(
  type="Full",
  nrow=nVariables,
  ncol=nFactors,
  free=T,
  values=.5,
  lbound=-.99,
  ubound=.99,
  name="facLoadings",
  labels=c("L1","L2","L3","L4","L5" )
  )

### Given restrictions on the loadings we can now free the factor variance in the second group:

FFacVariance <- mxMatrix(type="Symm",nrow=nFactors,ncol=nFactors, free=T, values=1, name="facVar",labels="G2P")


#### Build the Male and Female model, but now with the restricted Loadings:


maleModel <- mxModel("Male",ILoadings,MFacVariance,MFacMean,MUnitVector,MResiduals,MExpCovariance, MIntercepts,MThresholdDeviations,MUnitLower,MExpThresholds,MData,MExpMeans,MObjective)

femaleModel <- mxModel("Female",ILoadings,FFacVariance,FFacMean,FUnitVector,FResiduals,FExpCovariance, FIntercepts,FThesholdDeviations,FUnitLower,FExpThresholds,FData,FExpMeans,FObjective)

M2LL <-  mxAlgebra(
  Male.objective + Female.objective,
  name="minus2loglikelihood"
  )
Alg <-   mxAlgebraObjective("minus2loglikelihood")


#########################################################
#### Combine and run the metric invariance model:     ###
#########################################################


combinedModel <- mxModel("metric invariance" ,maleModel,femaleModel,M2LL,Alg)
metricFit <- mxRun(combinedModel)


########################################################################################
### our next step towards measurement invariance involves restricting the intercepts ###
### and freeing the factor mean. This is done because we need any difference between ###
### the sex groups to be only due to the factor   .                                 ###
########################################################################################

### We fix the female intercepts:
FIntercepts <- mxMatrix(
  type="Full",
  nrow=1,
  ncol=nVariables,
  free=F,
  values=0,
  name="varMeans",
  label=c("G2M1","G2M2","G2M3","G2M4","G2M5")               
  )

### We free the factor mean:
FFacMean <-mxMatrix(
  type="Full",
  nrow=nFactors,
  ncol=nFactors,
  free=T,
  values=0,
  name="facMeans",
  label="G2FacMean"
  )

### We rebuild the the two subgroup models withg the restricted means:
maleModel <- mxModel("Male",ILoadings,MFacVariance,MFacMean,MUnitVector,MResiduals,MExpCovariance, MIntercepts,MThresholdDeviations,MUnitLower,MExpThresholds,MData,MExpMeans,MObjective)

femaleModel <- mxModel("Female",ILoadings,FFacVariance,FFacMean,FUnitVector,FResiduals,FExpCovariance, FIntercepts,FThesholdDeviations,FUnitLower,FExpThresholds,FData,FExpMeans,FObjective)

M2LL <-  mxAlgebra(
  Male.objective + Female.objective,
  name="minus2loglikelihood"
  )
Alg <-   mxAlgebraObjective("minus2loglikelihood")

#########################################################
#### Combine and run the strong invariance model:     ###
#########################################################


combinedModel <- mxModel("strong invariance" ,maleModel,femaleModel,M2LL,Alg)
strongFit <- mxRun(combinedModel)

########################################################################################
### the last parameter that must be restricted to be equal over groups are the       ###
### residuals. These however are not free parameters in the models above they are    ###
### fully dependend on the labda's as you van see in this piece of code:             ###
###                                                                                  ###
### Residuals <-   mxAlgebra(                                                        ###
### expression=vectorofOnes - (diag2vec(facLoadings %*% facVar %*% t(facLoadings))), ###
### name="resVariances")                                                             ###
###                                                                                  ### 
### We cannot make the residuals free parameters, we can however fix them to         ###
### a value. This results in strict measurement invariance all differences between   ###
### groups are a function of differences in the common factor                        ###
###                                                                                  ###
### we first create a model in which the residuals are free in the second group      ###
### then we test wheter these free residuals van be resticted to be equal            ###
####to the first group                                                               ###
###                                                                                  ###
########################################################################################

####
#### step one is a restriction of the residuals in 1 group, this will allow us
#### to estimate them in the female group. This is strong invariance, but then with 
#### free residuals in the second female group:


IResiduals <- mxMatrix(type="Diag",
                       nrow=nVariables,
                       ncol=nVariables,
                       free=FALSE,
                       values=.4,
                       name="resVariances"
                       )
FResiduals <-  mxMatrix(type="Diag",
                        nrow=nVariables,
                        ncol=nVariables,
                        free=TRUE,
                        values=.4,
                        name="resVariances"
                        )
### Slightly tweak the expected covariance matrix because the residuals are no longer
### a vector:
MExpCovariance <-    mxAlgebra(
  expression=facLoadings %*% facVar %*% t(facLoadings) + (resVariances),
  name="expCovariances"
  )

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


### We rebuild the the two subgroup models withg the fixed residuals in males, 
### free in females:
maleModel <- mxModel("Male",ILoadings,MFacVariance,MFacMean,MUnitVector,
                     IResiduals,MExpCovariance, MIntercepts,MThresholdDeviations,
                     MUnitLower,MExpThresholds,MData,MExpMeans,MObjective)

femaleModel <- mxModel("Female",ILoadings,FFacVariance,FFacMean,FUnitVector,
                       FResiduals,FExpCovariance, FIntercepts,FThesholdDeviations,
                       FUnitLower,FExpThresholds,FData,FExpMeans,FObjective)

M2LL <-  mxAlgebra(
  Male.objective + Female.objective,
  name="minus2loglikelihood"
  )
Alg <-   mxAlgebraObjective("minus2loglikelihood")


#########################################################
#### Combine and run the strong invariance model:     ###
#########################################################


combinedModel <- mxModel("strong invariance 2" ,maleModel,femaleModel,M2LL,Alg)
strongFit2 <- mxRun(combinedModel)



################################################
### Last step! strict measurement invariance ###
################################################


### We rebuild the the two subgroup models withg the fixed residuals:
maleModel <- mxModel("Male",ILoadings,MFacVariance,MFacMean,MUnitVector,IResiduals,MExpCovariance, MIntercepts,MThresholdDeviations,MUnitLower,MExpThresholds,MData,MExpMeans,MObjective)

femaleModel <- mxModel("Female",ILoadings,FFacVariance,FFacMean,FUnitVector,IResiduals,FExpCovariance, FIntercepts,FThesholdDeviations,FUnitLower,FExpThresholds,FData,FExpMeans,FObjective)

M2LL <-  mxAlgebra(
  Male.objective + Female.objective,
  name="minus2loglikelihood"
  )
Alg <-   mxAlgebraObjective("minus2loglikelihood")


#########################################################
#### Combine and run the strict invariance model:     ###
#########################################################


combinedModel <- mxModel("strict invariance" ,maleModel,femaleModel,M2LL,Alg)
strictFit <- mxRun(combinedModel)