Manually Plot Path Diagrams for Structural Equation Models

Introduction

When plotting path diagrams for structural equation models (SEM) in R, users must choose between automated plotting tools and drawing diagrams entirely by hand. While automated SEM plotting tools are fast, they lack granular customisability and are often not suitable for use as a final product. Conversely, doing this work manually in an external drawing programme is tedious for large models and is completely disconnected from R, meaning estimates must be individually adjusted if the data or model changes.

diy.sem.plot bridges these approaches to give users the best of both worlds, enabling manual specification of the layout and design while handling the rendering automatically.

The package contains six functions: the primary function diyPaths(), four helpers (node(), path(), panel_title(), and panel_note()) that simplify specifying its arguments, and template(), which generates starting code from your fitted model.

The diyPaths() function takes your fitted lavaan model along with specifications of the layout, and returns a ggplot2 object of the diagram. It automatically centres estimates on each path’s midpoint and assigns the correct geometric shapes to each variable type. It includes extensive arguments that facilitate near-complete control for fine-tuning every visual detail.

The basic workflow of the package is simple:

  1. Fit your model as a lavaan object.
  2. Run template() on your fitted model. It prints code for a list specifying each node, a list specifying each path, lists specifying panel titles and notes, and a diyPaths() call, with every variable name and path in your model already filled in.
  3. Edit the node list to specify each node’s position on the diagram using x-y coordinates.
  4. Edit the path list to specify, for each path, where on the perimeter (top, bottom, left, right) of the relevant nodes a path should begin and end.
  5. Edit the title and note lists to give each panel of the diagram a title and, optionally, a note. Models without a grouping variable have a single panel, while multi-group models have one panel per group.
  6. Render the diagram with diyPaths(), then fine-tune and re-render until the layout looks right.

You can skip step 2 and write the lists by hand, but for larger models, generating them with template() saves a lot of time and avoids syntax errors.

Below is a series of fully worked examples. Through each, I will demonstrate this workflow and explore the arguments available.

Examples

Example 1: Making a Path Diagram for a Mediation Model

library(diy.sem.plot)
library(ggplot2)
library(lavaan)
#> This is lavaan 0.7-2
#> lavaan is FREE software! Please report any bugs.

The first step is to specify and fit the model as a lavaan object. Here, I use the mtcars dataset to specify a simple mediation model, where vehicle weight (wt) mediates the effect of horsepower (hp) on miles per gallon (mpg). I also include engine displacement (disp) as a covariate predicting mpg, and allow it to covary with horsepower.

data_cars <- mtcars


model_cars <- '
  wt  ~ a * hp
  mpg ~ b * wt + c * hp + disp
  disp ~~ hp

  indirect := a * b
  total := c + a * b
'

fit_cars <- sem(model = model_cars, data = data_cars)
#> Warning: lavaan->lav_data_full():  
#>    some observed variances are (at least) a factor 1000 times larger than 
#>    others; use varTable(fit) to investigate

Next, I need to specify the layout of the diagram. Rather than writing this from scratch, I will use template(). It takes the fitted model and prints ready-to-run starter code, with every variable name and path in the model already filled in. The code is organised into labelled sections, with comments grouping nodes by variable type and paths by type (regression, loading, covariance, and variance), making larger models easier to navigate. Since diyPaths() relies on these names matching the lavaan model exactly, generating them this way saves time and prevents errors. All I need to do is copy the output into the script and edit the values on each line marked #EDIT ME!.

template() has a few arguments that control how much code it prints. For example, setting include_variance_paths = FALSE leaves out the code for individual variance path customisation and specifying diyPaths_arguments = "min" prints only the arguments diyPaths() needs to run, while "max" prints every available argument with its default value.

Below is what the raw output of template() looks like for the fitted model. In the later examples, I will omit it to save space by setting the argument quiet to TRUE.

template(fit = fit_cars, include_variance_paths = FALSE, diyPaths_arguments = "min", quiet = FALSE)
#### BEGIN diyPaths CODE ####

node_list <- list(
    #node(s) for observed variable(s)
    node(name = "wt", x = 0, y = 0, label = "wt"), #EDIT ME!
    node(name = "mpg", x = 0, y = 0, label = "mpg"), #EDIT ME!
    node(name = "hp", x = 0, y = 0, label = "hp"), #EDIT ME!
    node(name = "disp", x = 0, y = 0, label = "disp") #EDIT ME!
)


