Suggestion: Workaround for `mark_sig()` in Multigroup Models

[jorgesinval]

The function mark_sig() does not directly support multigroup models when using the fit object from lavaan. However, I have developed a workaround: we can use the parameter table, as suggested in the mark_sig() function.

Here’s a reproducible example:

library(lavaan)
#> This is lavaan 0.6-19
#> lavaan is FREE software! Please report any bugs.
library(semptools)
library(semPlot)

# Define a CFA model
model <- "visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9"

# Fit the model with multigroup data
fit <- cfa(model = model, data = HolzingerSwineford1939, group = "school")

# Generate initial SEM plots
plots <- semPaths(object = fit,
                  thresholds = FALSE,
                  whatLabels = "std",
                  intercepts = FALSE,
                  exoCov = TRUE,
                  style = "ram",
                  curvePivot = TRUE,
                  edge.color = "black",
                  DoNotPlot = FALSE) # TRUE produces an error (reported on semPlot's GitHub)

Initial plots for multigroup models:

First Approach: Using mark_sig() Directly

Attempting to use mark_sig() directly results in an error:

#### first approach ####
#error produced
plot_g1 <- mark_sig(plots[[1]], fit) #individual approach; plot by plot
#> Error in `mark_sig()`:
#> ! Multiple-group models are not currently supported.

Second Approach: Workaround with Parameter Estimates

To bypass the limitation, we extract parameter estimates for each group and apply mark_sig() manually:

#### second approach ####
#extracting parameter estimates
par_ests <- parameterEstimates(fit)

n_groups <- length(names(inspect(object = fit)))

# Loop over each group
for (g in 1:n_groups) {
  # Extracting parameter estimates
  par_ests <- parameterEstimates(fit)
  
  # Selecting the ones required for the current group
  par_ests_group <- par_ests[par_ests$group == g, c("lhs", "op", "rhs", "pvalue")]
  
  # Mark significance and store the plot
  plots[[g]] <- mark_sig(semPaths_plot = plots[[g]], ests = par_ests_group)
}

for (i in 1:length(plots)) {
  plot(plots[[i]])
}

Updated plots with marked significance for each group:

^{Created on 2024-12-16 with reprex v2.1.1}

[marklhc]

Hi @jorgesinval, thank you for the suggestion and the nice example. Indeed, the mark_se() function already supported multiple-group models, and @sfcheung and I discussed doing the same for mark_sig(). I've just made a quick update to the mg branch 8e4ba88. You can try the following example by installing

remotes::install_github("sfcheung/semptools", ref = "mg")

library(lavaan)
#> This is lavaan 0.6-19
#> lavaan is FREE software! Please report any bugs.
library(semptools)
library(semPlot)
model <- "visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 "
fit <- cfa(model = model, data = HolzingerSwineford1939, group = "school")
plots <- semPaths(object = fit,
                  thresholds = FALSE,
                  whatLabels = "std",
                  intercepts = FALSE,
                  exoCov = TRUE,
                  style = "ram",
                  curvePivot = TRUE,
                  edge.color = "black")

# Add asterisks
plots2 <- mark_sig(plots, fit)
lapply(plots2, plot)

#> [[1]]
#> [[1]]$mar
#> [1] 0 0 0 0
#> 
#> [[1]]$bg
#> [1] "white"
#> 
#> 
#> [[2]]
#> [[2]]$mar
#> [1] 0 0 0 0
#> 
#> [[2]]$bg
#> [1] "white"

^{Created on 2024-12-15 with reprex v2.1.0}

What do you think, @sfcheung?

[jorgesinval]

Hi @marklhc! Confirmed, it is working! Nice approach. Let's hope that the semPath() function allows the argument DoNotPlot = F for multigroup objects in the future.

[sfcheung]

Thank @jorgesinval for opening this issue and @marklhc for implementing the support to mark_sig(). The main branch has been updated to incorporate @marklhc 's change:

remotes::install_github("sfcheung/semptools")

I hope we will soon have time to add support for multigroup models to other functions.

P.S.: @marklhc , I added a short test for multiple group models to these two functions.

[jorgesinval]

@sfcheung, you're most definitely welcome.

I think that pretty much all functions we might want to use with multigroup work. One example:

library(lavaan)
#> This is lavaan 0.6-19
#> lavaan is FREE software! Please report any bugs.
library(semPlot)
library(semptools)

model <- "
#measurement model
Vis =~ x1 + x2 + x3
Tex =~ x4 + x5 + x6

x7 ~ Vis + Tex

x8 ~ x7 + Vis + Tex
"

