---
title: "Meta-Analytic SEM Tutorial"
author: "Yuwei Li"
date: "2022-08-28"
output: 
  rmdformats::robobook:
  html_document: default
  word_document: default
editor_options:
  chunk_output_type: console
---
# Background

Meta-analytic structural equation is a technique that compiles cumulative empirical evidence from meta-analysis to form a matrix that can in turn be analyzed using SEM. I had developed a workflow for meta-analytic SEM for a forthcoming project, and adapted it to replicate Rains' analysis on the nature of psychological reactance to demonstrate this technique. Another, more recent example of this technique is Goodboy et al.'s (2020) meta-analysis on relational turbulence model. 

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.

Rains, S. A. (2013). The nature of psychological reactance revisited: A meta-analytic review. Human Communication Research, 39(1), 47-73.

# Prep work

Before you start, make sure you have the data file "meta-sem-data.xlsx" stored in the same directory as this .rmd file, so R can find it and save any output files to the same place. 

**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(readxl)      # For reading in excel spreadsheet
library(lavaan)      # For SEM
library(lavaanPlot)  # For plotting SEM
library(psych)       # For computing harmonic mean
library(psychmeta)   # For meta-analysis
library(reshape2)    # For data manipulation
library(reshape)     # For data manipulation
```

**Reading the data into the environment.**  

You'll notice that effect sizes dataset has 4 columns: 

+ `study`: study identifier
+ `var1` & `var2`: variable pair corresponding with the correlation effect size
+ `r`: correlation effect size

The study sample size dataset has the same study identifier column, in addition to sample size `n`. 
```{r}
# The two worksheets contain effect sizes and study sample size information respectively

ES <- read_xlsx("meta-sem-data.xlsx",1)
StudyN <- read_xlsx("meta-sem-data.xlsx",2)

# View the first 10 rows of these data frames
head(ES,10)
head(StudyN,10)

# Let's add the sample sizes to the effect sizes
ES <- merge(ES, StudyN, by = "study")
head(ES,10)
```

Before we jump into the meta-analysis, let's take a look at the numbers of studies and effect sizes available for each correlation. 
```{r}
ES %>%
  group_by(var1,var2) %>%
  summarize(no.studies = n(), total.n = sum(n))
```

There are 6 cells in the matrix of 4 variables (lower half, and excluding the diagonal), and we can see for each meta-analysis, we have 9-19 effect sizes, and the total sample size for each cell was between 2334 and 4821. 

# Create a matrix with meta-analysis results

First, let's review how to the psychmeta package to conduct meta-analysis. We will use the correlation between threat and anger as an example. 
```{r}
# Subset the data for the correlation of interest
ES.subset <- filter(ES, var1 == "threat" & var2 == "anger")

cor.threat.anger <- ma_r(# specifying the column of effect sizes
                         rxyi = r,
                         # the column for numbers
                         n = n,
                         # we will use the barebones method
                         ma_method = "bb",
                         # giving names to focal variables
                         construct_x = ES.subset$var1[1],
                         construct_y = ES.subset$var2[1],
                         # the input dataframe 
                         data = ES.subset)

# Let's view the results of the first meta 
first.meta<-as.data.frame(cor.threat.anger$meta_tables$`analysis_id: 1`$barebones)
first.meta
```

The barebones (because we didn't adjust for any artifacts) results indicate that on average, threat and anger are correlated at .217. 

Next, we write a loop function to run a meta-analysis for each cell in the matrix.

```{r}
# Let's add an identifier for each effect
ES <- ES %>%
  mutate(ESid = paste(var1, "~~", var2))

# Create an empty data frame to save the results into
results <- NULL

# Create a loop to run the same procedure for each cell
# Specify the elements used for the loop
for (cell in unique(ES$ESid)) {
  # Subset the effect size table
  ES.subset <- ES %>%
    filter(ESid == cell)
  # Run the meta-analysis
  meta.subset <- ma_r(# specify the column of effect sizes
                      rxyi = r,
                      # the column for numbers
                      n = n,
                      # we will use the barebones method
                      ma_method = "bb",
                      # give names to focal variables
                      construct_x = ES.subset$var1[1],
                      construct_y = ES.subset$var2[1],
                      # the input dataframe 
                      data = ES.subset)
  # Using the headers of the first meta-analysis we ran as a base,
  # save the results of each new meta-analysis to it
  subset.results <- as.data.frame(c(var1 = as.character(meta.subset$construct_x),
                  var2 = as.character(meta.subset$construct_y),
                  meta.subset$meta_tables$`analysis_id: 1`$barebones))
  
  results <- rbind(results,subset.results)
}

# View the results
results
```

Alright! We have all the effect sizes we need to build a matrix for SEM/path analysis. 
```{r}
# Select the mean effect size estimates
mean.cors <- select(results, var1, var2, mean_r)

