---
title: "Bifactor ESEM Tutorial"
author: "Yuwei Li"
date: "2022-08-25"
output: 
  rmdformats::robobook:
  html_document: default
  word_document: default
editor_options:
  chunk_output_type: console
---

# Background

Within relational turbulence theory (RTT; Solomon et al., 2016), relational uncertainty is conceptualized as three facets of doubts and insecurities one experiences within a romantic association related to the self, the partner, and the relationship. Conventionally, RTT recognizes self and partner uncertainty as antecedents of relationship uncertainty, which in turn leads to a global characterization of the relationship as turbulent, and this causal model was supported by meta-analytic evidence (Goodboy et al., 2020). However, the correlations among these variables are consistently strong and positive, indicating that they may be influenced by a common factor. In Goodboy et al. (2021), they constrasted multiple measurement models for relational uncertainty with two samples from prior studies, and found the bifactor exploratory structural equation model to be the superior measurement model for the construct (data and Mplus code are available as supplemental materials of this article). I replicate these models in R using the lavaan package with their college sample data. 

Goodboy, A. K., Bolkan, S., Brisini, K., & Solomon, D. H. (2021). Relational uncertainty within relational turbulence theory: The bifactor exploratory structural equation model. Journal of Communication, 71(3), 403-430. https://doi.org/10.1093/joc/jqab009 

Goodboy, A. K., Bolkan, S., Sharabi, L. L., Myers, S. A., & Baker, J. P. (2020). The relational turbulence model: A meta-analytic review. Human Communication Research, 46(2-3), 222-249.

Solomon, D. H., Knobloch, L. K., Theiss, J. A., & McLaren, R. M. (2016). Relational turbulence theory: Explaining variation in subjective experiences and communication within romantic relationships. Human Communication Research, 42(4), 507-532.

# Prep work

Before you start, make sure you placed this .rmd file and "College Sample JOC.dat" in the same folder, so R can find and save any output files to the same place. 

**Reading the data into the environment.**

```{r}
# Reading the college sample data
college_sample<-read.csv("College-Sample-JOC.csv", header = FALSE)

# Because the data set does not contain variable names, 
# we will label variables here
colnames(college_sample)<-c(paste0("Self",1:6), # Self uncertainty
                            paste0("Part",1:6), # Partner uncertainty
                            paste0("Rel",1:6),  # Relationship uncertainty
                            paste0("Disc",1:4), # Disclosure to network
                            paste0("Lone",1:4), # Loneliness
                            paste0("Turb",1:4)) # Relational turbulence

# Coding missing data
college_sample[college_sample == -99]<-NA

# Subset the data (we only need the 18 relational uncertainty items)
college_ru<-college_sample[1:18]

# View the first 10 rows of college_ru
head(college_ru,10)
```

**Loading R packages**

RStudio may prompt you at the top of this panel to install packages that you have not used before. You can also install packages manually by typing in the console: for example, install.packages("dplyr"). 
```{r message=FALSE}
library(devtools)    # For version control
library(dplyr)       # For data management 
library(GPArotation) # For rotation used in factor analysis
library(lavaan)      # For SEM
library(lavaanPlot)  # For plotting SEM
library(psych)       # For data descriptives and reliability statistics
library(reshape2)    # For data manipulation
library(rlang)       # For string manipulation
library(semTools)    # For factor analysis 
```

# Data descriptives and reliabilities

We can first take a look at the item-level descriptives and inter-item correlations. 
```{r}
# Running descriptives
describe(college_ru)
# Bivariate correlations among items
round(cor(college_ru, use = "pairwise.complete.obs"),2)
```
We can examine the coefficient omega for each subscale. All three scales were reliable. 
```{r warning=FALSE,message=FALSE}
# Self uncertainty
omega(college_ru[1:6], nfactors = 1)$omega.tot

# Partner uncertainty
omega(college_ru[7:12], nfactors = 1)$omega.tot

# Relationship uncertainty
omega(college_ru[13:18], nfactors = 1)$omega.tot
```

# Fitting measurement models

Goodboy et al. (2021) fitted five such models, and we will replicate them here:  

* Unidimensional CFA  
* Independent Clusters Model-CFA (3-Factor CFA)  
* Exploratory Structural Equation Model  
* Bifactor CFA  
* Bifactor ESEM  

## Unidimensional CFA

This model considers relational uncertainty to be one unidimensional construct, and all 18 items serve as equivalent indicators for this construct. 

First, we write out the model in lavaan's language.  

+ `=~` represents "measured by"  
+ `~` represents "regressed on"  
+ `~~` represents "correlated with"  