vars_col <- c("Vis","Tex","x7","x8")

latent_layout <- matrix(c(NA, "x7",  NA,
                          "Vis",NA, "x8",
                          "Tex",NA, NA),
                        byrow = TRUE,
                        ncol = 3,
                        nrow = 3)


plot_model <- semPlotModel_lavaanModel(model,auto.var = TRUE,
                                       auto.fix.first = F, #change it to avoid dashed line between first- and second-order latent variables (in hierarchical models)
                                       auto.cov.lv.x = TRUE)

plot_model <- keep_nodes(object = plot_model,
                         nodes = vars_col)


diagram_plot <- semPlot::semPaths(object = plot_model,
                                  thresholds = F,
                                  whatLabels = "none",
                                  layout = latent_layout,
                                  intercepts = F,
                                  exoCov = T,
                                  style = "ram",
                                  curvePivot=F,
                                  edge.color="black",
                                  DoNotPlot = T,
                                  residuals = F,
                                  edgeLabels = c(expression(H[1]),
                                                 expression(H[2]),
                                                 expression(H[3]),
                                                 expression(H[4]),
                                                 expression(H[5]),
                                                 ""),
                                  sizeMan = 8)

my_curve_list <- list(list(from = "Vis",
                           to = "Tex",
                           new_curve = -1))

diagram_plot <- set_curve(diagram_plot, my_curve_list)

my_position_list <- c("x7 ~ Vis" = .40,
                      "x7 ~ Tex" = .40)
diagram_plot <- set_edge_label_position(diagram_plot, my_position_list)

order_of_appearance <- match(tolower(vars_col),
                             tolower(diagram_plot$Arguments$labels))

model_latent_vars <- c("Visual",
                       "Textual",
                       "Index I",
                       "Index II")

model_plot_names <- gsub(" ", "\n", model_latent_vars[order_of_appearance])
diagram_plot_long_names <- change_node_label(semPaths_plot = diagram_plot,
                                             label_list = setNames(object = model_plot_names, #change if single indicator variables are added
                                                                   nm = diagram_plot$Arguments$labels),
                                             label.cex = 1.0)

plot(diagram_plot_long_names)

##### Modified Model

model <- "
#measurement model
Vis =~ x1 + x2 + x3
Tex =~ x4 + x5 + x6

x4 ~~ x6

x7 ~ Vis + Tex

x8 ~ x7 + Vis + Tex
"

fit_model <- sem(m=model,
                 d=HolzingerSwineford1939,
                 group = "school",
                 group.equal = c("loadings", "intercepts"))
#> Warning: lavaan->lav_object_post_check():  
#>    some estimated ov variances are negative


vars_names <- c("", paste0(model_latent_vars, " (",1:length(model_latent_vars),")"))

plot_model <- semPlotModel(fit_model)
plot_model <- keep_nodes(object = plot_model,
                                    nodes = vars_col)

plot_sem <- semPlot::semPaths(object = plot_model,
                     thresholds = F,
                     whatLabels = "std",
                     layout = latent_layout,
                     intercepts = F,
                     exoCov = T,
                     style = "ram",
                     curvePivot=TRUE,
                     edge.color="black",
                     DoNotPlot = F)

#based on the number of unique groups
n_groups <- length(names(inspect(object = fit_model)))


# Add asterisks
plot_sem <- mark_sig(plot_sem, fit_model)

#covariance path angle
plot_sem <- lapply(plot_sem, function(x) set_curve(x, my_curve_list))

#edge label position
plot_sem <- lapply(plot_sem, function(x) set_edge_label_position(x, my_position_list))

my_rotate_resid_list <- c("Vis" = 0,
                          "Tex" = 0,
                          "x7" = 0,
                          "x8" = 0)

#residuals rotation
plot_sem <- lapply(plot_sem, function(x) rotate_resid(x, my_rotate_resid_list))

#solution for multiple plots with loop
for (g in 1:n_groups) {
  plot_sem[[g]] <- change_node_label(plot_sem[[g]],
                         setNames(model_plot_names,
                                  plot_sem[[g]]$Arguments$labels),
                         label.cex = 1.0)
}


lapply(plot_sem, plot)

#> [[1]]
#> [[1]]$mar
#> [1] 0 0 0 0
#> 
#> [[1]]$bg
#> [1] "white"
#> 
#> 
#> [[2]]
#> [[2]]$mar
#> [1] 0 0 0 0
#> 
#> [[2]]$bg
#> [1] "white"

^{Created on 2024-12-17 with reprex v2.1.1}