# Manipulate the data into matrix form
cor.matrix <- cast(mean.cors, var2~var1)

# Complete the matrix
cor.matrix$attitude <- c(NA, NA, NA)
cor.matrix <- rbind(c(var2 = "threat", rep(NA,4)),
                    cor.matrix)

# Rearrange the matrix
cor.matrix <- cor.matrix %>% 
  arrange(factor(var2, levels = c("threat","anger","counterargue","attitude"))) %>%
  select(var2, threat, anger, counterargue, attitude)

# Fill in 1's on the diagonal 
cor.matrix[1,2]<-1
cor.matrix[2,3]<-1
cor.matrix[3,4]<-1
cor.matrix[4,5]<-1

# Fill in the matrix
cor.matrix <- cor.matrix %>%
  mutate(cor.input = paste(threat,anger,counterargue,attitude))

input.matrix <- gsub("NA","",as.character(cor.matrix$cor.input))

input.matrix <- getCov(input.matrix, names = colnames(cor.matrix[2:5]))

# View the matrix 
input.matrix
```

# Testing Models

Next, we will build alternative SEM models to determine the best configuration of reactance. 

Rains (2013) tested 3 such models. 

+ A dual-process model, where anger and counterarguing functioned in a parallel manner to convey the influence of threat to freedom to attitude.
+ An intertwined model, where anger and counterarguing are two components of a latent construct reactance 
+ A linear affective-cognitive model, where anger precededs counterarguing

Let's build these models: 

As a reminder, in lavaan language -- 

+ `=~` represents "measured by"  
+ `~` represents "regressed on"  
+ `~~` represents "correlated with"  
```{r}
dual.process.model <- "anger ~ threat
                       counterargue ~ threat
                       attitude ~ anger + counterargue"

intertwined.model <- "reactance =~ anger + counterargue
                      reactance ~ threat
                      attitude ~ reactance"

linear.model <- "anger ~ threat
                 counterargue ~ anger
                 attitude ~ counterargue"
```

In meta-analytic SEM, some scholars opt to use a sample size that represents the collection of studies, such as a mean, median, or minimum *n* of selected studies (such as Rains), while others compute an average of sorts of all total *n's* for each cell (Landis, 2013). Here, we take the latter approach by calculating a harmonic mean for all cells. 

A harmonic mean is the reciprocal of the arithmetic mean of the reciprocals of the given set of observations. For example, the harmonic mean of 1, 10, and 100 is 2.70. 

$$\begin{aligned}
H &= (\frac{\frac{1}{1} + \frac{1}{10} + \frac{1}{100}}{3})^{-1}  \\
&= \frac{3}{1+.1+.01} = \frac{3}{1.11}\\
&= 2.70
\end{aligned}$$

Landis, R. S. (2013). Successfully combining meta-analysis and structural equation modeling: Recommendations and strategies. Journal of Business and Psychology, 28, 251–261. doi: 10.1007/s 10869-013-9285-x.
``` {r}
harmonic.mean(c(1,10,100))

# Fit all three models
fit.dual.process <- sem(model = dual.process.model,
                        # Using the cor matrix as input
                        sample.cov = input.matrix,
                        # Computing the harmonic mean of n's corresponding 
                        # to each cell in the matrix
                        sample.nobs = harmonic.mean(results$N))

fit.intertwined <- sem(model = intertwined.model,
                       sample.cov = input.matrix,
                       sample.nobs = harmonic.mean(results$N))

fit.linear <- sem(model = linear.model,
                  sample.cov = input.matrix,
                  sample.nobs = harmonic.mean(results$N))

# Creating a subset of fit indices we will examine across all models
fit.subset<-c("chisq","df","pvalue",
              "rmsea","rmsea.pvalue",
              "rmsea.ci.lower","rmsea.ci.upper",
              "cfi","tli","srmr","aic","bic")

# Extracting fit indices from all three models
round(rbind(fitmeasures(fit.dual.process, fit.subset),
      fitmeasures(fit.intertwined, fit.subset),
      fitmeasures(fit.linear, fit.subset)),3)


# Let's take a look at the model in graphic form and the parameters
lavaanPlot(model = fit.intertwined,
           coefs = T,
           stand = T,
           sig = .05,
           digits = 3)

parameterestimates(fit.intertwined, standardized = T)
```

The fit statistics tell us Dillard & Shen's (2005) intertwined model was indeed the best fit for the data. 
  
You can also add to the analysis by providing reliability information to correct the measurement errors present in barebones meta-analysis. Doing so will result in a matrix of true correlations, and this will remove any attenuation in effects due to imperfect measurement of constructs.  

-----
### 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"))
```
