#############################################################################
#################### Basics of R and introduction to lavaan #################
#################### PSYn5440 - Introduction to Factor Analysis #############

## "All that exists in R is an object. All that happens in R is a function call" 
##                                      -- someone whose name I don't remember

#################################################################
###### Getting data into R

# First, we need R to know what our data is. We'll be working with a correlation matrix, but the logic
# applies to raw data as well. I created a .csv file in Excel with the correlation matrix in it (each cell
# contains one element), this time, the matrix is full (i.e., not the triangular form) to avoid potential 
# issues. The first row was also used to name the columns of the matrix, with the names corresponding
# to the names of the MVs. 

# We read the .csv file using the "read.csv()" function:

read.csv(file = "C:/FSS/PSY544/datamatrix.csv", header = TRUE)

# The function requires us to specify some of its arguments, so that it knows what to do. 
# In this case, I have specified the full location of the .csv file (file = ...) [beware! the slashes are /, not \],
# and that the file contains a header (the top row with variable names) by assigning a "TRUE" [can be abbreviated as "T"]
# value to the "header" argument.

# However, when we call the function, we just get the file content printed into the console, although what we want
# is to store the data in R. We need to assign the function output into an object that we name in some way, like "cormatrix":

cormatrix <- read.csv(file = "C:/FSS/PSY544/datamatrix.csv", header = T)

# There are different data structures [ways how data are represented] in R. The "read.csv" implicitly creates the
# so-called "data frame" structure. We want a matrix, so we tell R to convert the object "cormatrix" into a matrix and
# store the result into the same object (effectively, we are replacing the former object):

cormatrix <- as.matrix(cormatrix)

# You can again call the object to see if the content is right by printing it into the console:

cormatrix

#################################################################
###### Running lavaan

# All-right! Now that we have the data, we need to load the lavaan package into the R environment so that we can use it. 
# To load a package, you have to install it first if it was not previously installed on your machine:

install.packages("lavaan")

# Once the package is installed, you can simply load it into the environment by calling the library() function:

library(lavaan)

######## Okay, now, onto the lavaan syntax. 
# Since we're doing CFA, we require to have a hypothesis ready beforehand, which we will test by fitting a model to the data. 

# Our first hypothesis is that there is only a single factor affecting all nine MVs. We shall call the factor "F1" for
# the sake of the syntax, but remember that this is just a label. 

# We need to represent our hypothesis by setting up the model matrices Lambda, Phi (which lavaan calls Psi) and 
# D-psi (which lavaan calls Theta)
# So let's get it straight - we have 9 MVs and we hypothesize about a single factor, so our Lambda will be a 9 x 1 matrix.
# We only hypothesize about a single factor, so our Phi matrix will be a 1x1 matrix.
# The D-psi matrix will be a diagonal matrix, like we're used to. 

# Now, it's important to understand how the lavaan syntax works in this case. 
# In this simple CFA scenario, the lavaan defaults are as follows: All the Lambda elements are fixed to zero, all the
# Phi elements are free, all the diagonal elements of D-psi are free and its off-diagonal elements are fixed to zero. 

# We will formulate our model by changing some of these defaults so that the matrices conform to what we want.

### Lambda matrix: ##
# As said, the Lambda matrix is full of zeroes by default. Thus, we need to free up individual elements of this matrix
# in our syntax. Lavaan doesn't know how many columns we want in the Lambda matrix, so we will set the number of columns 
# by naming (and thus creating) the latent variables - if a latent variable isn't named in the syntax, it does not 
# exist in the eyes of lavaan and the respective column in the Lambda matrix is not created. 
# It's not enough to just *name* the factor, we right away need to tell lavaan what MVs we hypothesize load onto that factor,
# so that lavaan can free up all the elements in the column corresponding to the named factor and in the rows corresponding to 
# the named MVs. 

# This is done by naming the factor on the left-hand side, followed by '=~', followed by the names of the MVs separated
# by '+' signs, like so:

### F1 =~ WrdMean + SntComp + OddWrds
# -- this creates a column in the Lambda matrix labeled F1, and frees up elements in rows of the Lambda matrix that
# correspond to the MVs that are labeled "WrdMean", "SntComp" and "OddWrds". All other elements will be left fixed to 0 by default. 

### Phi: ##
# As said, the Phi matrix is full of free parameters by default. If we want to fix some of the parameters, we need to 
# do so in our syntax. Because we began with setting up the Lambda matrix, Lavaan now knows how many rows/columns does 
# the Phi matrix have. All that we pretty much have to do is tinker with the off-diagonal elements (factor correlations)
# in case we want to fix them. 

# We fix or free up elements in the Phi matrix by first naming a given factor on the left side, followed by '~~', followed
# by the name of the other factor involved in what we are about to fix / free up. 

### F1 ~~ 0*F2
# -- this fixes the element in the F1 row and F2 column to 0.

