Manually Plot Path Diagrams for Structural Equation Models

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

Introduction

When plotting path diagrams for structural equation models (SEM) in R, users must choose between automated plotting tools or 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 basic workflow of the package is simple. Users fit their model as a lavaan object, then create a list specifying each node’s position on the diagram using x-y coordinates. Next, users create another list which specifies, for each path, which nodes it connects to and where on the perimeter (top, bottom, left, right) of the relevant nodes the path should begin and end.

The primary function, diyPaths(), processes this information and returns a ggplot2 object of the diagram. It automatically centers estimates on each path’s midpoint and assigns the correct geometric shapes to each variable type. The package includes extensive arguments that facilitate near-complete control for fine-tuning every visual detail.

The package contains four functions: the primary diyPaths() function alongside three helper functions (node(), path(), and panel_title()) designed to simplify argument specification and prevent errors. The most intuitive way to learn how to use diy.sem.plot is through practical examples. Through each, I will demonstrate the basic workflow and explore the arguments available.

Examples

Example 1: Making a Path Diagram for a Simple Path Model

The first step is to specify and fit your model as a lavaan object. Here, we use the mtcars dataset to specify a multiple regression model where vehicle weight (wt) and horsepower (hp) simultaneously predict miles per gallon (mpg), while also estimating the covariance between the two predictors.

data_cars <- mtcars

model_cars <- '
  mpg ~  wt + hp
  wt ~~ hp
'

