Updated May 4, 2015 to show how to impose more constraints on the graph by editing the dot code.

This is a quick demo/proof-of-concept of how one may construct a path diagram for a structural equation model. The basic strategy is to convert a path-estimates data-frame from lavaan into node and edge data-frames for DiagrammeR.

First, let’s fit a model using lavaan’s example growth curve analysis.

library("stringr")
library("lavaan")
library("DiagrammeR")
library("dplyr")
library("semPlot")

# Linear growth model with a time-varying covariate
model <- "
  # intercept and slope with fixed coefficients
  i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  # regressions
  i ~ x1 + x2
  s ~ x1 + x2
  # time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
"
fit <- growth(model, data = Demo.growth)

Here’s the simplest way to plot this model in R using the semPlot package.

semPaths(fit, intercept = FALSE, whatLabel = "est",
         residuals = FALSE, exoCov = FALSE)

But for publication-quality graphics, we would want a manually created diagram (typically drawn with Powerpoint, I’m told) or a Graphviz diagram (that we can obsessively tweak till perfect). We’ll use the latter strategy.

Transforming a model summary into a graph description

To create a Graphiviz diagram, we convert the dataframe provided by parameterestimates into node and edge data-frames for DiagrammeR.

# Look at some sample rows
fit %>% parameterestimates %>% sample_n(5)
##   lhs op rhs    est    se     z pvalue ci.lower ci.upper
## 1  c2 ~~  c2  0.899 0.000    NA     NA    0.899    0.899
## 2   s =~  t3  2.000 0.000    NA     NA    2.000    2.000
## 3  x2 ~~  c2 -0.007 0.000    NA     NA   -0.007   -0.007
## 4  t4  ~  c4  0.330 0.058 5.655      0    0.216    0.445
## 5   i ~1      0.580 0.062 9.368      0    0.459    0.702
# To describe edges, we need the origin, path-type, endpoint, and label.
paths <- fit %>%
  parameterestimates %>%
  select(lhs, op, rhs, est)

# Latent variables are left-hand side of "=~" lines
latent <- paths %>%
  filter(op == "=~") %>%
  select(nodes = lhs) %>%
  distinct %>%
  mutate(shape = "circle")

# This is a node data-frame now
latent
##   nodes  shape
## 1     i circle
## 2     s circle
# Manifest variables are not latent variables
`%not_in%` <- Negate(`%in%`)
manifest <- paths %>%
  filter(op != "~1", lhs %not_in% latent$nodes) %>%
  select(nodes = lhs) %>%
  distinct %>%
  mutate(shape = "square")

# Nodes are prepared
node_set <- combine_nodes(latent, manifest)
node_set
##    nodes  shape
## 1      i circle
## 2      s circle
## 3     t1 square
## 4     t2 square
## 5     t3 square
## 6     t4 square
## 7     x1 square
## 8     x2 square
## 9     c1 square
## 10    c2 square
## 11    c3 square
## 12    c4 square
# Edges will be labeled by the parameter estimates
all_paths <- paths %>%
  filter(op != "~1") %>%
  mutate(label = round(est, 2)) %>%
  select(-est)

# Factor loadings are the paths in the "=~" lines
loadings <- all_paths %>%
  filter(op == "=~") %>%
  mutate(edge_from = lhs, edge_to = rhs, style = "dashed") %>%
  select(edge_from, edge_to, style, label)

# This is now an edge dataframe
loadings
##   edge_from edge_to  style label
## 1         i      t1 dashed     1
## 2         i      t2 dashed     1
## 3         i      t3 dashed     1
## 4         i      t4 dashed     1
## 5         s      t1 dashed     0
## 6         s      t2 dashed     1
## 7         s      t3 dashed     2
## 8         s      t4 dashed     3
# Regressions are the paths in the "~" lines
regressions <- all_paths %>%
  filter(op == "~") %>%
  rename(edge_to = lhs, edge_from = rhs) %>%
  mutate(style = "solid") %>%
  select(edge_from, edge_to, style, label)

