其他
独家 R 语言小工具：结构方程模型从拟合到作图，完全使用 R 解决！
引言
完全在 R 语言中，拟合，并使用辅助小工具（一个简单的 shiny APP，见视频演示）绘制结构方程模型图示。
这个 APP 小工具可以保存结构方程模型图示为 PDF 或 PPT 格式（以便后续进一步编辑）。
帮助本公众号传播的（包括转发到票圈、群聊等），或打赏的同学，请私信我截图和邮箱，我可以将此 APP 的代码发送到您的邮箱。
1、工具
library(tidyverse)
library(lavaan)
library(igraph)
library(ggraph)
library(tidygraph)
library(eoffice)2、数据
来自这篇论文的公开数据，其代码也是公开且完全可运行的，
Land management shapes drought responses of dominant soil microbial taxa across grasslands注意，dat_sem_t0 数据较长，此处无法展示。不过 dat_sem_t0 可以运行该论文的代码后获取，
dat_sem_t0 <- read_csv("dat_sem_t0.csv")dat_sem_t0 即论文代码中的 dat_sem |> filter(day == "0")，
dat_sem_t0 <- dat_sem |> filter(day == "0")3、SEM 的结构
"## soil properties and pH separately ----
soil_var ~ lat + management
pH ~ lat + management
## soil moisture ----
soil_moisture ~ treatment + management + lat + soil_var
## bacterial groups ----
index_bact_opp ~ treatment + management + pH + lat + soil_moisture + soil_var
index_bact_resis ~ treatment + management + pH + lat + soil_moisture + soil_var
index_bact_sens ~ treatment + management + pH + lat + soil_moisture + soil_var
## extrinsic covariances ----
management ~~ 0 * treatment
pH ~~ soil_var
pH ~~ soil_moisture" -> model_sem_B4、拟合 SEM
fit_model_sem_B0 <- sem(
model_sem_B,
data = dat_sem_t0,
std.lv = TRUE,
std.ov = TRUE,
fixed.x = FALSE
)拟合效果，
summary(
fit_model_sem_B0,
fit.measures = TRUE,
standardized = T,
rsq = T
)## lavaan 0.6.17 ended normally after 25 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 40
##
## Number of observations 60
##
## Model Test User Model:
##
## Test statistic 2.531
## Degrees of freedom 5
## P-value (Chi-square) 0.772
##
## Model Test Baseline Model:
##
## Test statistic 314.687
## Degrees of freedom 36
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.064
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -523.442
## Loglikelihood unrestricted model (H1) NA
##
## Akaike (AIC) 1126.883
## Bayesian (BIC) 1210.657
## Sample-size adjusted Bayesian (SABIC) 1084.846
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.121
## P-value H_0: RMSEA <= 0.050 0.823
## P-value H_0: RMSEA >= 0.080 0.119
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.040
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## soil_var ~
## lat 0.569 0.106 5.361 0.000 0.569 0.569
## management -0.071 0.210 -0.340 0.734 -0.071 -0.036
## pH ~
## lat -0.478 0.100 -4.772 0.000 -0.478 -0.478
## management 0.816 0.199 4.107 0.000 0.816 0.411
## soil_moisture ~
## treatment -1.392 0.158 -8.816 0.000 -1.392 -0.722
## management 0.226 0.159 1.427 0.154 0.226 0.117
## lat 0.244 0.097 2.509 0.012 0.244 0.251
## soil_var -0.255 0.097 -2.626 0.009 -0.255 -0.263
## index_bact_opp ~
## treatment 0.724 0.224 3.228 0.001 0.724 0.369
## management 0.428 0.169 2.533 0.011 0.428 0.218
## pH 0.548 0.097 5.671 0.000 0.548 0.553
## lat 0.274 0.106 2.581 0.010 0.274 0.277
## soil_moisture -0.237 0.121 -1.957 0.050 -0.237 -0.233
## soil_var -0.033 0.096 -0.345 0.730 -0.033 -0.033
## index_bact_resis ~
## treatment -0.010 0.301 -0.032 0.975 -0.010 -0.005
## management 0.439 0.227 1.938 0.053 0.439 0.221
## pH 0.483 0.130 3.725 0.000 0.483 0.482
## lat 0.057 0.143 0.396 0.692 0.057 0.056
## soil_moisture -0.100 0.162 -0.617 0.538 -0.100 -0.097
## soil_var -0.146 0.129 -1.133 0.257 -0.146 -0.146
## index_bact_sens ~
## treatment -1.360 0.274 -4.954 0.000 -1.360 -0.698
## management -0.066 0.207 -0.321 0.748 -0.066 -0.034
## pH 0.060 0.118 0.505 0.614 0.060 0.061
## lat 0.243 0.130 1.866 0.062 0.243 0.247
## soil_moisture -0.042 0.148 -0.283 0.777 -0.042 -0.041
## soil_var -0.181 0.117 -1.542 0.123 -0.181 -0.184
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## management ~~
## treatment 0.000 0.000 0.000
## .soil_var ~~
## .pH -0.021 0.081 -0.255 0.798 -0.021 -0.033
## .pH ~~
## .soil_moisture 0.037 0.061 0.603 0.547 0.037 0.078
## .index_bact_opp ~~
## .index_bact_rss 0.301 0.069 4.368 0.000 0.301 0.683
## .index_bact_sns 0.097 0.053 1.811 0.070 0.097 0.241
## .index_bact_resis ~~
## .index_bact_sns 0.297 0.080 3.738 0.000 0.297 0.551
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .soil_var 0.664 0.121 5.477 0.000 0.664 0.675
## .pH 0.592 0.108 5.477 0.000 0.592 0.602
## .soil_moisture 0.376 0.069 5.477 0.000 0.376 0.405
## .index_bact_opp 0.329 0.060 5.477 0.000 0.329 0.341
## .index_bact_rss 0.592 0.108 5.477 0.000 0.592 0.602
## .index_bact_sns 0.492 0.090 5.477 0.000 0.492 0.519
## management 0.250 0.046 5.477 0.000 0.250 1.000
## treatment 0.250 0.046 5.477 0.000 0.250 1.000
## lat 0.983 0.180 5.477 0.000 0.983 1.000
##
## R-Square:
## Estimate
## soil_var 0.325
## pH 0.398
## soil_moisture 0.595
## index_bact_opp 0.659
## index_bact_rss 0.398
## index_bact_sns 0.4815、SEM 重抽样
耗时取决于不同电脑配置。应该不到一分钟，
system.time(
B0_boot <- bootstrapLavaan(
fit_model_sem_B0,
R = 1000L,
type = "bollen.stine",
verbose = FALSE,
FUN = fitMeasures,
fit.measures = c("chisq")
)
)## user system elapsed
## 36.31 0.89 36.92重抽样 P 值，
# Compute a bootstrap based p-value
B0_fit_orig <- fitMeasures(fit_model_sem_B0, "chisq")
B0_boot_pvalue <- length(which(B0_boot > B0_fit_orig)) / length(B0_boot); B0_boot_pvalue## [1] 0.786、获取 SEM 结果，转换为 graph
小技巧，relocate 很好用，可以调整列的位置，同时重命名该列
sem_graph <- standardizedSolution(fit_model_sem_B0) |>
filter(op == "~") |>
relocate(from = rhs, to = lhs) |>
as_tbl_graph(directed = TRUE)绘制 SEM，
sem_graph |>
ggraph(layout = "circle") +
geom_edge_link(
aes(
edge_width = abs(est.std),
edge_color = est.std > 0,
edge_alpha = pvalue < 0.05
),
start_cap = circle(10, "pt"),
end_cap = circle(10, "pt"),
arrow = arrow(45, unit(5, "pt"), type = "closed")
) +
scale_edge_width_continuous(range = c(0.2, 2)) +
# geom_node_point(shape = 1, size = 9) +
geom_node_text(aes(label = name)) +
coord_equal(clip = "off") +
theme_void()## Warning: Using alpha for a discrete variable is not advised.## Warning: Using the `size` aesthetic in this geom was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` in the `default_aes` field and elsewhere instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.7、进阶！调整 SEM 图示
首先，保存 SEM 结果表格，将这个表格输入到以下 APP 中，即可开始调整，
standardizedSolution(fit_model_sem_B0) |>
filter(op == "~") |>
relocate(from = rhs, to = lhs, weight = est.std, p = pvalue) |>
write_csv("sem_table.csv")然后，使用 shiny APP 编辑 SEM 图示，
见开头的视频演示
帮助本公众号传播的同学——包括但不限于转发到票圈、群聊等，请私信我一个有效的截图及邮箱，我可以将此应用的代码发送到您的邮箱
对于某些较为简单的 SEM，此 APP 还是可用的
source("shiny_app_sem_graph_advanced.R")
shinyApp(ui = ui, server = server)调整完成后，保存得到一个 PPT（或 PDF）文件，
同时，也会保存最终布局为
sem_layout.csv以便重复使用图中，默认忽略了不显著路径
若要显示不显著路径，注释掉
shiny_app_sem_graph_advanced.R里的第 187 行代码，然后再次运行，重新调整布局