How can one go about modelling the interaction between a categorical independent variable and a continous moderator (created through a CFA) in a SEM model using the lavaan package in R?

In particular, in my real dataset I am essentially interested in re-creating a two-way ANOVA in SEM, and also want to include a moderating variable to test with each factor variable.

Example data and problem:

### load packages: ###

### Create some data: ###

# Dependent variable: Taken example data from psychTools
DV1 <- bfi$A1 # item 1 
DV2 <- bfi$A2 # item 2
DV3 <- bfi$A3 # item 3

# Moderating variable: Taken example data from psychTools
MOD1 <- bfi$C1 # item 1
MOD2 <- bfi$C2 # item 2
MOD3 <- bfi$C3 # item 3

# Create example factor variables
x1 <- c("A","B")
x2 <- c("C","D")
FAC1 <- as.factor(sample(x1, 200, replace = TRUE)) # Factor 1, with two levels "A" and "B"
FAC2 <- as.factor(sample(x2, 200, replace = TRUE)) # Factor 2, with two levels "C" and "D"
FAC12 <- interaction(FAC1,FAC2) # Factor 12, interaction of FAC1 and FAC2, with four levels "A.C" "B.C" "A.D" "B.D"

# Combine to data frame
StudyData <- data.frame(DV1,DV2,DV3, 
                        FAC1, FAC2, FAC12)

### Make all categorical variables numeric for use in SEM (orthogonal contrast coded as in ANOVA): ###

StudyData$FAC1 <- recode(StudyData$FAC1, "c('A')='-1';
StudyData$FAC1 <- as.numeric(levels(StudyData$FAC1))[StudyData$FAC1]

StudyData$FAC2 <- recode(StudyData$FAC2, "c('C')='-1';
StudyData$FAC2 <- as.numeric(levels(StudyData$FAC2))[StudyData$FAC2]

StudyData$FAC12 <- recode(StudyData$FAC12, "c('A.D','B.C')='-1';
StudyData$FAC12 <- as.numeric(levels(StudyData$FAC12))[StudyData$FAC12]

### SEM Model One: ###
Model.one <- '
# cfa
DV =~ DV1 + DV2 + DV3
MOD =~ MOD1 + MOD2 + MOD3
# regressions
DV ~ FAC1 + FAC2 + FAC12

Modelone <- sem(Model.one, StudyData, estimator="MLM", effect.coding=TRUE, meanstructure=TRUE) 
fitMeasures(Modelone, c("chisq","cfi","rmsea","srmr","nfi","gfi"))

### SEM Model Two: ###
Model.two <- '
# cfa
DV =~ DV1 + DV2 + DV3
MOD =~ MOD1 + MOD2 + MOD3
# regressions
DV ~ FAC1 + FAC1:MOD + FAC2 + FAC12

Modeltwo <- sem(Model.two, StudyData, estimator="MLM", effect.coding=TRUE, meanstructure=TRUE) 
fitMeasures(Modeltwo, c("chisq","cfi","rmsea","srmr","nfi","gfi"))

### EDIT ###

### SEM Model Three: ###
Model.three <- '
# cfa
DV =~ DV1 + DV2 + DV3
MOD =~ MOD1 + MOD2 + MOD3
# regressions

Modelthree <- sem(Model.three, StudyData, estimator="MLM", effect.coding=TRUE, meanstructure=TRUE, group="FAC1") 
fitMeasures(Modelthree, c("chisq","cfi","rmsea","srmr","nfi","gfi"))

Model one runs fine. I can run my "ANOVA" in the SEM environment.

However, when I want to run Model two, which includes an interaction term between FAC1 and MOD (as created via CFA in the SEM model), I receive the error:

"lavaan WARNING: The variance-covariance matrix of the estimated parameters (vcov) does not appear to be positive definite! The smallest eigenvalue (= -3.458498e-20) is smaller than zero. This may be a symptom that the model is not identified."


  1. Is it not possible to create a factor:continuous interaction in lavaan in this manner?
  2. Are there any work arounds & how to do them? (For example, extract the values calculated during the CFA for MOD and calculate FAC1:MOD interaction outside of the SEM, then re-use the variable in the path analysis (regressions) part of the SEM)
  3. Can Mplus do this without the need for work arounds?

1 Answer 1

  1. No, the : operator only works on observed variables. It triggers lavaan to actually calculate the product term and include it in the covariance matrix and mean vector to which the model is fitted. That cannot be what happens with a latent variable, which is not part of the observed variables' summary statistics.

  2. Moderation is symmetric, so you could use a multigroup model, with the categorical IV as the grouping variable. Differences in the DV~MOD simple effects across groups would be moderation by FAC1. Differences in the DV~1 intercepts across groups would be the simple effect of FAC1, which could be probed by centering MOD's mean at different values. Or you might be able to use the emmeans utility in the semTools package to probe the interaction; see the ?lavaan2emmeans help-page examples. I suggest another possibility below.

    • Note that using this approach would make the Fac12 variable redundant because Fac2's intercept and effects are likewise moderated by Fac1 by virtue of parameters differing across groups.
  3. Yes, Mplus can simply use LMS estimation, but that is fraught with some restrictive assumptions. My student's PhD research (also here) has revealed that the product-indicator approach is less restrictive, and that can be implemented in lavaan (see her tutorial about using this method for invariance testing, available for download from my faculty page).

  • $\begingroup$ On 2: I understand I can add a single categorical grouping variable e.g., 'group ="FAC1"', referring to my code this could be FAC1, FAC2, or the interaction term FAC12. However, in factorial ANOVA it is common to include both main effects and the interaction term. Thus, the grouping variable method is not suitable for this, correct? I won't be able to estimate main effects and the interaction in an ANOVA-like fashion (as code for 'SEM Model One' currently does). Isn't it misleading to not include main effects and only examine the interaction term levels with a grouping var? $\endgroup$
    – JoeyyyFunk
    Commented Nov 27, 2021 at 8:26
  • 1
    $\begingroup$ In your OP, mod only interacted with Fac1, so interpret my (2) in that context. But using group="Fac1" still includes its interaction with Fac2 because Fac2's effects can differ across Fac1 groups. But that means you can (must) exclude the redundant Fac12 from your model. I will edit my post above to state that. $\endgroup$
    – Terrence
    Commented Nov 29, 2021 at 11:52
  • $\begingroup$ But by doing it the way you mention, I would be unable to estimate main effects for FAC2 or MOD right? As they always depend on the grouping variable. To refer to my new example (see SEM model 3, original post edited), this basically shows the interaction between FAC1 and FAC2, and FAC1 and MOD. If I wanted to estimate main effects of FAC1, I could instead use DV ~ 0*FAC2 + 0*MOD… but there is no way to estimate FAC2 or MODs main effects while still using a grouping variable, correct? $\endgroup$
    – JoeyyyFunk
    Commented Nov 29, 2021 at 19:45
  • 1
    $\begingroup$ If you want main effects, you can either define new parameters (with the := operator) in your Model 3 script that are the average intercepts across groups, or you can fit a model with their effects constrained to equality across groups. With balanced/equal Ns, both should provide the same result. With unequal Ns, the latter provides estimates that are weighted (by N) averages, as opposed the unweighted average effects in the former. Either way, the MOD effects are only comparable across groups if you impose measurement-invariance constraints on MOD's loadings and intercepts. $\endgroup$
    – Terrence
    Commented Nov 30, 2021 at 11:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.