## D-psi:
# This matrix is, by default, diagonal with the diagonal elements free and all the other elements fixed to 0. 
# That's what we want, so we don't have to specify the D-psi matrix any further in the syntax, since we are just using
# the default setting of the matrix.

##### Setting up the syntax for the one-factor model #####

# Now, we need to write the syntax and store it somewhere, inside an object, so that lavaan/R can work with it. 
# We write the syntax as a string of characters, by enclosing the whole of it in quotation marks, ' ' or " ". 
# Below, we are going to write the model syntax and store it into an object called model1F (arbitrary name). 
# The pound signs "#" ("hashtags") denote comments, and the syntax parser ignores whatever follows these symbols. 

model1F <- '   # <- this opening quotation mark says "the syntax begins here"

### Set elements in Lambda matrix.  

# Our hypothesis is that all nine MVs load onto a single factor that we are going to label F1:
# (the first loading on each factor is set to 1 by default. Override with NA*, like shown below)

F1 =~ WrdMean + SntComp + OddWrds + MxdArit + Remndrs + MissNum + Gloves + Boots + Hatchts

### No need to set anything in the Phi matrix, we are using the default
### No need to set anything in the D-psi matrix, we are using the default

' # <- this closing quotation mark says "the syntax ends here"

# That's it, the model matrices are set up. The syntax containing the instructions is stored inside an object
# called model1F

##### Fitting the model #####

# Syntax isn't enough - we also need to choose our data, our discrepancy function (and possibly other things),
# and we need to call the model fitting function cfa() to fit the model, which we will store inside an object
# called fit1F. 

# the cfa() function has many potential arguments, separated by commas. Every argument has a name. 
# Inside the function call, we name the respective argument and assign something to it - either an R object,
# or a number, or a string. 

# the arguments we are going to use are the following:
# "model" -- what is the object where the model syntax is stored. A name of an existing object.
# "sample.cov" -- what is the object where the correlation/covariance matrix is stored. A name of an existing object.
# "sample.nobs" -- how many persons did we have in our sample? A number.
# "estimator" -- what discrepancy function do we want to use? A string [enclosed by quotation marks] with the name of the estimator. 
#                we are going to use the Unweighted Least Squares (also known as OLS, Ordinary Least Squares) discrepancy function,
#                which is abbreviated ULS in lavaan
# "std.lv" -- do we want the variances of the common factors to be fixed to 1? We set this to TRUE because that's the approach we have been following.

# We call the function with the respective arguments, and store the result of fitting the function into an object
# that we will call fit1F:

fit1F <- cfa(model = model1F, sample.cov = cormatrix, sample.nobs = 696, estimator = "ULS", std.lv = TRUE)

### After a little while, the model is estimated - now, we have an object (fit1F) with the fitted model at our disposal.

# First, let's check if the model matrices seem right. The lavInspect() function takes the object as an argument,
# and gives out some stuff - like an overview of the model matrices. We can look at the entire thing:

lavInspect(fit1F)

# Or individual matrices:

lavInspect(fit1F)$lambda
lavInspect(fit1F)$psi
lavInspect(fit1F)$theta

# Doing this, we can see what elements are fixed to zero and what elements were free (the numbered ones)

# We can also ask for the actual estimates within the matrices:

lavInspect(fit1F, what = "est")

lavInspect(fit1F, what = "est")$lambda
lavInspect(fit1F, what = "est")$psi
lavInspect(fit1F, what = "est")$theta

# The summary() function will give us more information about the fitted model:

summary(fit1F)

# The residuals() function outputs the residual matrix:

residuals(fit1F)
round(residuals(fit1F)$cov,2)

library(corrplot)
corrplot(cormatrix)


# ...we can also ask for some fit statistics, fit indices and other fit-related information using the fitmeasures() function
# (there's a lot of things the fitmeasures() function can give us, here I'm just showing a couple)
# If you drop the 'fit.measures' argument, you will get the full list of various measures of fit and related information

fitmeasures(fit1F, fit.measures = c("rmsea", 
                                     "tli", "df","npar",
                                     "chisq",
                                     "rmsea.ci.lower","rmsea.ci.upper", "cfi", "srmr"))

### We can also try to estimate the model parameters using Maximum Likelihood. 
### Simply change the value of the "estimator" argument to "ML", like so:

fit1F <- cfa(model = model1F, sample.cov = cormatrix, sample.nobs = 696, estimator = "ML", std.lv = TRUE)

### As can be seen, the rest of the arguments is the same. If you would want to compare the outputs more easily,
### I would suggest putting the resulting model into a different object than where the ULS-estimated model lives,
### so that you don't replace the contents of the already-created object. Let's name this one "fit1F_ML":

fit1F_ML <- cfa(model = model1F, sample.cov = cormatrix, sample.nobs = 696, estimator = "ML", std.lv = TRUE)

