引言

又有一个 SEM 代码学习好案例！最新出炉的好文章，公开了代码：

• 代码地址见原文 Code availability 一节。

Larger nations benefit more than smaller nations from the stabilizing effects of crop diversity

一、加载工具包

library(nlme)
library(piecewiseSEM)
library(igraph)
library(ggraph)
library(tidygraph)
library(dplyr)

二、原文所附 SEM 代码

# 加载数据
load("NATFOOD-23060800_data.RData")

# 拟合 SEM 的代码
figure1b_SEM <- psem(
  # 回归：
  gls(
    log10(pre.cv) ~ log10(Area.h),
    correlation = corAR1(form = ~ 1 | Country),
    Nation_CV
  ),
  gls(
    log10(LPre.n) ~ log10(Area.h),
    correlation = corAR1(form = ~ 1 | Country),
    Nation_CV
  ),
  gls(
    log10(asy) ~ EH + log10(tmp.cv),
    correlation = corAR1(form = ~ 1 | Country),
    Nation_CV
  ),
  gls(
    log10(Avg.crop.S) ~ log10(Area.h) + log10(pre.cv) + log10(tmp.cv) + log10(LPre.n) + log10(Irrigation),
    correlation = corAR1(form = ~ 1 | Country),
    Nation_CV
  ),
  lm(
    log10(Yield.S) ~ log10(Avg.crop.S) + log10(asy),
    Nation_CV
  ),
  
  # 相关：
  log10(Avg.crop.S) %~~% log10(asy),
  log10(tmp.cv) %~~% log10(pre.cv),
  log10(Irrigation) %~~% log10(pre.cv),
  log10(LPre.n) %~~% log10(pre.cv),
  EH %~~% log10(pre.cv)
)

三、保存 SEM 数据

• 以 CSV 格式保存。

# 结果整理为 from、to、weight、p 数据框
coefs(figure1b_SEM) |>
  select(-9) |> # 去掉星号列
  filter(
    str_detect(Response, "~~", negate = TRUE) # 去掉相关性部分
  ) |> 
  dplyr::relocate(
    from = Predictor,
    to = Response,
    weight = Std.Estimate,
    p = P.Value
  ) |> 
  write.csv("nature-food-sem-table.csv", row.names = FALSE)

四、SEM 结果变成图

• SEM 是一个有向无环图：DAG（Directed Acyclic Graph）。

# 变为一个 graph
graph <- read.csv('nature-food-sem-table.csv') |>
  graph_from_data_frame()
print(graph)

IGRAPH d7deb86 DNW- 9 11 -- 
+ attr: name (v/c), weight (e/n), p (e/n), Estimate (e/n), Std.Error
| (e/n), DF (e/n), Crit.Value (e/n)
+ edges from d7deb86 (vertex names):
 [1] log10(Area.h)    ->log10(pre.cv)     log10(Area.h)    ->log10(LPre.n)    
 [3] EH               ->log10(asy)        log10(tmp.cv)    ->log10(asy)       
 [5] log10(Area.h)    ->log10(Avg.crop.S) log10(pre.cv)    ->log10(Avg.crop.S)
 [7] log10(tmp.cv)    ->log10(Avg.crop.S) log10(LPre.n)    ->log10(Avg.crop.S)
 [9] log10(Irrigation)->log10(Avg.crop.S) log10(Avg.crop.S)->log10(Yield.S)   
[11] log10(asy)       ->log10(Yield.S)

# 是一个 DAG 吗？
igraph::is_dag(graph)

[1] TRUE

五、Style 1：设置默认布局

绘图

font <- 'DejaVu Sans'

p0 <- ggraph(graph, layout = "sugiyama") +
  geom_edge_diagonal(
    strength = 1,
    aes(
      color = ifelse(weight > 0, "+", "-"),
      edge_width = abs(weight),
      edge_linetype = ifelse(p < 0.05, 'SIG', 'NS'),
      label = round(weight, 2),
      start_cap = label_rect(node1.name),
      end_cap = label_rect(node2.name)
    ),
    angle_calc = 'along',
    label_size = 3,
    label_push = unit(10, 'pt'),
    label_dodge = unit(10, 'pt'),
    family = font,
    arrow = arrow(60, length = unit(7, 'pt'), type = 'closed'),
    linejoin = "mitre",
    linemitre = 2
  ) +
  geom_node_label(
    aes(label = name),
    size = 3.5,
    label.r = unit(0, "pt"),
    alpha = 0.5,
    family = font
  ) +
  scale_color_brewer(name = NULL, palette = 'Paired') +
  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(clip = 'off') +
  theme_graph(base_family = 'DejaVu Sans')

# Show plot
p0

存图为 SVG 格式，可以进一步在 PPT 里编辑

• 修改变量名、添加 R2 等等。

# Now you can edit the SVG in PPT
ggsave('Graph0.svg', p0, width = 10, height = 5)

六、Style 2：手动指定布局（节点位置）

• 以下手动布局的坐标，来自一个 SEM 可视化工具，详情见另几篇推送：

👉 R 语言 SEM｜鼠标修改布局｜结构方程模型可视化工具教程｜全球独家
👉 R 语言 SEM｜结构方程模型直接间接效应｜可视化

（点击 👉 即可跳转。）

• 空间想象能力好的话，可以自定义坐标。

# Define the layout coordinates
layout_coord <- data.frame(
  name = c('log10(Area.h)', 'EH', 'log10(tmp.cv)', 'log10(pre.cv)', 'log10(LPre.n)', 'log10(Irrigation)', 'log10(Avg.crop.S)', 'log10(asy)', 'log10(Yield.S)'),
  x = c(5, 19, 12, 1, 5, 15, 8, 19, 14),
  y = c(19, 19, 19, 13, 13, 19, 7, 7, 1)
)

# Create the layout manually
layout <- create_layout(
  graph,
  layout = 'manual',
  x = layout_coord$x,
  y = layout_coord$y
)

绘图

p1 <- ggraph(layout) +
  geom_edge_diagonal(
    strength = 1,
    aes(
      color = ifelse(weight > 0, "+", "-"),
      edge_width = abs(weight),
      edge_linetype = ifelse(p < 0.05, 'SIG', 'NS'),
      label = round(weight, 2),
      start_cap = label_rect(node1.name),
      end_cap = label_rect(node2.name)
    ),
    angle_calc = 'along',
    label_size = 3,
    label_push = unit(10, 'pt'),
    label_dodge = unit(10, 'pt'),
    family = font,
    arrow = arrow(60, length = unit(7, 'pt'), type = 'closed'),
    linejoin = "mitre",
    linemitre = 2
  ) +
  geom_node_label(
    aes(label = name),
    size = 3.5,
    label.r = unit(0, "pt"),
    alpha = 0.5,
    family = font
  ) +
  scale_color_brewer(name = NULL, palette = 'Paired') +
  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(clip = 'off') +
  theme_graph(base_family = 'DejaVu Sans')

# Show plot
p1

存图为 SVG 格式，可以进一步在 PPT 里编辑

• 修改变量名、添加 R2 等等。

# Now you can edit the SVG in PPT
ggsave('Graph1.svg', p1, width = 10, height = 5)