fit_cars <- sem(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, we need to specify each node’s position on the diagram. For this,we will create an object called node_list and specify it as a list of node objects created using the helper function node(). The node() function takes four key arguments: name (the variable name in your lavaan model), x and y (its coordinates on the diagram) and an optional label that specifies the text actually displayed on the node in the final diagram (defaulting to the variable name if left unspecified).

Since we have not rendered the diagram yet, we will make an initial estimate for node positions and fine-tune them afterwards.

node_list <- list(
  node(name = "wt", x = 1, y = 2, label = "Weight"),
  node(name = "hp", x = 1, y = 1, label = "Horse Power"),
  node(name = "mpg", x = 3, y = 1, label = "Miles Per\n Gallon")
)

Next, we need to specify each path’s position relative to its respective starting and ending nodes. We create an object called path_list and specify it as a list of path specifications created using the path() helper function.

The path() function takes four main arguments alongside four optional fine-tuning arguments. Each 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, you can use nudge_text_x and nudge_text_y to adjust where a path’s parameter estimate appears on the diagram. For covariances or correlations, the numeric cov_curve argument controls the scale of the arced path and the sign of the number dictates direction it bends. By default, variance and residual loops are toggled off (set show_variances to TRUE within diyPaths() to display them). Explicitly defining a path for a node’s variance/residual lets you override its default top position by setting the variance_position argument to the desired location on the node’s perimeter (top, bottom, left, or right).

Again as we have not rendered the diagram yet, we’ll make an initial estimate for path positions.


path_list <- list(
  path(from = "wt", to = "mpg", side_from = "right", side_to = "left"),
  path(from = "hp",  to = "mpg", side_from = "right", side_to = "left"),
  path(from = "wt",  to ="hp", side_from = "left", side_to = "left", cov_curve = -0.4) 
)

Finally, we can render our diagram as a ggplot object using the function diyPaths(). This function only requires three arguments - fit, node_positions, and path_positions - to successfully render the diagram with everything else optional. fit takes the fitted lavaan object of your model, node_positions will take the object containing the list of node() in our case node_list, and path_positions will take the object containing the list of path() in our case path_list.

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

diyPaths(fit = fit_cars, 
         node_positions = node_list, 
         path_positions = path_list, 
         show_grid = TRUE)

Now that we can see our render, it is time to do some fine-tuning. First, most notably, our covariance curve is facing the wrong way. We can fix this by adjusting the sign of the value within the path’s cov_curve argument. For a mostly vertical path, a positive curvature bends it left and a negative curvature bends it right. For a mostly horizontal path, positive curvature bends it down and negative curvature bends it up. We will then use the grid to center the boxes and give them a bit of space. We can then remove the grid by deleting the show_grid argument.

We should add an indication of path significance. diyPaths() has a few options for this. For our plot, we will be maximal and include significance stars, 95% confidence intervals and p-values for the estimates by changing the arguments est_stars, est_ci, and est_p to equal TRUE. We will also change the line type to represent significance using the classic convention of dashed lines for non-significant paths and solid lines for significant paths by changing the argument sig_linetype equal to TRUE.

This is quite a small model, so the default text size is not proportionate to the scale of the boxes. We can do some easy fine-tuning by scaling the path estimates’ text size up via path_text_size, and since the nodes are observed variables, their text size via observed_node_text_size. The lines are a bit thin and the arrowheads are a bit small, so we can adjust them by increasing line_thickness (default = 0.6) and arrow_size (default = 0.2).

Lastly, the plot did not render properly within the margins (i.e. path estimates are cut off at edges). We can fix this by increasing the margins via margin_x and margin_y arguments, which are currently set to 0.5.

node_list <- list(
  node(name = "wt", x = 1, y = 3, label = "Weight"),
  node(name = "hp", x = 1, y = 1, label = "Horse Power"),
  node(name = "mpg", x = 4, y = 2, label = "Miles Per \nGallon"))

path_list <- list(
           path(from = "wt", to = "mpg", side_from = "right", side_to = "left"),
           path(from = "hp",  to = "mpg", side_from = "right", side_to = "left"),
           path(from = "wt",  to ="hp", side_from = "left", side_to = "left", cov_curve = 0.4))

diyPaths(fit_cars, 
         node_positions = node_list,
         path_positions = path_list,
         est_stars = TRUE,
         est_ci = TRUE,
         est_p = TRUE,
         sig_linetype = TRUE,
         observed_node_text_size = 8,
         path_text_size = 6,
         line_thickness = 0.8,
         arrow_size = 0.3,
         margin_x = 1,
         margin_y =1
         )

Example 2: Making Path Diagram for a Full SEM

Here, we 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(sem_model, data = HolzingerSwineford1939)

Next, we will draft our path diagram with rough estimations of where each node should go. As it is classic convention in SEM to have observed variables smaller than latent, we 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, we 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. We will also enable show_grid to be TRUE and set grid_axis_scale to a finer 0.5 rather than 1, giving us a more precise reference to plot against.

For this diagram I want to use standardised estimates, so I set the argument standardised to TRUE. I also want to see the variances and residuals, so I set the argument show_variances to TRUE, which automatically applies them to all nodes.


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

   # observed variables that visual perception ability loads onto
   node("x1", x = 0.5, y = 0, label = "Visual\nPerception"),
   node("x2", x = 1, y = 0, label = "Cubes"),
   node("x3", x = 1.5, y = 0, label = "Lozenges"),

   # observed variables that textual ability loads onto
   node("x4", x = 0.5, y = 3.5, label = "Paragraph\nComprehension"),
   node("x5", x = 1, y = 3.5, label = "Sentence\nCompletion"),
   node("x6", x = 1.5, y = 3.5, label = "Word\nMeaning"),

   # observed variables that speeded cognitive processing loads onto
   node("x7", x = 5.5, y = 1, label = "Speeded\nAddition"),
   node("x8", x = 5.5, y = 1.75, label = "Speeded\nCounting"),
   node("x9", x = 5.5, y = 2.5, label = "Speeded\nDiscrimination")
)

# Specify the paths