######## Okay, that was the basics - a one-factor model estimated using ULS or ML. 
######## Let's take a look at two more models with the same data. 

######## Three-factor model with uncorrelated factors ####

### This model represents the hypothesis that there are three common factors out there:
### A verbal ability factor (associated MVs are WrdMean, SntComp and OddWrds)
### A numerical ability factor (associated MVs are MxdArit, Remndrs and MissNum)
### A spatial ability factor (associated MVs are Gloves, Boots and Hatchts)

### In addition, in this case we would hypothesize the factors to be uncorrelated. 

# Remember, the model matrices are as follows: Lambda, Phi (which lavaan calls Psi) and D-psi (which lavaan calls Theta)
# We have 9 MVs and we hypothesize about three factors, so our Lambda will be a 9 x 3 matrix.
# Our Phi matrix will be a 3x3 matrix, with zeroes off the diagonal
# The D-psi matrix will be a diagonal matrix, like we're used to. 

# Let's store the model syntax in an object, and let's call it "model3F":

model3F <- '

###Set elements in Lambda matrix

Verbal =~ WrdMean + SntComp + OddWrds
Numerical =~ MxdArit + Remndrs + MissNum
Spatial =~ Gloves + Boots + Hatchts

###Set the off-diagonal elements in Phi matrix to zero

Verbal ~~ 0*Numerical
Verbal ~~ 0*Spatial
Numerical ~~ 0*Spatial
'

# Fit the model using the cfa() function and store the result in an object called "fit3F". Here we are using ULS:

fit3F <- cfa(model3F, sample.cov = cormatrix, sample.nobs = 696, estimator = "ML", std.lv = TRUE)

# Actually, you could have omitted setting the off-diagonal elements in Phi matrix to zero in the model syntax. 
# Instead, you can fit the model using the cfa() function and utilize the "orthogonal" argument. The argument is set to "FALSE"
# by default [meaning the off-diagonal elements in Phi are by default free], but you can set it to "TRUE" instead, which 
# overrides the default and fixes the off-diagonal elements to zero:

fit3F <- cfa(model3F, sample.cov = cormatrix, sample.nobs = 696, estimator = "ULS", std.lv = TRUE, orthogonal = TRUE)

# Both work the same - feel free to compare:

summary(fit3F)

lavInspect(fit3F)$lambda
lavInspect(fit3F)$psi
lavInspect(fit3F)$theta

lavInspect(fit3F, what = "est")$lambda
lavInspect(fit3F, what = "est")$psi
lavInspect(fit3F, what = "est")$theta

residuals(fit3F)

fitmeasures(fit3F, fit.measures = c("rmsea", 
                                    "tli", "df","npar",
                                    "chisq",
                                    "rmsea.ci.lower","rmsea.ci.upper", "cfi", "srmr"))

######## Three-factor model with correlated factors ####

### This model represents the hypothesis that there are three common factors out there:
### A verbal ability factor (associated MVs are WrdMean, SntComp and OddWrds)
### A numerical ability factor (associated MVs are MxdArit, Remndrs and MissNum)
### A spatial ability factor (associated MVs are Gloves, Boots and Hatchts)

### In addition, in this case we would hypothesize the factors to be CORRELATED. 

# This model is identical to the previous model, except for not being orthogonal
# As such, the Phi matrix will be a 3x3 matrix, with free elements off the diagonal

# Let's store the model syntax in an object, and let's call it "model3Fc":

model3Fc <- '

###Set elements in Lambda matrix

Verbal =~ WrdMean + SntComp + OddWrds
Numerical =~ MxdArit + Remndrs + MissNum
Spatial =~ Gloves + Boots + Hatchts

###Because the off-diagonal elements in Phi are by default free, we do not have to do much else.

'
fit3Fc <- cfa(model3Fc, sample.cov = cormatrix, sample.nobs = 696, estimator = "ML", std.lv = TRUE)

summary(fit3Fc)

lavInspect(fit3Fc)$lambda
lavInspect(fit3Fc)$psi
lavInspect(fit3Fc)$theta

lavInspect(fit3Fc, what = "est")$lambda
lavInspect(fit3Fc, what = "est")$psi
lavInspect(fit3Fc, what = "est")$theta

residuals(fit3Fc)

fitmeasures(fit3Fc, fit.measures = c("rmsea", 
                                    "tli", "df","npar",
                                    "chisq",
                                    "rmsea.ci.lower","rmsea.ci.upper", "cfi", "srmr"))

#### Comparing models directly - a Chi square difference test

## You may have realized by now that the model stored in fit3F is nested within the model stored in fit3Fc - 
## in other words, the uncorrelated model is a special version of the more general correlated model. 

## Since the models are nested, we can conduct a direct comparison using the chi-square difference test.
## For these purposes, we will utilize the anova() function, where the two model objects will act as arguments:

anova(fit3F, fit3Fc)

## According to the difference test, the correlated model yields better model-data fit. 