require(OpenMx)
require(MASS)
dataset <- read.table("Measurement_invariance_data.dat",header=T,sep="\t" )
source("GenEpiHelperFunctions.r")



#########################
#### Descriptives    ####
#########################

correlations <- cor(dataset[,2:5])
means <- colMeans(dataset[,2:5])

correlations
means

#########################
#### Saturated Model ####
#########################

# Create starting values for a covariance matrix, 
# this can be handy as the starting values should be "positive definite"
# so all 1's won't work in this case, the covariance as estimated by R will work
startcov <- cov(as.matrix(dataset[,2:5]))
# read in data:
dataSat <-  mxData( observed=dataset[,2:5], type="raw" )
#define a free covariance matrix:
covSat <- mxMatrix(type="Symm",nrow=4,ncol=4,free=T,values=startcov,
                   name="cov")
#define a free means matrix:
meanSat <- mxMatrix(type="Full",nrow=1,ncol=4,labels=c("m1","m2","m3","m4"),
                    values=10,free=T,name="M")
# create a objective, tell openMX how to match de data to the parameters:
objSat <- mxFIMLObjective( covariance="cov", means="M", dimnames=names(dataset[,2:5]))
#  combine everything into the model:
modelSat <- mxModel("Saturated Model",dataSat,covSat,meanSat,objSat)
# fit the model:
fitSat <- mxRun(modelSat)


#########################
#### 1-Factor Model  ####
#########################

# read in the data:
data <-  mxData( observed=dataset[,2:5], type="raw" )
# Matrix holding the factor loadings (Labda's in greek/algebra)
loadings <-  mxMatrix( type="Full", nrow=4, ncol=1, values=1, free=c(F,T,T,T),
                       labels=c("l1","l2","l3","l4"),
                       name="facLoadings" )
# Matrix holding the factor variance (Psy in greek/algebra)
FactorVariance <- mxMatrix( type="Symm", nrow=1, ncol=1, values=1, free=T,
                            labels="varF1",
                            name="facVariances" )
# Matrix holding the factor means (Kappa in greek/algebra)
FactorMean <-  mxMatrix( type="Full", nrow=1, ncol=1, values=0, free=F,
                         name="facMeans" )
# Matrix holding the residual variances (Theta in greek/algebra)
Residuals <- mxMatrix( type="Diag", nrow=4, ncol=4, free=T, values=1,
                       labels=c("e1","e2","e3","e4"),
                       name="resVariances" )
# Matrix holding the intercepts (Tau in greek/algebra)
Intercepts <- mxMatrix( type="Full", nrow=1, ncol=4, values=1, free=T,
                        labels=c("meanx1","meanx2","meanx3","meanx4"),
                        name="varMeans" )


ExpMean <- mxAlgebra(expression= varMeans + t(facLoadings %*% facMeans),
                     name="expMean" )
##############################################################################
### Algebra to calculate the expected COVARIANCE under the 1-factor model, ###
### ASSIGNMENT you have to write the expression (formula) yourself!        ###
##############################################################################

ExpCov <-mxAlgebra( expression= facLoadings %&% facVariances + resVariances,
                    name="expCov" )

Objective <- mxFIMLObjective( covariance="expCov", means="expMean", dimnames=names(dataset[,2:5]))

FactorModel <- mxModel("single factor model",data,loadings,FactorVariance,FactorMean,Residuals,
                       Intercepts,ExpMean,ExpCov,Objective)
 
FactorFit <- mxRun(FactorModel)
tableFitStatistics(fitSat,FactorFit)