path_list <- list(
  #Structural 
    path(from = "visual",  to = "speed", side_from = "right", side_to = "left"),
   path(from = "textual", to = "speed", side_from = "right", side_to = "left"),

   path(from = "visual",  to = "textual", side_from = "left", side_to = "left", cov_curve = -0.6),
  
   #visual loadings
   path(from = "visual", to = "x1", side_from = "bottom", side_to = "top"), 
   path(from = "visual", to = "x2", side_from = "bottom", side_to = "top"),
   path(from = "visual", to = "x3", side_from = "bottom", side_to = "top"),

   #textual loadings
   path(from = "textual", to = "x4", side_from = "top", side_to = "bottom"),
   path(from = "textual", to = "x5", side_from = "top", side_to = "bottom"),
   path(from = "textual", to = "x6", side_from = "top", side_to = "bottom"),

   #speed loadings
   path(from = "speed",   to = "x7", side_from = "right", side_to = "left"),
   path(from = "speed",   to = "x8", side_from = "right", side_to = "left"),
   path(from = "speed",   to = "x9", side_from = "right", side_to = "left")
)

# creating the diagram

p <- diyPaths(
   fit = fit_sem,
   node_positions = node_list,
   path_positions = path_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, we can see the boxes need at least 0.9 x units of space, not 0.5. We will also give the observed variables more space from their latent variables, and increase the spacing between the textual and visual latent variables. Lastly, variance loops are also awkwardly placed, as show_variances defaults them to the top of the node. We can manually adjust this by specifying the variance loops we want to change as a path within path_list, and setting each variance loop’s variance_position to the side of the node (top, bottom, left, or right) we want it to sit on.

As with the previous example, now that the layout is fixed, we can move on to the final touches: adding an indication of statistical significance, adjusting text sizes, and removing the grid overlay.


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

   # observed variables that visual perception ability loads onto
   node("x1", x = -0.2, y = -0.5, label = "Visual\nPerception"),
   node("x2", x = 1, y = -0.5, label = "Cubes"),
   node("x3", x = 2.2, y = -0.5, label = "Lozenges"),

   # observed variables that textual ability loads onto
   node("x4", x = -0.2, y = 4.5, label = "Paragraph\nComprehension"),
   node("x5", x = 1, y = 4.5, label = "Sentence\nCompletion"),
   node("x6", x = 2.2, y = 4.5, label = "Word\nMeaning"),

   # observed variables that speeded cognitive processing loads onto
   node("x7", x = 6, y = 0.8, label = "Speeded\nAddition"),
   node("x8", x = 6, y = 2, label = "Speeded\nCounting"),
   node("x9", x = 6, y = 3.2, label = "Speeded\nDiscrimination")
)

# Specify the paths

path_list <- list(
  # Structural
  path(from = "visual",  to = "speed",   side_from = "right", side_to = "left"),
  path(from = "textual", to = "speed",   side_from = "right", side_to = "left"),
  path(from = "visual",  to = "textual", side_from = "left",  side_to = "left", cov_curve = -0.6),

  # Visual loadings
  path(from = "visual", to = "x1", side_from = "bottom", side_to = "top"),
  path(from = "visual", to = "x2", side_from = "bottom", side_to = "top"),
  path(from = "visual", to = "x3", side_from = "bottom", side_to = "top"),

  # Textual loadings
  path(from = "textual", to = "x4", side_from = "top", side_to = "bottom"),
  path(from = "textual", to = "x5", side_from = "top", side_to = "bottom"),
  path(from = "textual", to = "x6", side_from = "top", side_to = "bottom"),

  # Speed loadings
  path(from = "speed", to = "x7", side_from = "right", side_to = "left"),
  path(from = "speed", to = "x8", side_from = "right", side_to = "left"),
  path(from = "speed", to = "x9", side_from = "right", side_to = "left"),

  # Latent variance/residual
  path(from = "visual",  to = "visual",  variance_position = "top"),
  path(from = "textual", to = "textual", variance_position = "bottom"),
  path(from = "speed",   to = "speed",   variance_position = "top"),

  # Measurement variances/residuals 
  path(from = "x1", to = "x1", variance_position = "bottom"),
  path(from = "x2", to = "x2", variance_position = "bottom"),
  path(from = "x3", to = "x3", variance_position = "bottom"),

  path(from = "x4", to = "x4", variance_position = "top"),
  path(from = "x5", to = "x5", variance_position = "top"),
  path(from = "x6", to = "x6", variance_position = "top"),

  path(from = "x7", to = "x7", variance_position = "right"),
  path(from = "x8", to = "x8", variance_position = "right"),
  path(from = "x9", to = "x9", variance_position = "right")
)

