Plot interaction effect in sem model with observed variables in R

Question

I am estimating an SEM model that has observed variables. I am using SEM to handle missing data using FIML. My model has an interaction term to test for moderation. Here is a toy example that illustrates the issue.

library(lavaan)
library(car)
library(dplyr)

data(starwars)

sw2 <- starwars %>% mutate(
  male = Recode(sex, "'male' = 1; NA=NA; else = 0"),
  human = Recode(species, "'Human' = 1; NA=NA; else = 0"),
  maleXby = male * birth_year,
)

mod <- 'mass ~ height + human + male + birth_year + maleXby'
fit <- sem(mod, data = sw2, missing="fiml.x")
summary(fit)

What I want to do is plot the interaction term like a margin plot, to visualize the interaction effect. But package like library(interactions) does not work with an object of class lavaan. How could I visualize this? Is there a package (like interactions) that makes this easier.

Answer 1

You could fit this model using lm(), but I think you want to be able to use FIML estimates, yes? In that case, you could use the emmeans package, which can work on lavaan-class objects if you have the semTools package loaded.

You didn't say which predictor was focal vs. moderator, but I assume you want to treat male as moderator because it is a grouping variable. The example below can be adapted by switching their roles in the pairs() function, as well as by selecting different birth_year levels at= which to probe the effect of male. When birth_year is the focal predictor, its linear effect will be the same regardless of which levels are chosen, so I chose the full range() below.

library(emmeans)
library(semTools)

## for ease of use, fit model using colon operator
mod <- 'mass ~ height + human + male + birth_year + male:birth_year'
fit <- sem(mod, data = sw2, missing = "fiml.x")

## calculate expected marginal means for multiple 
## levels of male (1:0) and birth_year
BYrange <- range(sw2$birth_year, na.rm = TRUE)
em.mass <- emmeans(fit, specs = ~ birth_year | male, 
                   at = list(male = 1:0, birth_year = BYrange),
                   # because SEMs can have multiple DVs:
                   lavaan.DV = "mass")
em.mass
## probe effect of year across sex
rbind(pairs(em.mass))
## plot effect of year across sex
emmip(em.mass, male ~ birth_year)   # 2 lines in same plot
emmip(em.mass, ~ birth_year | male) # in separate panels