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.
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.
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)