Stat009
Stat009

Reputation: 45

Path diagram in r

I am trying to plot a path diagram of a Structural Equation Model(SEM) in R. I was able to plot it using semPlot::semPaths(). The output is similar to enter image description here The SEM was modeled using lavaan package.

I want a plot similar to enter image description here. with estimates and p values. Can anyone help me out?

Upvotes: 3

Views: 9314

Answers (1)

Sinval
Sinval

Reputation: 1417

My suggestion would be lavaanPlot (see more of it in the author's personal website):

library(lavaan)
library(lavaanPlot)
# path model
model <- 'mpg ~ cyl + disp + hp
          qsec ~ disp + hp + wt'

fit1 <- sem(model, data = mtcars)
labels1 <- list(mpg = "Miles Per Gallon", cyl = "Cylinders", disp = "Displacement", hp = "Horsepower", qsec = "Speed", wt = "Weight") #define labels
lavaanPlot(model = fit1, labels = labels1, coefs = TRUE, stand = TRUE, sig = 0.05) #standardized regression paths, showing only paths with p<= .05

enter image description here

Edit

After trying several options, lavaanPlot included, I believe that the best existing solution for publish-ready plots is the combination of semPlot with semptools. See this example:

library(lavaan)
library(semPlot)
#remotes::install_github("sfcheung/semptools")
library(semptools)

model <- "
#regression paths
x2 ~ x1
x3 ~ x1
x4 ~ x2 + x3 +x5
x5 ~ x3
x6 ~ x3 + x5
x7 ~ x5
x8 ~ x4 + x5 + x6 + x7

#covariance paths
x2 ~~ x3
"

fit_model <- sem(model, data = HolzingerSwineford1939)

layout <- matrix(c(NA,   "x2", "x4", NA,  NA,
                          "x1", NA,   NA,   NA,  "x8",
                          NA,   "x3", "x5", "x7", NA,
                          "SA", NA,   "x6", NA,   NA),
                        byrow = TRUE,
                        ncol = 5,
                        nrow = 4)



plot_model <- semPlotModel_lavaanModel(model,auto.var = TRUE,
                                       auto.fix.first = F,
                                       auto.cov.lv.x = TRUE)

diagram_plot <- semPlot::semPaths(object = plot_model,
                                  thresholds = F,
                                  whatLabels = "none",
                                  layout = layout,
                                  intercepts = F,
                                  exoCov = T,
                                  style = "ram",
                                  curvePivot=F,
                                  edge.color="black",
                                  DoNotPlot = T,
                                  residuals = F,
                                  edgeLabels = c("a",
                                                 "d",
                                                 "b",
                                                 "e",
                                                 "f",
                                                 "g",
                                                 "h",
                                                 "i",
                                                 "j",
                                                 "c",
                                                 "m",
                                                 "k",
                                                 "l",
                                                 ""),
                                  sizeMan = 8)


my_curve_list <- list(list(from = "x2",
                           to = "x3",
                           new_curve = -1),
                      list(from = "x6",
                           to = "x8",
                           new_curve = -1.5))

diagram_plot <- set_curve(diagram_plot, my_curve_list)

vars_names <- c("Sanitation coverage",
                "Piped water coverage",
                "Household sanitation use",
                "Water availability",
                "Handwashing station coverage",
                "Drinking water storage",
                "HAZ",
                "Intervention")


model_plot_names <- gsub(" ", "\n", vars_names)

diagram_plot_long_names <- change_node_label(semPaths_plot = diagram_plot,
                                            label_list = setNames(object = model_plot_names,
                                                                   nm = diagram_plot$Arguments$labels),
                                             label.cex = .60, label.scale=F)

plot(diagram_plot_long_names)

enter image description here

Now, the plot with estimated values...

plot_model <- semPlotModel(fit_model)


plot_path_analysis <- semPlot::semPaths(object = plot_model,
                              thresholds = F,
                              whatLabels = "std",
                              layout = latent_layout,
                              intercepts = F,
                              exoCov = T,
                              style = "ram",
                              curvePivot=F,
                              edge.color="black",
                              DoNotPlot = T,
                              sizeMan = 7)

plot_path_analysis <- set_curve(plot_path_analysis, my_curve_list)
plot(plot_path_analysis)

my_rotate_resid_list <- c(x3 = 270,
                          x5 = 315,
                          x6 = 180,
                          x7 = 180)

plot_path_analysis <- rotate_resid(plot_path_analysis, my_rotate_resid_list)

plot_path_analysis <- mark_sig(plot_path_analysis,
                       fit_model)

plot_path_analysis <- change_node_label(plot_path_analysis,
                         setNames(model_plot_names, 
                                  plot_path_analysis$Arguments$labels),
                         label.cex = .60,
                         label.scale=F)

plot(plot_path_analysis)

enter image description here

Upvotes: 3

Related Questions