path_list <- list(
    # regression path(s)
    path(from = "hp", to = "wt", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0), #EDIT ME!
    path(from = "wt", to = "mpg", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0), #EDIT ME!
    path(from = "hp", to = "mpg", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0), #EDIT ME!
    path(from = "disp", to = "mpg", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0), #EDIT ME!


    # covariance/correlation path(s)
    path(from = "disp", to = "hp", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0, cov_curve = 0.4) #EDIT ME!
)

title_list <- list(
    panel_title(panel_num = 1, title = "My SEM Plot") #EDIT ME!
)

note_list <- list(
    panel_note(panel_num = 1, note_text = NULL) #EDIT ME!
)

my_path_diagram <- diyPaths(fit = fit_cars,
                            node_positions = node_list,
                            path_positions = path_list,
                            panel_titles = title_list,
                            panel_notes = note_list)
my_path_diagram
#### END diyPaths CODE :) ####

The template is made up of five components: a node list, a path list, a title list, a note list, and the diyPaths() call that renders the diagram. In this example, I will break down each component in turn.

The first component is the node list, which specifies where each node sits on the diagram. Each node is created using the helper function node(), which takes four arguments: name (the variable name in the lavaan model), x and y (its coordinates on the diagram), and an optional label that specifies the text displayed on the node in the final diagram, defaulting to the variable name if left unspecified.

template() sets every node to x = 0, y = 0, so rendering it as-is would stack every node on top of one another. Since I have not rendered the diagram yet, I will make an initial estimate for each node’s position and fine-tune them afterwards. I will also give each node a more readable label, using \n to add a line break where a label is long.

node_list <- list(
  #node(s) for observed variable(s)
  node(name = "wt", x = 2, y = 2.5, label = "Weight"),
  node(name = "mpg", x = 4, y = 1, label = "Miles Per\nGallon"),
  node(name = "disp", x = 0, y = -0.5, label = "Engine\nDisplacement"),
  node(name = "hp", x = 0, y = 1, label = "Horsepower")
)

The second component is the path list, which specifies where each path begins and ends. Each path is created using the helper function path(). A path requires a starting node (from) and an ending node (to), matched to their lavaan variable names. The side_from argument dictates which perimeter edge (top, bottom, left, or right) the path anchors to on the starting node, and side_to dictates the same for the ending node.

For further control, nudge_text_x and nudge_text_y adjust where a path’s estimate appears along the x and y axes. For covariances or correlations, the numeric cov_curve argument controls how strongly the path arcs, and its sign dictates which way it bends. Variance and residual paths, which I left out of this example, also take a variance_position argument that sets which side of the node (top, bottom, left, or right) the loop sits on.

As with the nodes, I will make an initial estimate for each path’s position, then fine-tune after rendering.

