Mediated moderation analysis with two mediators in R

This is my df. There is a significant interaction between IV1 and IV2. I'd like to perform a multiple mediation analysis where the two mediators (M1 and M2) explain the association between IV1 and DV. IV2 moderates the association between the DV and the mediators.

id <- rep(c(1, 2, 3), each = 4)
IV1 <- rep(c(2, 1), each = 2, times = 3)
IV2 <- rep(c(1, 2), times = 6)
DV <- c(4,3,2,4,4,4,5,4,4,4,3,5)
M1 <- c(1,1,1,3,3,5,3,4,3,4,2,4)
M2 <- c(3,4,4,4,4,5,4,5,4,4,4,4)
df2 <- data.frame(id, IV1, IV2, DV, M1, M2)
df2
> df2
# A tibble: 696 × 6
      ID IV1   IV2      DV    M1    M2
   <dbl> <fct> <fct> <dbl> <dbl> <dbl>
 1     1 2     1         4     1     3
 2     1 2     2         3     1     4
 3     1 1     1         2     1     4
 4     1 1     2         4     3     4
 5     2 2     1         4     3     4
 6     2 2     2         4     5     5
 7     2 1     1         5     3     4
 8     2 1     2         4     4     5
 9     3 2     1         4     3     4
10     3 2     2         4     4     4
# … with 686 more rows
# ℹ Use `print(n = ...)` to see more rows

I tried mediation and laavan packages but failed to get results.

model <- "
  # Regression paths from IV to mediators
  M1 ~ IV1 * IV2
  M2 ~ IV1 * IV2
  
  # Regression paths from mediators to DV
  DV ~ IV1 + IV2
  
  # Indirect effects
  indirect1 := IV1 * IV2 * M1
  indirect2 := IV1 * IV2 * M2
"

fit <- lavaan::cfa(model, data = df2)

> Error in lav_partable_constraints_def(partable, con = LIST, debug = debug) : 
  lavaan ERROR: unknown label(s) in variable definition(s): IV2 M1 M2

# Extract mediation results
mediation_results <- lavaan::parameterEstimates(fit)

I think the problem is how you're specifying the interactions. Rather than using IV1*IV2, you need IV1 + IV2 + IV1:IV2. However, it may be easier to calculate the indirect effects if you supply parameter names for the paths. You can do this by putting the parameter name before the variable and then separating the parameter and variable name with an asterisk (e.g., a1*IV1 + a2*IV2 + a3*IV1:IV2). The second problem is how you're calculating the indirect effect. If you think about the mathematical notation of the model, you're proposing something like this (as far as I can tell)

In the code below, I use three different levels of IV2 to calculate the indirect effect of IV1: the mean (1.5), the mean minus a SD (1) and the mean plus a SD (2). I did the same for the values of IV1 in the indirect effect of IV2.

library(lavaan)
#> This is lavaan 0.6-15
#> lavaan is FREE software! Please report any bugs.
id <- rep(c(1, 2, 3), each = 4)
IV1 <- rep(c(2, 1), each = 2, times = 3)
IV2 <- rep(c(1, 2), times = 6)
DV <- c(4,3,2,4,4,4,5,4,4,4,3,5)
M1 <- c(1,1,1,3,3,5,3,4,3,4,2,4)
M2 <- c(3,4,4,4,4,5,4,5,4,4,4,4)
df2 <- data.frame(id, IV1, IV2, DV, M1, M2)
df2
#>    id IV1 IV2 DV M1 M2
#> 1   1   2   1  4  1  3
#> 2   1   2   2  3  1  4
#> 3   1   1   1  2  1  4
#> 4   1   1   2  4  3  4
#> 5   2   2   1  4  3  4
#> 6   2   2   2  4  5  5
#> 7   2   1   1  5  3  4
#> 8   2   1   2  4  4  5
#> 9   3   2   1  4  3  4
#> 10  3   2   2  4  4  4
#> 11  3   1   1  3  2  4
#> 12  3   1   2  5  4  4

model <- "
  # Regression paths from IV to mediators
  M1 ~ a1*IV1 + a2*IV2 + a3*IV1:IV2
  M2 ~ b1*IV1 + b2*IV2 + b3*IV1:IV2
  
  # Regression paths from mediators to DV
  DV ~ c1*M1 + c2*M2
  
  # Indirect effects
  indirect1_lo := (a1 + a3*1)*c1 + (b1 + b3*1)*c2
  indirect1_ave := (a1 + a3*1.5)*c1 + (b1 + b3*1.5)*c2
  indirect1_lo := (a1 + a3*2)*c1 + (b1 + b3*2)*c2
  indirect2_lo := (a2 + a3*1)*c1 + (b2 + b3*1)*c2
  indirect2_ave := (a2 + a3*1.5)*c1 + (b2 + b3*1.5)*c2
  indirect2_lo := (a2 + a3*2)*c1 + (b2 + b3*2)*c2
"

fit <- lavaan::cfa(model, data = df2)
fit
#> lavaan 0.6.15 ended normally after 1 iteration
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                        11
#> 
#>   Number of observations                            12
#> 
#> Model Test User Model:
#>                                                       
#>   Test statistic                                 7.563
#>   Degrees of freedom                                 4
#>   P-value (Chi-square)                           0.109

^{Created on 2023-07-03 with reprex v2.0.2}