Here, we specify a latent variable `relunc`, and that it is measured by all 18 variables in the scale.
```{r}
# Building the model
UNI.CFA <- 'relunc =~ Self1 + Self2 + Self3 + Self4 + Self5 + Self6
                    + Part1 + Part2 + Part3 + Part4 + Part5 + Part6 
                    + Rel1 + Rel2 + Rel3 + Rel4 + Rel5 + Rel6'
```

```{r}
# We run the model and save the results into fit.UNI.CFA
fit.UNI.CFA <- cfa(model = UNI.CFA,
                   data = college_ru,
                   # We use the maximum likelihood with robust SE estimator
                   estimator = "MLR",
                   # Full-information maximum likelihood for missingness
                   missing = "FIML") 

# We can also look at the graphic form of the model
lavaanPlot(name = "UnidimensionalCFA",
           model = fit.UNI.CFA)
```

Next, we extract the fit statistics from the object. Because it gives a lot of different measures, we can prespecify a subset of fit measures and use them to contrast all the models we will run. 

We can see the unidimensional CFA does not fit the data well.
```{r}
# Creating a subset of fit indices we will examine across all models
fit.subset<-c("chisq.scaled","df","pvalue.scaled",
              "rmsea.scaled","rmsea.pvalue.scale",
              "rmsea.ci.lower.scaled","rmsea.ci.upper.scaled",
              "cfi","tli","srmr","aic","bic")

# Extracting fit indices from the first model
fitmeasures(fit.UNI.CFA, fit.subset)
```

We will use the same workflow to fit the remaining models. 

## Independent Clusters Model-CFA (3-Factor CFA)  

The second model specifies self, partner, and relationship uncertainty as three correlated latent variables, and each is measured by six indicators with no cross-loading. It fits the data a little bit better compared to the unidimensional model. 
```{r}
# Building the model
FA3.CFA <- 'self =~ Self1 + Self2 + Self3 + Self4 + Self5 + Self6
            part =~ Part1 + Part2 + Part3 + Part4 + Part5 + Part6 
            rel  =~ Rel1 + Rel2 + Rel3 + Rel4 + Rel5 + Rel6'

# We run the model and save the results into fit.UNI.CFA
fit.FA3.CFA <- cfa(model = FA3.CFA,
                   data = college_ru,
                   # We use the maximum likelihood with robust SE estimator
                   estimator = "MLR",
                   # Full-information maximum likelihood for missingness
                   missing = "FIML") 

# We can also look at the graphic form of the model
lavaanPlot(name = "IndependentClustersModelCFA",
           model = fit.FA3.CFA)

# Extracting fit indices from the model
fitmeasures(fit.FA3.CFA, fit.subset)
```

## Exploratory Structural Equation Model 

The third model specifies self, partner, and relationship uncertainty as three correlated latent variables, and all items are allowed to cross-load onto all factors, though the target rotation will keep cross-loadings as small as possible. Fitting an ESEM in R is not as straightforward as in Mplus, and there is not a simple function for it. However, it is doable. 

Mateus Silvestrin's post explains the workflow for running ESEM in R: https://msilvestrin.me/post/esem/. Overall, there are three steps involved:

+ Running an exploratory factor analysis to obtain the loadings
+ Specifying the ESEM model using start values
+ Fitting the ESEM model

### Exploratory factor analysis with target rotation (oblique)
```{r}
# First, we conduct an EFA without any rotation
efa <- efaUnrotate(# Specifying the dataset
                   data = college_ru,
                   # 3 factors
                   nf = 3)

# Next, we construct a target matrix for the rotation
target.mat <- as.matrix(cbind(c(rep(NA,6),rep(0,12)),
                              c(rep(0,6),rep(NA,6),rep(0,6)),
                              c(rep(0,12),rep(NA,6))))

# View the target matrix
target.mat

# Rotate the matrix of EFA loadings 
rotated <- funRotate(efa, 
                     # Oblique target rotation
                     fun = "targetQ", 
                     # Specifying the target matrix
                     Target = target.mat)

# View the standardized loadings
rotated.loadings <- rotated@loading

# To use loadings as starting values for the ESEM, we need to 
# revert them back to the metric of the items. We do so by 
# multiplying the loadings by the items' standard deviations
rotated.loadings.unstd <- rotated.loadings * describe(college_ru)$sd
```

### Writing a function to generate model syntax