path_list <- list(
  # regression path(s)
  path(from = "hp", to = "wt", side_from = "top", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
  path(from = "wt", to = "mpg", side_from = "right", side_to = "top", nudge_text_x = 0, nudge_text_y = 0),
  path(from = "hp", to = "mpg", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
  path(from = "disp", to = "mpg", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),

  # covariance/correlation path(s)
  path(from = "disp", to = "hp", side_from = "left", side_to = "left", nudge_text_x = 0, nudge_text_y = 0, cov_curve = -0.4)
)

The next components are the title list and the note list. Titles are created using the panel_title() helper function, which takes two arguments: panel_num, identifying which panel of the diagram the title belongs to, and title, the text to display. Since this model has no grouping variable, the diagram has a single panel, so panel_num is always 1. Notes are created using the panel_note() helper function, which works the same way but takes note_text, the text displayed as a caption beneath the panel. Both titles and notes support basic Markdown and HTML formatting, such as **bold**, *italics*, and <br> for line breaks.

I will title the diagram and use the note to report the indirect effect I defined in the model.

title_list <- list(
      panel_title(panel_num = 1, title = "**Figure 1**<br>*Mediation of Horsepower on MPG by Weight*")
)

note_list <- list(
  panel_note(panel_num = 1, note_text = "*Note.* Indirect effect = -0.036 (SE = 0.009, p < .001)")
)

The final component is the diyPaths() call, which renders the diagram as a ggplot object. It only requires three arguments: fit takes the fitted lavaan model, node_positions takes the node list, and path_positions takes the path list. The title and note lists are passed to panel_titles and panel_notes; if omitted, no title or note is shown.

I will now run the first draft render. To help with fine-tuning the layout, I will overlay a grid on the diagram by setting the optional argument show_grid to TRUE.

diyPaths(fit = fit_cars,
         node_positions = node_list,
         path_positions = path_list,
         panel_titles = title_list,
         panel_notes = note_list,
         show_grid = TRUE)

Now that I can see the render, it is time to do some fine-tuning. The layout works well, so I can remove the grid by deleting the show_grid argument. However, the covariance curve bends in the wrong direction; I can fix this by flipping the sign of the path’s cov_curve value.

I should add an indication of path significance. diyPaths() has a few options for this. For the plot, I will be maximal and include significance stars, 95% confidence intervals and p-values for the estimates by setting the arguments est_stars, est_ci, and est_p to TRUE. I will also change the line type to represent significance, using the common convention of dashed lines for non-significant paths and solid lines for significant paths, by setting sig_linetype to TRUE. By default, estimates are displayed to three decimal places, which can be changed via the digits argument.

This is quite a small model, so the default text size is not proportionate to the scale of the boxes. I can scale up the path estimates’ text size via path_text_size, and, since all the nodes are observed variables, their text size via observed_node_text_size. Likewise, the title and note use ggplot2’s default text sizes, which look small next to the diagram, so I can increase them via title_text_size and note_text_size. The lines are a bit thin and the arrowheads are a bit small, so I can increase line_thickness (default = 0.6) and arrow_size (default = 0.2). I can also specify the fonts used in the diagram. The font_family argument sets the font for all text in the diagram, while title_note_font_family sets a separate font for the title and note.

Lastly, the plot did not render properly within the margins (e.g. path estimates are cut off at the edges). I can fix this by increasing the horizontal margin via the margin_x argument, which defaults to 0.5.

path_list <- list(
  # regression path(s)
  path(from = "hp", to = "wt", side_from = "top", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
  path(from = "wt", to = "mpg", side_from = "right", side_to = "top", nudge_text_x = 0, nudge_text_y = 0),
  path(from = "hp", to = "mpg", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
  path(from = "disp", to = "mpg", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),

  # covariance/correlation path(s)
     path(from = "disp", to = "hp", side_from = "left", side_to = "left", nudge_text_x = 0, nudge_text_y = 0, cov_curve = 0.4))

diyPaths(fit = fit_cars,
         node_positions = node_list,
         path_positions = path_list,
         panel_titles = title_list,
         panel_notes = note_list,
         est_stars = TRUE,
         est_ci = TRUE,
         est_p = TRUE,
         sig_linetype = TRUE,
         observed_node_text_size = 6,
         path_text_size = 4,
         title_text_size = 18,
         note_text_size = 14,
         font_family = "sans",
         title_note_font_family = "serif",
         line_thickness = 0.8,
         arrow_size = 0.3,
         margin_x = 1,
         margin_y = 0.5)

Example 2: Making a Path Diagram for a Full SEM

Here, I use the HolzingerSwineford1939 dataset from the lavaan package to specify a full SEM. The model has three latent variables, visual perception ability, textual ability, and speeded cognitive processing, each measured by three observed variables. Structurally, speeded cognitive processing is regressed on both visual perception ability and textual ability, with the latter two allowed to covary.

sem_model <- '
   visual  =~ x1 + x2 + x3
   textual =~ x4 + x5 + x6
   speed   =~ x7 + x8 + x9

   speed ~ visual + textual
   visual ~~ textual
'

# fit the model

fit_sem <- sem(model = sem_model, data = HolzingerSwineford1939)

As before, I start by running template(fit = fit_sem, include_variance_paths = TRUE, quiet = TRUE, diyPaths_arguments = "min") and editing its output. The raw output is omitted here to save space. As I want to include variance loops, I will set include_variance_paths to TRUE to give me access to the specific paths for later fine-tuning.

template(fit = fit_sem, include_variance_paths = TRUE, quiet = TRUE, diyPaths_arguments = "min")

Next, I will draft the path diagram with rough estimations of where each node should go. As it is common convention in SEM to draw observed variables smaller than latent, I will scale them smaller by setting observed_node_size_adjust in diyPaths() to be a decimal. To see the exact width and height of these scaled boxes, I set look_up_table to TRUE, which returns the width and height of latent and observed nodes on the x and y scale, along with the text sizes used for latent, observed, and path labels. I will also set show_grid to TRUE and grid_axis_scale to a finer 0.5 rather than 1, giving a more precise reference to plot against.

For this diagram, I want to use standardised estimates and include the variance loops, so I set both standardised and show_variances to TRUE.

I will also give the diagram a title, following the same format as Example 1. Since I don’t need a note, I do not specify a note_list and simply leave panel_notes out of the call.

node_list <- list(
   #node(s) for latent variable(s)
   node(name = "visual", x = 1, y = 1, label = "Visual"),
   node(name = "textual", x = 1, y = 2.5, label = "Textual"),
   node(name = "speed", x = 4, y = 1.75, label = "Speed"),

   #node(s) for observed variables loading onto latent variable visual
   node(name = "x1", x = 0.5, y = 0, label = "Visual\nPerception"),
   node(name = "x2", x = 1, y = 0, label = "Cubes"),
   node(name = "x3", x = 1.5, y = 0, label = "Lozenges"),

   #node(s) for observed variables loading onto latent variable textual
   node(name = "x4", x = 0.5, y = 3.5, label = "Paragraph\nComprehension"),
   node(name = "x5", x = 1, y = 3.5, label = "Sentence\nCompletion"),
   node(name = "x6", x = 1.5, y = 3.5, label = "Word\nMeaning"),

   #node(s) for observed variables loading onto latent variable speed
   node(name = "x7", x = 5.5, y = 1, label = "Speeded\nAddition"),
   node(name = "x8", x = 5.5, y = 1.75, label = "Speeded\nCounting"),
   node(name = "x9", x = 5.5, y = 2.5, label = "Speeded\nDiscrimination")
)

path_list <- list(
   # regression path(s)
   path(from = "visual", to = "speed", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "textual", to = "speed", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),

   #loading path(s) from visual
   path(from = "visual", to = "x1", side_from = "bottom", side_to = "top", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "visual", to = "x2", side_from = "bottom", side_to = "top", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "visual", to = "x3", side_from = "bottom", side_to = "top", nudge_text_x = 0, nudge_text_y = 0),

   #loading path(s) from textual
   path(from = "textual", to = "x4", side_from = "top", side_to = "bottom", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "textual", to = "x5", side_from = "top", side_to = "bottom", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "textual", to = "x6", side_from = "top", side_to = "bottom", nudge_text_x = 0, nudge_text_y = 0),

   #loading path(s) from speed
   path(from = "speed", to = "x7", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "speed", to = "x8", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "speed", to = "x9", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),

   # covariance/correlation path(s)
   path(from = "textual", to = "visual", side_from = "left", side_to = "left", nudge_text_x = 0, nudge_text_y = 0, cov_curve = -0.6),

   # variance path(s)
   path(from = "x1", to = "x1", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x2", to = "x2", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x3", to = "x3", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x4", to = "x4", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x5", to = "x5", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x6", to = "x6", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x7", to = "x7", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x8", to = "x8", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x9", to = "x9", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "visual", to = "visual", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "textual", to = "textual", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "speed", to = "speed", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top")
)


title_list <- list(
   panel_title(panel_num = 1, title = "**Figure 2**<br>*Visual and Textual Ability Predicting Speeded Cognitive Processing*")
)

p <- diyPaths(
   fit = fit_sem,
   node_positions = node_list,
   path_positions = path_list,
   panel_titles = title_list,
   standardised = TRUE,
   observed_node_size_adjust = 0.6,
   show_grid = TRUE,
   grid_axis_scale = 0.5,
   look_up_table = TRUE,
   show_variances = TRUE
)

print(p)
#> $plot

#> 
#> $look_up_table
#>       type width height text_size
#> 1   latent   1.5    1.0       4.0
#> 2 observed   0.9    0.6       4.0
#> 3     path    NA     NA       3.5

This looks decent for a first draft. However, the observed variable nodes for textual and visual clearly do not have enough room on the x-axis. Looking at the grid and look-up table, I can see the boxes need at least 0.9 x units of space, not 0.5. I will also give the observed variables more space from their latent variables, and increase the spacing between the textual and visual latent variables. The estimates on the outer loading paths sit quite close together, so I nudge them apart along the x-axis using nudge_text_x.

Lastly, the variance loops are awkwardly placed, as the template positions every loop at the top of its node. I can fix this by editing each variance path’s variance_position to the side of the node (top, bottom, left, or right) I want it to sit on.

As with the previous example, now that the layout is fixed, I can move on to the final touches: adding an indication of statistical significance, adjusting text sizes, removing the grid overlay, and applying the same title size and fonts as Example 1.

node_list <- list(
   #node(s) for latent variable(s)
   node(name = "visual", x = 1, y = 1, label = "Visual"),
   node(name = "textual", x = 1, y = 3, label = "Textual"),
   node(name = "speed", x = 4, y = 2, label = "Speed"),

   #node(s) for observed variables loading onto latent variable visual
   node(name = "x1", x = -0.2, y = -0.5, label = "Visual\nPerception"),
   node(name = "x2", x = 1, y = -0.5, label = "Cubes"),
   node(name = "x3", x = 2.2, y = -0.5, label = "Lozenges"),

   #node(s) for observed variables loading onto latent variable textual
   node(name = "x4", x = -0.2, y = 4.5, label = "Paragraph\nComprehension"),
   node(name = "x5", x = 1, y = 4.5, label = "Sentence\nCompletion"),
   node(name = "x6", x = 2.2, y = 4.5, label = "Word\nMeaning"),

   #node(s) for observed variables loading onto latent variable speed
   node(name = "x7", x = 6, y = 0.8, label = "Speeded\nAddition"),
   node(name = "x8", x = 6, y = 2, label = "Speeded\nCounting"),
   node(name = "x9", x = 6, y = 3.2, label = "Speeded\nDiscrimination")
)

path_list <- list(
   # regression path(s)
   path(from = "visual", to = "speed", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "textual", to = "speed", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),

   #loading path(s) from visual
   path(from = "visual", to = "x1", side_from = "bottom", side_to = "top", nudge_text_x = -0.1, nudge_text_y = 0),
   path(from = "visual", to = "x2", side_from = "bottom", side_to = "top", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "visual", to = "x3", side_from = "bottom", side_to = "top", nudge_text_x = 0.1, nudge_text_y = 0),

   #loading path(s) from textual
   path(from = "textual", to = "x4", side_from = "top", side_to = "bottom", nudge_text_x = -0.1, nudge_text_y = 0),
   path(from = "textual", to = "x5", side_from = "top", side_to = "bottom", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "textual", to = "x6", side_from = "top", side_to = "bottom", nudge_text_x = 0.1, nudge_text_y = 0),

   #loading path(s) from speed
   path(from = "speed", to = "x7", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "speed", to = "x8", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),
   path(from = "speed", to = "x9", side_from = "right", side_to = "left", nudge_text_x = 0, nudge_text_y = 0),

   # covariance/correlation path(s)
   path(from = "textual", to = "visual", side_from = "left", side_to = "left", nudge_text_x = 0, nudge_text_y = 0, cov_curve = -0.6),

   # variance path(s)
   path(from = "x1", to = "x1", nudge_text_x = 0, nudge_text_y = 0, variance_position = "bottom"),
   path(from = "x2", to = "x2", nudge_text_x = 0, nudge_text_y = 0, variance_position = "bottom"),
   path(from = "x3", to = "x3", nudge_text_x = 0, nudge_text_y = 0, variance_position = "bottom"),
   path(from = "x4", to = "x4", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x5", to = "x5", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x6", to = "x6", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "x7", to = "x7", nudge_text_x = 0, nudge_text_y = 0, variance_position = "right"),
   path(from = "x8", to = "x8", nudge_text_x = 0, nudge_text_y = 0, variance_position = "right"),
   path(from = "x9", to = "x9", nudge_text_x = 0, nudge_text_y = 0, variance_position = "right"),
   path(from = "visual", to = "visual", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top"),
   path(from = "textual", to = "textual", nudge_text_x = 0, nudge_text_y = 0, variance_position = "bottom"),
   path(from = "speed", to = "speed", nudge_text_x = 0, nudge_text_y = 0, variance_position = "top")
)

title_list <- list(
   panel_title(panel_num = 1, title = "**Figure 2**<br>*Visual and Textual Ability Predicting Speeded Cognitive Processing*")
)

p <- diyPaths(
   fit = fit_sem,
   node_positions = node_list,
   path_positions = path_list,
   panel_titles = title_list,
   standardised = TRUE,
   sig_linetype = TRUE,
   observed_node_size_adjust = 0.6,
   observed_node_text_size = 3.5,
   latent_node_text_size = 5,
   est_stars = TRUE,
   est_ci = TRUE,
   show_variances = TRUE,
   title_text_size = 16,
   font_family = "sans",
   title_note_font_family = "serif"
)

print(p)

Example 3: Making a Path Diagram for a Multi-Group SEM

In this example, I will demonstrate how to create a multi-group path diagram using diyPaths(). I will use the model from Example 2 with the HolzingerSwineford1939 dataset, but this time fit it separately by the school the children attended, either “Grant-White” or “Pasteur”.

fit_sem_groups <- sem(model = sem_model, data = HolzingerSwineford1939, group = "school")

Since the model is the same, running template(fit = fit_sem_groups) produces the same node and path lists as Example 2, so I can reuse node_list and path_list as they are. The difference is that, because the model now has two groups, the template prints one panel_title() and one panel_note() line for each panel.

Rendering this model produces two SEM diagrams, separated into two panels side by side, with no indication of which school each one represents. I can figure this out by setting show_group_labels to TRUE, which adds a small label to the bottom corner of each panel, showing its group name and panel number.

p <- diyPaths(
   fit = fit_sem_groups,
   node_positions = node_list,
   path_positions = path_list,
   standardised = TRUE,
   sig_linetype = TRUE,
   observed_node_size_adjust = 0.6,
   observed_node_text_size = 3.5,
   latent_node_text_size = 5,
   est_stars = TRUE,
   est_ci = TRUE,
   show_variances = TRUE,
   show_group_labels = TRUE
)

print(p)

I can use the labels rendered by show_group_labels to fill in the template’s title and note lists. In both panel_title() and panel_note(), panel_num identifies which diagram the title or note belongs to, matching the panel numbers shown in the labels. I will title each panel with its school, and use the notes to report each group’s sample size.

title_list <- list(
  panel_title(panel_num = 1, title = "**Figure 3a**<br>*SEM Diagram: Pasteur School*"),
  panel_title(panel_num = 2, title = "**Figure 3b**<br>*SEM Diagram: Grant-White School*")
)

note_list <- list(
  panel_note(panel_num = 1, note_text = "*Note.* n = 156"),
  panel_note(panel_num = 2, note_text = "*Note.* n = 145")
)

I will apply the standard fine-tuning, as well as set show_group_labels to FALSE, since I now have proper titles instead.

Currently, the panels are cluttered because they are crammed to fit the width of the document. To fix this, I will stack the plots atop each other rather than side by side by adjusting the panel_cols argument, which controls the number of columns used when arranging multi-group panels. For two groups, setting this to 1 will give vertical stacking.

p <- diyPaths(
   fit = fit_sem_groups,
   node_positions = node_list,
   path_positions = path_list,
   panel_titles = title_list,
   panel_notes = note_list,
   standardised = TRUE,
   sig_linetype = TRUE,
   observed_node_size_adjust = 0.6,
   observed_node_text_size = 3.5,
   path_text_size = 3.5,
   latent_node_text_size = 5,
   est_stars = TRUE,
   est_ci = TRUE,
   show_variances = TRUE,
   title_text_size = 16,
   note_text_size = 12,
   font_family = "sans",
   title_note_font_family = "serif",
   panel_cols = 1
)

print(p)

Tips for Rendering

How a diagram renders can vary depending on where you view it, for example in RStudio’s Plots pane, or in knitted output from R Markdown (which can be previewed via Notebook output when running chunks individually). A layout that looks well-proportioned in one context may appear cramped or misaligned in another. Before fine-tuning, I highly recommend staying consistent with the method you use to view your diagram, and making sure it’s rendering at the actual size you need the final image to be.

Viewing your diagram in RStudio’s Plots pane has the advantage that, when you export it, you can specify the exact proportions of the image. However, the pane’s proportion adjusters have fixed limits, so larger or more complex diagrams can end up cramped and difficult to read.

In the examples above, diagrams were rendered by knitting this R Markdown file. To fine-tune the diagrams, I previewed them by running each code chunk individually and viewing the resulting Notebook output. R Markdown has the advantage of letting you manually specify plot size directly in the code chunk, with no fixed limit on the dimensions you can set. For example, the single-group diagrams were specified with {r, fig.width = 10, fig.height = 8, out.width = "100%"}, while the multi-group diagram used {r, fig.width = 10, fig.height = 16, out.width = "100%"}. This made it possible to render larger images for bigger plots without cramping, and the resulting image can also be directly copied or saved.