--- title: "Manually Plot Path Diagrams for Structural Equation Models" output: rmarkdown::html_vignette: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{Manually Plot Path Diagrams for Structural Equation Models} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(diy.sem.plot) library(ggplot2) library(lavaan) ``` ## 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. ```{r} data_cars <- mtcars model_cars <- ' mpg ~ wt + hp wt ~~ hp ' fit_cars <- sem(model_cars, data = data_cars) ``` 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. ```{r} 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. ```{r} 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`. ```{r, fig.width = 10, fig.height = 8,out.width = "100%"} 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. ```{r, fig.width = 10, fig.height = 8,out.width = "100%"} 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. ```{r} 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. ```{r, fig.width = 10, fig.height = 8,out.width = "100%"} 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) ``` 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. ```{r, fig.width = 10, fig.height = 8,out.width = "100%"} 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". ```{r} 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. ```{r, fig.width = 10, fig.height = 8,out.width = "100%"} 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()`. ```{r} 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. ```{r, fig.width = 10, fig.height = 16, out.width = "100%"} 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.