Because the model becomes really complex with all the items cross-loading onto each latent variable, we will write a few lines to generate our lavaan model. 
```{r}
# Manipulating the loadings matrix
loadings.esem <- melt(round(rotated.loadings.unstd,3))

# Naming the columns
colnames(loadings.esem) <- c("item", "latent", "value")

# Labeling our latent variables
loadings.esem$latent <- c(rep("self",18),
                          rep("part",18),
                          rep("rel",18))

# Generating the model syntax

# First, we set two anchors for each latent variables: one being the first item 
# in a corresponding subscale, and another not in the same subscale
anchors <- as.data.frame(cbind(latent = c("self","part","rel","self","part","rel"),
                                 item = c("Self1","Part1","Rel1","Part1","Rel1","Self1")))
# View the anchors dataframe
anchors

# I created a function for generating the model syntax "esem.lavaan.syntax"
esem.lavaan.syntax <- function (loadings.dt, anchors){
  
  # Create the is_anchor variable in the anchors dataframe
  anchors$is.anchor <- 1
  
  # Merge the anchors dataframe with the loadings dataframe, and rearrange it
  loadings.dt <- merge(anchors, loadings.dt, by = c("item","latent"), all.y = TRUE)
  loadings.dt$is.anchor[is.na(loadings.dt$is.anchor)]<-0
  loadings.dt <- arrange(loadings.dt, desc(is.anchor))  
  loadings.dt <- arrange(loadings.dt, latent)

  # Make syntax column per item; syntax is different depending on is.anchor
  loadings.dt <- loadings.dt %>%
    mutate(syntax = ifelse(is.anchor == 1, 
                           paste0(value,"*", item),
                           paste0("start(",value,")*", item)))
  
  # Combine starting/fix values and item names to generate a line
  # for each latent variable
  loadings.combined <- loadings.dt %>%
    group_by(latent) %>%
    summarize(syntaxline = paste(syntax, collapse = "+"))
  syntaxlines <- paste(loadings.combined$latent, "=~", loadings.combined$syntaxline)

  # Add lines to fix latent variable variance to 1
    c(syntaxlines,
    paste(loadings.combined$latent, "~~ 1*", loadings.combined$latent))
  
}

# Let's take a look at the generated syntax
esem.lavaan.syntax(loadings.esem, anchors)
```

The first three lines specify the measurement portion of the model. For example, `self =~ 1.375*Self1` fixes Self1's unstandardized loading on self uncertainty at that value, while `self =~ start(0.96)*Self2` gives the loading for Self2 a starting value without fixing it. Lines 4-6 sets the variance of each latent variable at 1. 

### Fitting the ESEM model 

Now that the new function has taken care of the syntax, we can just plug the function along with the input data frames we specified above into the CFA function. 
```{r}
fit.FA3.ESEM <- cfa(model = esem.lavaan.syntax(loadings.esem, anchors),
                    data = college_ru,
                    # We use the maximum likelihood with robust SE estimator
                    estimator = "MLR",
                    # Full-information maximum likelihood for missingness
                    missing = "FIML")

# We can also look at the graphic form of the model
lavaanPlot(name = "ESEM",
           model = fit.FA3.ESEM)

# Let's take a look at the standardized loadings 
select(arrange(arrange(parameterestimates(fit.FA3.ESEM, standardized = T)[1:54,], rhs),lhs),
       # only view a subset of output of interest
       lhs:est,pvalue,std.all)

# Extracting fit indices from the model
fitmeasures(fit.FA3.ESEM, fit.subset)
```

We can see from the standardized loadings that items mainly loaded onto the latent variables they belong to, while cross-loadings are permitted in this model. Note that the loadings for partner uncertainty are negative due to factor indeterminacy. We can reverse the signs on those loadings when we generate the model syntax. 

## Bifactor CFA

The fourth model specifies a bifactor configuration: a general relational uncertainty factor, and three residualized specific factors corresponding to self, partner, and relationship. This model requires factors be orthogonal to one another. 
```{r}
# Building the model
BIFA.CFA <- 'general =~ NA*Self1 + Self2 + Self3 + Self4 + Self5 + Self6
                      + Part1 + Part2 + Part3 + Part4 + Part5 + Part6 
                      + Rel1 + Rel2 + Rel3 + Rel4 + Rel5 + Rel6
                self =~ NA*Self1 + Self2 + Self3 + Self4 + Self5 + Self6
                part =~ NA*Part1 + Part2 + Part3 + Part4 + Part5 + Part6 
                rel  =~ NA*Rel1 + Rel2 + Rel3 + Rel4 + Rel5 + Rel6
                # Factors are orthogonal
                general ~~ 0*self + 0*part + 0*rel
                self ~~ 0*part + 0*rel
                part ~~ 0*rel
                # Setting latent variable variance at 1
                general ~~ 1*general
                self ~~ 1*self
                part ~~ 1*part
                rel ~~ 1*rel'

# We run the model and save the results into fit.UNI.CFA
fit.BIFA.CFA <- cfa(model = BIFA.CFA,
                   data = college_ru,
                   # We use the maximum likelihood with robust SE estimator
                   estimator = "MLR",
                   # Full-information maximum likelihood for missingness
                   missing = "FIML") 

# We can also look at the graphic form of the model
lavaanPlot(name = "BifactorCFA",
           model = fit.BIFA.CFA)

# Extracting fit indices from the model
fitmeasures(fit.BIFA.CFA, fit.subset)
```

