GG SEM example
Thomas V. Pollet
2021-11-11
Background.
Some students occasionally have issues with Semplot due to the .png package. Here we provide an alternative. This is based on this webpage.
Load the packages.
Install the packages below.
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.4 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(tidygraph)##
## Attaching package: 'tidygraph'## The following object is masked from 'package:stats':
##
## filterlibrary(ggraph)
library(lavaan)## This is lavaan 0.6-9
## lavaan is FREE software! Please report any bugs.Sample model.
Model_1 <- ' # Three factors.
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9 'fit <- cfa(Model_1, data=HolzingerSwineford1939)
summary(fit, fit.measures=TRUE)## lavaan 0.6-9 ended normally after 35 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 21
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 85.306
## Degrees of freedom 24
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 918.852
## Degrees of freedom 36
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.931
## Tucker-Lewis Index (TLI) 0.896
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3737.745
## Loglikelihood unrestricted model (H1) -3695.092
##
## Akaike (AIC) 7517.490
## Bayesian (BIC) 7595.339
## Sample-size adjusted Bayesian (BIC) 7528.739
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.092
## 90 Percent confidence interval - lower 0.071
## 90 Percent confidence interval - upper 0.114
## P-value RMSEA <= 0.05 0.001
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.065
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## visual =~
## x1 1.000
## x2 0.554 0.100 5.554 0.000
## x3 0.729 0.109 6.685 0.000
## textual =~
## x4 1.000
## x5 1.113 0.065 17.014 0.000
## x6 0.926 0.055 16.703 0.000
## speed =~
## x7 1.000
## x8 1.180 0.165 7.152 0.000
## x9 1.082 0.151 7.155 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## visual ~~
## textual 0.408 0.074 5.552 0.000
## speed 0.262 0.056 4.660 0.000
## textual ~~
## speed 0.173 0.049 3.518 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .x1 0.549 0.114 4.833 0.000
## .x2 1.134 0.102 11.146 0.000
## .x3 0.844 0.091 9.317 0.000
## .x4 0.371 0.048 7.779 0.000
## .x5 0.446 0.058 7.642 0.000
## .x6 0.356 0.043 8.277 0.000
## .x7 0.799 0.081 9.823 0.000
## .x8 0.488 0.074 6.573 0.000
## .x9 0.566 0.071 8.003 0.000
## visual 0.809 0.145 5.564 0.000
## textual 0.979 0.112 8.737 0.000
## speed 0.384 0.086 4.451 0.000GGsem
Below sets up a function.
GGsem function
ggsem <- function(fit, layout = "sugiyama") {
# Extract standardized parameters
params <- lavaan::standardizedSolution(fit)
# Edge properties
param_edges <- params %>%
filter(op %in% c("=~", "~", "~~"), lhs != rhs, pvalue < .10) %>%
transmute(to = lhs,
from = rhs,
val = est.std,
type = dplyr::case_when(
op == "=~" ~ "loading",
op == "~" ~ "regression",
op == "~~" ~ "correlation",
TRUE ~ NA_character_))
# Identify latent variables for nodes
latent_nodes <- param_edges %>%
filter(type == "loading") %>%
distinct(to) %>%
transmute(metric = to, latent = TRUE)
# Node properties
param_nodes <- params %>%
filter(lhs == rhs) %>%
transmute(metric = lhs, e = est.std) %>%
left_join(latent_nodes) %>%
mutate(latent = if_else(is.na(latent), FALSE, latent))
# Complete Graph Object
param_graph <- tidygraph::tbl_graph(param_nodes, param_edges)
# Plot
ggraph(param_graph, layout = layout) +
# Latent factor Nodes
geom_node_point(aes(alpha = as.numeric(latent)),
shape = 16, size = 5) +
geom_node_point(aes(alpha = as.numeric(latent)),
shape = 16, size = 4, color = "white") +
# Observed Nodes
geom_node_point(aes(alpha = as.numeric(!latent)),
shape = 15, size = 5) +
geom_node_point(aes(alpha = as.numeric(!latent)),
shape = 15, size = 4, color = "white") +
# Regression Paths (and text)
geom_edge_link(aes(color = val, label = round(val, 2),
alpha = as.numeric(type == "regression")),
linetype = 1, angle_calc = "along", vjust = -.5,
arrow = arrow(20, unit(.3, "cm"), type = "closed")) +
# Factor Loadings (no text)
geom_edge_link(aes(color = val, alpha = as.numeric(type == "loading")),
linetype = 3, angle_calc = "along",
arrow = arrow(20, unit(.3, "cm"), ends = "first", type = "closed")) +
# Correlation Paths (no text)
geom_edge_link(aes(color = val, alpha = as.numeric(type == "correlation")),
linetype = 2, angle_calc = "along",
arrow = arrow(20, unit(.3, "cm"), type = "closed", ends = "both")) +
# Node names
geom_node_text(aes(label = metric),
nudge_y = .25, hjust = "inward") +
# Node residual error
geom_node_text(aes(label = sprintf("%.2f", e)),
nudge_y = -.1, size = 3) +
# Scales and themes
scale_alpha(guide = FALSE, range = c(0, 1)) +
scale_edge_alpha(guide = FALSE, range = c(0, 1)) +
scale_edge_colour_gradient2(guide = FALSE, low = "red", mid = "darkgray", high = "green") +
scale_edge_linetype(guide = FALSE) +
scale_size(guide = FALSE) +
theme_graph()
}Apply function
This is a ggplot2 object which means that one can change everything.
ggsem(fit)## Joining, by = "metric"## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
Session Info
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lavaan_0.6-9 ggraph_2.0.5 tidygraph_1.2.0 forcats_0.5.1
## [5] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.0.0
## [9] tidyr_1.1.3 tibble_3.1.4 ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] viridis_0.6.1 httr_1.4.2 sass_0.4.0 jsonlite_1.7.2
## [5] viridisLite_0.4.0 tmvnsim_1.0-2 modelr_0.1.8 bslib_0.2.5.1
## [9] assertthat_0.2.1 highr_0.9 stats4_4.1.1 cellranger_1.1.0
## [13] yaml_2.2.1 ggrepel_0.9.1 pbivnorm_0.6.0 pillar_1.6.2
## [17] backports_1.2.1 glue_1.4.2 digest_0.6.27 polyclip_1.10-0
## [21] rvest_1.0.1 colorspace_2.0-2 htmltools_0.5.2 pkgconfig_2.0.3
## [25] broom_0.7.9 haven_2.4.3 scales_1.1.1 tweenr_1.0.2
## [29] tzdb_0.1.2 ggforce_0.3.3 generics_0.1.0 farver_2.1.0
## [33] ellipsis_0.3.2 withr_2.4.2 cli_3.0.1 mnormt_2.0.2
## [37] magrittr_2.0.1 crayon_1.4.1 readxl_1.3.1 evaluate_0.14
## [41] fs_1.5.0 fansi_0.5.0 MASS_7.3-54 xml2_1.3.2
## [45] tools_4.1.1 hms_1.1.0 lifecycle_1.0.0 munsell_0.5.0
## [49] reprex_2.0.1 compiler_4.1.1 jquerylib_0.1.4 rlang_0.4.11
## [53] grid_4.1.1 rstudioapi_0.13 igraph_1.2.6 labeling_0.4.2
## [57] rmarkdown_2.10 gtable_0.3.0 DBI_1.1.1 graphlayouts_0.7.1
## [61] R6_2.5.1 gridExtra_2.3 lubridate_1.7.10 knitr_1.33
## [65] fastmap_1.1.0 utf8_1.2.2 stringi_1.7.4 Rcpp_1.0.7
## [69] vctrs_0.3.8 dbplyr_2.1.1 tidyselect_1.1.1 xfun_0.25