# creating the diagram

p <- diyPaths(
   fit = fit_sem,
   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
)

print(p)

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

In this example we will explore how to create a multi-group path diagram using diyPaths(). We will use the model we created above 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(sem_model, data = HolzingerSwineford1939, group = "school")

We reuse the diagram structure from Example 2 as is. Rendering this new model produces two SEM diagrams, separated into two panels side by side, with no indication of which school each one represents. We 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.

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

   # observed variables that visual perception ability loads onto
   node("x1", x = -0.2, y = -0.5, label = "Visual\nPerception"),
   node("x2", x = 1, y = -0.5, label = "Cubes"),
   node("x3", x = 2.2, y = -0.5, label = "Lozenges"),

   # observed variables that textual ability loads onto
   node("x4", x = -0.2, y = 4.5, label = "Paragraph\nComprehension"),
   node("x5", x = 1, y = 4.5, label = "Sentence\nCompletion"),
   node("x6", x = 2.2, y = 4.5, label = "Word\nMeaning"),

   # observed variables that speeded cognitive processing loads onto
   node("x7", x = 6, y = 0.8, label = "Speeded\nAddition"),
   node("x8", x = 6, y = 2, label = "Speeded\nCounting"),
   node("x9", x = 6, y = 3.2, label = "Speeded\nDiscrimination")
)

# Specify the paths

path_list <- list(
  # Structural
  path(from = "visual",  to = "speed",   side_from = "right", side_to = "left"),
  path(from = "textual", to = "speed",   side_from = "right", side_to = "left"),
  path(from = "visual",  to = "textual", side_from = "left",  side_to = "left", cov_curve = -0.6),

  # Visual loadings
  path(from = "visual", to = "x1", side_from = "bottom", side_to = "top"),
  path(from = "visual", to = "x2", side_from = "bottom", side_to = "top"),
  path(from = "visual", to = "x3", side_from = "bottom", side_to = "top"),

  # Textual loadings
  path(from = "textual", to = "x4", side_from = "top", side_to = "bottom"),
  path(from = "textual", to = "x5", side_from = "top", side_to = "bottom"),
  path(from = "textual", to = "x6", side_from = "top", side_to = "bottom"),

  # Speed loadings
  path(from = "speed", to = "x7", side_from = "right", side_to = "left"),
  path(from = "speed", to = "x8", side_from = "right", side_to = "left"),
  path(from = "speed", to = "x9", side_from = "right", side_to = "left"),

  # Latent variance/residual
  path(from = "visual",  to = "visual",  variance_position = "top"),
  path(from = "textual", to = "textual", variance_position = "bottom"),
  path(from = "speed",   to = "speed",   variance_position = "top"),

  # Measurement variances/residuals 
  path(from = "x1", to = "x1", variance_position = "bottom"),
  path(from = "x2", to = "x2", variance_position = "bottom"),
  path(from = "x3", to = "x3", variance_position = "bottom"),

  path(from = "x4", to = "x4", variance_position = "top"),
  path(from = "x5", to = "x5", variance_position = "top"),
  path(from = "x6", to = "x6", variance_position = "top"),

  path(from = "x7", to = "x7", variance_position = "right"),
  path(from = "x8", to = "x8", variance_position = "right"),
  path(from = "x9", to = "x9", variance_position = "right")
)

# creating the diagram

p <- diyPaths(
   fit = fit_sem_groups,
   node_positions = node_list,
   path_positions = path_list,
   standardised = TRUE,
   sig_linetype = TRUE,
   observed_node_size_adjust = 0.5, 
   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)