## Bifactor ESEM

Here, we will follow the same steps as the previous ESEM. The difference between the two mostly stems from how latent factors are related in these models. In regular ESEM, factors correlate, so we use an oblique rotation for the EFA portion of the process, while bifactor ESEM sets factor correlations at 0, so orthogonal rotation is used in the EFA step.  

### Exploratory factor analysis with target rotation (orthogonal)
```{r}
# First, we conduct an EFA without any rotation
efa.bifac <- efaUnrotate(data = college_ru,
                         # 4 factors
                         nf = 4)

# Constructing a target matrix
target.mat.bifac <- as.matrix(cbind(rep(NA,18),
                              c(rep(NA,6),rep(0,12)),
                              c(rep(0,6),rep(NA,6),rep(0,6)),
                              c(rep(0,12),rep(NA,6))))

#View the target matrix
target.mat.bifac

# Using the target matrix to rotate the loadings
rotated.bifac <- funRotate(# The EFA results to be rotated
                           efa.bifac, 
                           # Target rotation (orthogonal)
                           fun = "targetT", 
                           # Specify the target matrix
                           Target = target.mat.bifac)

# To use loadings as starting values for the ESEM, we need to 
# revert them back to the metric of the items. We do so by 
# multiplying the loadings by the items' standard deviations
rotated.bifac.loadings.unstd <- rotated.bifac@loading*describe(college_ru)$sd


# Manipulate the loading matrix, name columns, label factors
loadings.esem.bifac <- melt(round(rotated.bifac.loadings.unstd,3))
colnames(loadings.esem.bifac) <- c("item","latent","value")
loadings.esem.bifac$latent<-c(rep("general",18),
                                 rep("self",18),
                                 rep("part",18),
                                 rep("rel",18))

# Setting anchors for the bifactor ESEM 
anchors.bifac <- as.data.frame(cbind(latent = c("general","general","general",
                                                "self","part","rel"),
                                     item = c("Self1","Part1","Rel1",
                                              "Self1","Part1","Rel1")))

#Examine the generated syntax
esem.lavaan.syntax(loadings.esem.bifac,anchors.bifac)

fit.BIFA.ESEM <- cfa(model = c(esem.lavaan.syntax(loadings.esem.bifac,anchors.bifac),
                               # Adding orthogonal constraints on factors
                               'general ~~ 0*self + 0*part + 0*rel
                                self ~~ 0*part + 0*rel
                                part ~~ 0*rel'),
                    data =  college_ru,
                    # We use the maximum likelihood with robust SE estimator
                    estimator = "MLR",
                    # Full-information meximum likelihood for missingness
                    missing = "FIML")

# We can take a look at the graphic form of the model
lavaanPlot(name = "BifactorESEM",
           model = fit.BIFA.ESEM)

# Examine the loadings 
select(arrange(arrange(parameterestimates(fit.BIFA.ESEM,standardized = T)[1:72,], rhs),lhs),
  # only view a subset of output of interest
  lhs:est,pvalue,std.all)

# Extracting fit indices from the model
fitmeasures(fit.BIFA.ESEM, fit.subset)
```

Notice that these loadings are within rounding errors with the parameters reported in Goodboy et al. (2021). We have successfully replicated the bifactor ESEM model. 

## Model Evaluation

Next, we can juxtapose the fit statistics of all five models. Again, some of these values slightly differ from Goodboy et al. (2021; Table 1, top panel), but this does not impact the substantive conclusions we can draw from this set of results. 

```{r}
fit.summary <- round(rbind(fitmeasures(fit.UNI.CFA, fit.subset),
                           fitmeasures(fit.FA3.CFA, fit.subset),
                           fitmeasures(fit.FA3.ESEM, fit.subset),
                           fitmeasures(fit.BIFA.CFA, fit.subset),
                           fitmeasures(fit.BIFA.ESEM, fit.subset)),3)

rownames(fit.summary)<-c("Uni CFA",
                         "ICM-CFA",
                         "ESEM",
                         "Bifactor CFA",
                         "Bifactor ESEM")

fit.summary
```

-----
### Additional Information

I created this tutorial with a system environment and versions of R and packages that might be different from yours. If R reports errors when you attempt to run this tutorial, running the code chunk below and comparing your output and the tutorial posted on the website may be helpful. 
  
```{r}
session_info(pkgs = c("attached"))
```