edge_set <- combine_edges(loadings, regressions)
edge_set
##    edge_from edge_to  style label
## 1          i      t1 dashed     1
## 2          i      t2 dashed     1
## 3          i      t3 dashed     1
## 4          i      t4 dashed     1
## 5          s      t1 dashed     0
## 6          s      t2 dashed     1
## 7          s      t3 dashed     2
## 8          s      t4 dashed     3
## 9         x1       i  solid  0.61
## 10        x2       i  solid   0.6
## 11        x1       s  solid  0.26
## 12        x2       s  solid  0.52
## 13        c1      t1  solid  0.14
## 14        c2      t2  solid  0.29
## 15        c3      t3  solid  0.33
## 16        c4      t4  solid  0.33
# What's missing? Covariance/variance (~~) and intercept (~1) paths

# Combine edges and nodes
my_graph <- graphviz_graph(
  nodes = node_set,
  edges_df = edge_set,
  graph_attrs = c("ranksep = 1"))

# We can plot the graph directly
graphviz_render(my_graph)

The time variables are out of order and the covariates would look best on the bottom, pointing upwards, but we’re off to a good start. We can modify to the dot code to tweak the final plot.

Tweaking the dot code

First let’s extract the dot code from the diagram we just made.

# Copy/edit the dot code directly to make the diagram just right
cat(my_graph$dot_code)
## digraph {
## 
## graph [ranksep = 1]
## 
## 
##   'i' [shape = 'circle'] 
##   's' [shape = 'circle'] 
##   't1' [shape = 'square'] 
##   't2' [shape = 'square'] 
##   't3' [shape = 'square'] 
##   't4' [shape = 'square'] 
##   'x1' [shape = 'square'] 
##   'x2' [shape = 'square'] 
##   'c1' [shape = 'square'] 
##   'c2' [shape = 'square'] 
##   'c3' [shape = 'square'] 
##   'c4' [shape = 'square'] 
## 'i'->'t1' [style = 'dashed', label = '1'] 
## 'i'->'t2' [style = 'dashed', label = '1'] 
## 'i'->'t3' [style = 'dashed', label = '1'] 
## 'i'->'t4' [style = 'dashed', label = '1'] 
## 's'->'t1' [style = 'dashed', label = '0'] 
## 's'->'t2' [style = 'dashed', label = '1'] 
## 's'->'t3' [style = 'dashed', label = '2'] 
## 's'->'t4' [style = 'dashed', label = '3'] 
## 'x1'->'i' [style = 'solid', label = '0.61'] 
## 'x2'->'i' [style = 'solid', label = '0.6'] 
## 'x1'->'s' [style = 'solid', label = '0.26'] 
## 'x2'->'s' [style = 'solid', label = '0.52'] 
## 'c1'->'t1' [style = 'solid', label = '0.14'] 
## 'c2'->'t2' [style = 'solid', label = '0.29'] 
## 'c3'->'t3' [style = 'solid', label = '0.33'] 
## 'c4'->'t4' [style = 'solid', label = '0.33'] 
## }

We will update this plot by adding hidden edges to enforce the time ordering (t1 before t2, etc.). We will also force the covariates to fall in the bottom row. We do this by simple copy/pasting the original code plus some new statements to impose those constraints. For more intricate diagrams, we would create a .gv file and use the edit/preview facilities provided by RStudio.