We can use the text rendered by show_group_labels to get information needed to construct proper titles for these diagrams. Diagram titles are specified in a similar way as nodes and paths. We first create an object called title_list, specifying it as a list of titles using the panel_title() helper function. This function takes two arguments: panel_num, which takes the panel number and identifies which diagram we are titling, and title, which is the actual text we want for that title. For models with no grouping variable, the panel argument can be omitted. The title_list object is then passed to the panel_titles argument in diyPaths().

title_list <- list(
  panel_title(panel_num = 1, title = "SEM Diagram: Pasteur School"),
  panel_title(panel_num = 2, title = "SEM Diagram: Grant-White School")
)

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

Currently, the panels are cluttered because they are crammed to fit the width of the document. To fix this, we’ll 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 models, setting this to 1 will give vertical stacking.

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

   # observed variables that visual perception ability loads onto
   node("x1", x = -0.2, y = -0.5, label = "Visual\nPerception"),
   node("x2", x = 1, y = -0.5, label = "Cubes"),
   node("x3", x = 2.2, y = -0.5, label = "Lozenges"),

   # observed variables that textual ability loads onto
   node("x4", x = -0.2, y = 4.5, label = "Paragraph\nComprehension"),
   node("x5", x = 1, y = 4.5, label = "Sentence\nCompletion"),
   node("x6", x = 2.2, y = 4.5, label = "Word\nMeaning"),

   # observed variables that speeded cognitive processing loads onto
   node("x7", x = 6, y = 0.8, label = "Speeded\nAddition"),
   node("x8", x = 6, y = 2, label = "Speeded\nCounting"),
   node("x9", x = 6, y = 3.2, label = "Speeded\nDiscrimination")
)

# Specify the paths

path_list <- list(
  # Structural
  path(from = "visual",  to = "speed",   side_from = "right", side_to = "left"),
  path(from = "textual", to = "speed",   side_from = "right", side_to = "left"),
  path(from = "visual",  to = "textual", side_from = "left",  side_to = "left", cov_curve = -0.6),

  # Visual loadings
  path(from = "visual", to = "x1", side_from = "bottom", side_to = "top"),
  path(from = "visual", to = "x2", side_from = "bottom", side_to = "top"),
  path(from = "visual", to = "x3", side_from = "bottom", side_to = "top"),

  # Textual loadings
  path(from = "textual", to = "x4", side_from = "top", side_to = "bottom"),
  path(from = "textual", to = "x5", side_from = "top", side_to = "bottom"),
  path(from = "textual", to = "x6", side_from = "top", side_to = "bottom"),

  # Speed loadings
  path(from = "speed", to = "x7", side_from = "right", side_to = "left"),
  path(from = "speed", to = "x8", side_from = "right", side_to = "left"),
  path(from = "speed", to = "x9", side_from = "right", side_to = "left"),

  # Latent variance/residual
  path(from = "visual",  to = "visual",  variance_position = "top"),
  path(from = "textual", to = "textual", variance_position = "bottom"),
  path(from = "speed",   to = "speed",   variance_position = "top"),

  # Measurement variances/residuals 
  path(from = "x1", to = "x1", variance_position = "bottom"),
  path(from = "x2", to = "x2", variance_position = "bottom"),
  path(from = "x3", to = "x3", variance_position = "bottom"),

  path(from = "x4", to = "x4", variance_position = "top"),
  path(from = "x5", to = "x5", variance_position = "top"),
  path(from = "x6", to = "x6", variance_position = "top"),

  path(from = "x7", to = "x7", variance_position = "right"),
  path(from = "x8", to = "x8", variance_position = "right"),
  path(from = "x9", to = "x9", variance_position = "right")
)

# creating the diagram

p <- diyPaths(
   fit = fit_sem_groups,
   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,
   path_text_size = 3.5,
   latent_node_text_size = 5,
   est_stars = TRUE,
   est_ci = TRUE,
   show_variances = TRUE,
   panel_cols = 1
)

print(p)

Considerations

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.

Notes & Future Development

diy.sem.plot currently only supports models fitted as lavaan objects. Support for other objects, such as blavaan, is planned. See the NEWS document for updates.