引言
经常看到,非常复杂的、节点十几个、路径数十条的结构方程模型(SEM)。怎么更简明地可视化这种 SEM 的路径图?
有办法!之前的推送介绍过,所谓 Manu 风格的 SEM 路径图
。例如以下一篇最近的 PNAS 论文中所示,
❗️ Positive associations fuel soil biodiversity and ecological
networks worldwide
👉 近期关于结构方程模型的 10 篇推送|Manu 风格的 SEM 路径图
(点击 👉 即可跳转。)
那么,怎么才能方便、快速地可视化这种 SEM?
我最近研究了一个 R 语言 sniny APP
小工具,可以在得到 SEM 的结果以后(需要的数据格式见后文),快速可视化 SEM 路径图。
该 APP 的主要功能特点包括(详见以下视频教程),
1、
上传 csv 数据,快速可视化 SEM 路径图
;2、
鼠标点击,手动移动 SEM 路径图节点的位置
;3、
调整完成后,自动生成 R 代码,以确保可视化结果可复现
。
演示视频如下(展示了 APP 在 RStudio 中运行的效果),
1、
打赏 16 元以上
;2、
或者,转发到票圈、集赞 10 个以上、并保持可见 12 小时,私信截图
;3、
最后,记得留下您的邮箱
,即可获得。4、本文后附该
APP 生成的 R 代码
,效果可供参考。
此外,我在之前也做过多个 SEM 的不写代码、拟合、可视化工具(lavaan、piecewiseSEM),另见以下几篇推送,
👉 复现 Science 论文 piecewiseSEM 结构方程模型|但不用写 R 代码
👉 同时做 20 个 lavaan SEM 的模型选择|但不用写 R 代码
👉 结构方程模型的可视化|Nature 子刊 SEM
👉 使用鼠标交互式可视化 SEM 路径图|独家 shiny APP 工具
👉 近期关于结构方程模型的 10 篇推送|Manu 风格的 SEM 路径图
(点击 👉 即可跳转。)
代码
以下,是该 APP 生成的 R 代码
(稍加微调)。
工具包
library(igraph)
library(ggraph)
library(tidygraph)
library(dplyr)
library(svglite)
数据
得到的 SEM 结果存在 csv 文件中,数据格式如下,
另存为
(manu) sem-table.csv
。
read.csv('(manu) sem-table.csv') |>
dplyr::select(from, to, weight, p)
from to weight p
1 fragility biodiversity 0.8182 0.0000
2 cycfac_total biodiversity 0.6357 0.0000
3 modularity biodiversity -0.5062 0.0000
4 Distance_from_equator biodiversity 0.0123 0.6906
5 Elevation biodiversity 0.1824 0.0000
6 TSEA biodiversity -0.1648 0.0000
7 MAP biodiversity 0.1431 0.0000
8 pH biodiversity -0.0345 0.2782
9 Texture biodiversity -0.0135 0.6913
10 Soil_C biodiversity -0.0807 0.0069
11 cycfac_total fragility -0.2172 0.0000
12 modularity fragility 0.3899 0.0000
13 Distance_from_equator fragility 0.0096 0.7731
14 Elevation fragility 0.0601 0.0974
15 TSEA fragility -0.1916 0.0000
16 MAP fragility -0.0793 0.0236
17 pH fragility -0.0734 0.0309
18 Texture fragility 0.3382 0.0000
19 Soil_C fragility 0.0446 0.1619
20 cycfac_total modularity -0.2402 0.0000
21 Distance_from_equator modularity 0.3103 0.0000
22 Elevation modularity 0.1090 0.0116
23 TSEA modularity -0.2543 0.0000
24 MAP modularity -0.2919 0.0000
25 pH modularity 0.0182 0.6544
26 Texture modularity 0.1997 0.0000
27 Soil_C modularity -0.0998 0.0087
28 Distance_from_equator cycfac_total -0.0617 0.1130
29 Elevation cycfac_total -0.0933 0.0365
30 TSEA cycfac_total 0.5755 0.0000
31 MAP cycfac_total 0.2301 0.0000
32 pH cycfac_total 0.0370 0.3796
33 Texture cycfac_total -0.0466 0.2553
34 Soil_C cycfac_total -0.0623 0.1130
35 Distance_from_equator pH 0.0171 0.6744
36 Elevation pH 0.1621 0.0001
37 TSEA pH 0.0641 0.0993
38 MAP pH -0.4652 0.0000
39 Distance_from_equator Soil_C 0.0308 0.4521
40 Elevation Soil_C 0.1538 0.0010
41 TSEA Soil_C 0.2732 0.0000
42 MAP Soil_C 0.0927 0.0293
43 pH Soil_C -0.3633 0.0000
44 Texture Soil_C 0.2147 0.0000
45 Distance_from_equator Texture -0.0398 0.3217
46 Elevation Texture 0.4339 0.0000
47 TSEA Texture 0.3736 0.0000
48 MAP Texture -0.1233 0.0008
49 Distance_from_equator TSEA 0.1617 0.0002
50 Elevation TSEA -0.2323 0.0000
51 Distance_from_equator MAP -0.0750 0.0983
52 Elevation MAP -0.2043 0.0000
53 TSEA MAP -0.2622 0.0000
选择节点
# Function to merge selected nodes
mergeSelectedNodes <- function(data, nodeGroups) {
for (nodes in nodeGroups) {
mergedNode <- paste(nodes, collapse = '_')
data <- data |>
mutate(
from = ifelse(from %in% nodes, mergedNode, from),
to = ifelse(to %in% nodes, mergedNode, to)
)
}
return(data)
}
# Node groups
nodeGroups <- list(
c(
"Distance_from_equator",
"Elevation",
"TSEA",
"MAP",
"pH",
"Texture",
"Soil_C"
)
)
# Make graph
graph <- read.csv('(manu) sem-table.csv') |>
mutate(
FROM_TO_WEIGHT_P = paste0(from, '→', to, ' (', round(weight, 3), ', ', round(p, 3), ')'),
WEIGHT_FROM_TO = paste0(round(weight, 3), gtools::stars.pval(p), ' (', from, '→', to, ')'),
WEIGHT_FROM = paste0(round(weight, 3), gtools::stars.pval(p), ' (', from, '→)'),
sign = ifelse(weight > 0, '+', '-'),
weight = abs(weight)
) |>
mergeSelectedNodes(nodeGroups) |>
graph_from_data_frame()
指定新布局坐标
# Define the layout coordinates
layoutCoordinates <- data.frame(
name = c(
"fragility",
"cycfac_total",
"modularity",
"Distance_from_equator_Elevation_TSEA_MAP_pH_Texture_Soil_C",
"biodiversity"
),
x = c(19, 12, 1, 1, 12),
y = c(10, 19, 16, 4, 1)
)
# Create the layout manually
layout <- create_layout(
graph,
layout = 'manual',
x = layoutCoordinates$x,
y = layoutCoordinates$y
)
新布局可视化
风格 0
This figure cannot be used directly, but serve as the basis for the following figures.
p1 <- ggraph(layout) +
geom_edge_parallel(
sep = unit(8, "pt"),
aes(
color = sign,
edge_width = abs(weight),
edge_linetype = ifelse(p < 0.05, "SIG", "NS"),
label = WEIGHT_FROM
),
angle_calc = "along",
label_size = 0,
label_push = unit(10, "pt"),
label_dodge = unit(10, "pt"),
arrow = arrow(60, length = unit(7, "pt"), type = "closed"),
start_cap = circle(0.5),
end_cap = circle(0.5)
) +
geom_node_label(
aes(label = name),
size = 3,
alpha = 0.5,
) +
scale_color_brewer(
name = NULL,
palette = "Dark2"
) +
scale_edge_color_brewer(
name = NULL,
palette = "Set2",
direction = -1
) +
scale_edge_width_continuous(
guide = "none",
range = c(0.2, 2)
) +
scale_edge_linetype_manual(
name = NULL,
values = c("SIG" = 1, "NS" = 2)
) +
coord_cartesian(
xlim = c(0, 20),
ylim = c(0, 20),
clip = "off"
) +
theme_graph()
# Show plot
p1
之所以保存为
SVG
格式,是因为SVG 能在 PPT 里编辑
,方便进一步美化。
# Now you can edit the SVG in PPT
ggsave('Graph1.svg', p1, width = 10, height = 8)
风格 1
# Edge data
biuldData <- ggplot_build(p1)
nodeNames <- igraph::V(graph)$name
edgeData <- biuldData$data[[1]] |>
mutate(
from0 = nodeNames[from],
to0 = nodeNames[to]
)
# Summarize edge data
edgeDataSumm <- edgeData |>
summarise(
x = mean(x),
y = mean(y),
.by = c(label, from0, to0)
) |>
mutate(
label.vjust = seq(-dplyr::n(), dplyr::n(), len = dplyr::n()) / 2,
.by = c(from0, to0)
)
# Plot
p2 <- layout |>
ggraph() +
geom_edge_parallel(
sep = unit(8, 'pt'),
aes(
color = sign,
edge_width = weight,
edge_linetype = ifelse(p < 0.05, 'SIG', 'NS')
),
start_cap = circle(0.5),
end_cap = circle(0.5),
arrow = arrow(60, length = unit(7, 'pt'), type = 'closed')
) +
geom_label(
data = edgeDataSumm,
aes(
x,
y,
label = label,
color = to0,
vjust = label.vjust,
),
size = 3,
hjust = 0,
inherit.aes = FALSE
) +
geom_node_label(
aes(label = name),
size = 3,
label.r = unit(0, 'pt'),
alpha = 0.5,
) +
scale_color_brewer(
name = NULL,
palette = "Dark2"
) +
scale_edge_color_brewer(
name = NULL,
palette = "Set2",
direction = -1
) +
scale_edge_width_continuous(
guide = "none",
range = c(0.2, 2)
) +
scale_edge_linetype_manual(
name = NULL,
values = c("SIG" = 1, "NS" = 2)
) +
coord_cartesian(
xlim = c(0, 20),
ylim = c(0, 20),
clip = "off"
) +
theme_graph()
# Show plot
p2
之所以保存为
SVG
格式,是因为SVG 能在 PPT 里编辑
,方便进一步美化。
# Now you can edit the SVG in PPT
ggsave('Graph2.svg', p2, width = 10, height = 8)
风格 2
# Slice, select required columns, and calculate slope and angle in one step
edgeDataJoin <- edgeData |>
slice(1:2, .by = label) |>
mutate(
slope = (y - lead(y)) / (x - lead(x)),
angle = atan(slope) * (180 / pi)
) |>
slice(1, .by = c(label, from0)) |>
select(label, angle, slope) |>
full_join(edgeDataSumm, by = 'label') |>
mutate(
label.vjust = ifelse(
slope > 0,
seq(dplyr::n(), -dplyr::n(), len = dplyr::n()) / 2 + 0.5,
seq(-dplyr::n(), dplyr::n(), len = dplyr::n()) / 2 + 0.5
),
.by = c(from0, to0)
) |>
select(label, from0, to0, x, y, label.vjust, angle)
# Plot
p3 <- layout |>
ggraph() +
geom_edge_parallel(
sep = unit(8, 'pt'),
aes(
edge_width = weight,
edge_linetype = ifelse(p < 0.05, 'SIG', 'NS')
),
start_cap = circle(0.5),
end_cap = circle(0.5),
arrow = arrow(60, length = unit(7, 'pt'), type = 'closed')
) +
geom_label(
data = edgeDataJoin,
aes(
x,
y,
label = label,
color = from0,
vjust = label.vjust,
angle = angle
),
size = 2.5,
hjust = 0.5,
inherit.aes = FALSE
) +
geom_node_label(aes(label = name), size = 3.5) +
scale_color_brewer(
name = NULL,
palette = "Dark2"
) +
scale_edge_color_brewer(
name = NULL,
palette = "Set2",
direction = -1
) +
scale_edge_width_continuous(
guide = "none",
range = c(0.2, 2)
) +
scale_edge_linetype_manual(
name = NULL,
values = c("SIG" = 1, "NS" = 2)
) +
coord_cartesian(
xlim = c(0, 20),
ylim = c(0, 20),
clip = "off"
) +
theme_graph() +
theme(legend.position = 'top')
# Show plot
p3
之所以保存为
SVG
格式,是因为SVG 能在 PPT 里编辑
,方便进一步美化。
# Now you can edit the SVG in PPT
ggsave('Graph3.svg', p3, width = 10, height = 8)