grViz("
digraph {

  graph [ranksep = 1]
  'i' [shape = 'circle']
  's' [shape = 'circle']
  't1' [shape = 'square']
  't2' [shape = 'square']
  't3' [shape = 'square']
  't4' [shape = 'square']
  'x1' [shape = 'square']
  'x2' [shape = 'square']
  'c1' [shape = 'square']
  'c2' [shape = 'square']
  'c3' [shape = 'square']
  'c4' [shape = 'square']
  'i'->'t1' [style = 'dashed', label = '1']
  'i'->'t2' [style = 'dashed', label = '1']
  'i'->'t3' [style = 'dashed', label = '1']
  'i'->'t4' [style = 'dashed', label = '1']
  's'->'t1' [style = 'dashed', label = '0']
  's'->'t2' [style = 'dashed', label = '1']
  's'->'t3' [style = 'dashed', label = '2']
  's'->'t4' [style = 'dashed', label = '3']
  'x1'->'i' [style = 'solid', label = '0.61']
  'x2'->'i' [style = 'solid', label = '0.6']
  'x1'->'s' [style = 'solid', label = '0.26']
  'x2'->'s' [style = 'solid', label = '0.52']
  'c1'->'t1' [style = 'solid', label = '0.14']
  'c2'->'t2' [style = 'solid', label = '0.29']
  'c3'->'t3' [style = 'solid', label = '0.33']
  'c4'->'t4' [style = 'solid', label = '0.33']

  # Additional constraints on the graph
  t1->t2->t3->t4 [style='invis']
  {rank = 'max'; c1; c2; c3; c4;}
  {rank = 'same'; t1; t2; t3; t4;}
}
")

devtools::session_info()
##  setting  value                       
##  version  R version 3.2.0 (2015-04-16)
##  system   x86_64, mingw32             
##  ui       RTerm                       
##  language (EN)                        
##  collate  English_United States.1252  
##  tz       America/Chicago             
## 
##  package      * version  date       source        
##  acepack      * 1.3-3.3  2013-05-03 CRAN (R 3.1.1)
##  assertthat   * 0.1      2013-12-06 CRAN (R 3.0.2)
##  car          * 2.0-25   2015-03-03 CRAN (R 3.1.2)
##  cluster      * 2.0.1    2015-01-31 CRAN (R 3.2.0)
##  colorspace   * 1.2-6    2015-03-11 CRAN (R 3.1.3)
##  corpcor      * 1.6.7    2014-09-29 CRAN (R 3.1.1)
##  d3Network    * 0.5.2.1  2015-01-31 CRAN (R 3.1.2)
##  DBI          * 0.3.1    2014-09-24 CRAN (R 3.1.1)
##  devtools     * 1.7.0    2015-01-17 CRAN (R 3.1.2)
##  DiagrammeR     0.6      2015-04-30 CRAN (R 3.2.0)
##  digest       * 0.6.8    2014-12-31 CRAN (R 3.1.2)
##  dplyr          0.4.1    2015-01-14 CRAN (R 3.1.2)
##  ellipse      * 0.3-8    2013-04-13 CRAN (R 3.0.2)
##  evaluate     * 0.7      2015-04-21 CRAN (R 3.2.0)
##  fdrtool      * 1.2.14   2015-03-09 CRAN (R 3.1.3)
##  foreign      * 0.8-63   2015-02-20 CRAN (R 3.2.0)
##  formatR      * 1.2      2015-04-21 CRAN (R 3.2.0)
##  Formula      * 1.2-1    2015-04-07 CRAN (R 3.2.0)
##  ggm          * 2.3      2015-01-21 CRAN (R 3.1.2)
##  ggplot2      * 1.0.1    2015-03-17 CRAN (R 3.1.3)
##  glasso       * 1.8      2014-07-22 CRAN (R 3.1.1)
##  gridExtra    * 0.9.1    2012-08-09 CRAN (R 3.1.2)
##  gtable       * 0.1.2    2012-12-05 CRAN (R 3.0.1)
##  gtools       * 3.4.2    2015-04-10 CRAN (R 3.1.3)
##  Hmisc        * 3.16-0   2015-04-30 CRAN (R 3.2.0)
##  htmltools    * 0.2.6    2014-09-08 CRAN (R 3.1.1)
##  htmlwidgets  * 0.3.2    2014-12-09 CRAN (R 3.1.2)
##  huge         * 1.2.6    2014-02-28 CRAN (R 3.0.2)
##  igraph       * 0.7.1    2014-04-22 CRAN (R 3.1.0)
##  jpeg         * 0.1-8    2014-01-23 CRAN (R 3.1.0)
##  knitr        * 1.10     2015-04-23 CRAN (R 3.2.0)
##  lattice      * 0.20-31  2015-03-30 CRAN (R 3.2.0)
##  latticeExtra * 0.6-26   2013-08-15 CRAN (R 3.0.2)
##  lavaan         0.5-18   2015-04-04 CRAN (R 3.1.3)
##  lazyeval     * 0.1.10   2015-01-02 CRAN (R 3.1.2)
##  lisrelToR    * 0.1.4    2013-05-08 CRAN (R 3.1.0)
##  lme4         * 1.1-7    2014-07-03 CRAN (R 3.1.0)
##  magrittr     * 1.5      2014-11-22 CRAN (R 3.1.2)
##  MASS         * 7.3-40   2015-03-21 CRAN (R 3.2.0)
##  Matrix       * 1.2-0    2015-04-04 CRAN (R 3.2.0)
##  matrixcalc   * 1.0-3    2012-09-15 CRAN (R 3.0.0)
##  mgcv         * 1.8-6    2015-03-31 CRAN (R 3.2.0)
##  minqa        * 1.2.4    2014-10-09 CRAN (R 3.1.1)
##  mnormt       * 1.5-2    2015-04-03 CRAN (R 3.1.3)
##  munsell      * 0.4.2    2013-07-11 CRAN (R 3.0.1)
##  nlme         * 3.1-120  2015-02-20 CRAN (R 3.2.0)
##  nloptr       * 1.0.4    2014-08-04 CRAN (R 3.1.1)
##  nnet         * 7.3-9    2015-02-11 CRAN (R 3.2.0)
##  pbivnorm     * 0.6.0    2015-01-23 CRAN (R 3.1.2)
##  pbkrtest     * 0.4-2    2014-11-13 CRAN (R 3.1.2)
##  plyr         * 1.8.2    2015-04-21 CRAN (R 3.2.0)
##  png          * 0.1-7    2013-12-03 CRAN (R 3.0.2)
##  proto        * 0.3-10   2012-12-22 CRAN (R 3.0.0)
##  psych        * 1.5.4    2015-04-27 CRAN (R 3.2.0)
##  qgraph       * 1.3.1    2015-03-10 CRAN (R 3.1.3)
##  quadprog     * 1.5-5    2013-04-17 CRAN (R 3.0.0)
##  quantreg     * 5.11     2015-01-11 CRAN (R 3.1.2)
##  RColorBrewer * 1.1-2    2014-12-07 CRAN (R 3.1.2)
##  Rcpp         * 0.11.6   2015-05-01 CRAN (R 3.2.0)
##  reshape2     * 1.4.1    2014-12-06 CRAN (R 3.1.2)
##  rjson        * 0.2.15   2014-11-03 CRAN (R 3.1.2)
##  RJSONIO      * 1.3-0    2014-07-28 CRAN (R 3.1.1)
##  rmarkdown    * 0.5.1    2015-01-26 CRAN (R 3.1.2)
##  rockchalk    * 1.8.92   2015-02-13 CRAN (R 3.1.2)
##  rpart        * 4.1-9    2015-02-24 CRAN (R 3.2.0)
##  rstudioapi   * 0.3.1    2015-04-07 CRAN (R 3.2.0)
##  scales       * 0.2.4    2014-04-22 CRAN (R 3.1.0)
##  sem          * 3.1-5    2014-08-03 CRAN (R 3.1.1)
##  semPlot        1.0.1    2014-08-15 CRAN (R 3.1.1)
##  sna          * 2.3-2    2014-01-14 CRAN (R 3.0.2)
##  SparseM      * 1.6      2015-01-05 CRAN (R 3.1.2)
##  stringi      * 0.4-1    2014-12-14 CRAN (R 3.1.2)
##  stringr        1.0.0    2015-04-30 CRAN (R 3.2.0)
##  survival     * 2.38-1   2015-02-24 CRAN (R 3.2.0)
##  tables       * 0.7.79   2014-07-24 CRAN (R 3.1.2)
##  whisker      * 0.3-2    2013-04-28 CRAN (R 3.0.1)
##  XML          * 3.98-1.1 2013-06-20 CRAN (R 3.0.1)
##  yaml         * 2.1.13   2014-06-12 CRAN (R 3.1.2)