R 语言 SEM|结构方程模型直接间接效应|可视化
![](/skin/images/mar_09.png)
![](/skin/images/mar_10.png)
![](/skin/images/mar_07.png)
![](/skin/images/mar_08.png)
📢 预告:新版 R 语言 SEM 可视化 Shiny APP 将于近期推送。
往期
👉 鼠标修改布局!全球独家 R 语言 SEM 可视化工具👉 完全用 R 语言实现 Manuel Delgado-Baquerizo 风格的 SEM 图👉️ 尝试用 R 语言画一个 PNAS 文章的 SEM 图👉 R 语言 SEM 之多组分析(Multigroup)👉 R 语言 SEM 之复合变量(Composite)的构造案例👉 R 语言 SEM 之使用 ggplot2 衍生包画 Nature 结构方程模型图👉 R 语言 SEM 之贝叶斯(Bayesian)结构方程模型👉 R 语言 SEM 合集(点击 👉 即可跳转。)目录
引言1 需要的包2 需要的数据3 使用 piecewiseSEM 拟合4 SEM 可视化:计算直接和间接效应 4.0 用鼠标绘制 SEM 图示 APP(更新) 4.1 第一步:将 SEM 结果变成网络 4.2 绘图:结构方程模型路径图 4.3 绘图:直接和间接效应6 SEM:重抽样的直接和间接效应 6.1 重抽样:直接和间接效应 6.2 重抽样:直接和间接效应的置信区间7 附:getEffects 函数引言
探索使用代码(而不是动手)计算 SEM 的直接和间接效应。
1 需要的包
library(tidyverse)library(piecewiseSEM)library(tidygraph)library(ggraph)library(igraph)2 需要的数据
数据标准化,
# 来自 piecewiseSEM 包data(keeley)# 数据标准化keeleyScaled <- keeley |> dplyr::mutate_if( is.numeric, ~ as.numeric(scale(.x)) )str(keeleyScaled)data.frame: 90 obs. of 8 variables:
$ distance: num 0.473 -1.381 0.505 0.505 0.309 ...
$ elev : num 3.1 -1.41 -0.87 -0.87 2.11 ...
$ abiotic : num 1.489 -1.08 0.228 1.552 -0.335 ...
$ age : num 1.1486 -0.0451 -0.8409 -0.8409 -0.2043 ...
$ hetero : num 0.6423 -1.672 1.4037 0.0656 -1.1992 ...
$ firesev : num -0.645 -0.312 -1.189 -1.008 -0.16 ...
$ cover : num 1.096 -0.673 0.812 1.588 1.913 ...
$ rich : num 0.117 -1.207 1.441 0.978 1.242 ...3 使用 piecewiseSEM 拟合
# 拟合 SEMpSem <- psem( lm(age ~ distance, data = keeleyScaled), lm(hetero ~ distance, data = keeleyScaled), lm(abiotic ~ distance, data = keeleyScaled), lm(firesev ~ age, data = keeleyScaled), lm(cover ~ firesev, data = keeleyScaled), lm(rich ~ distance + hetero + abiotic + cover, data = keeleyScaled))# 路径系数# coefs(pSem) |> print()4 SEM 可视化:计算直接和间接效应
整理 SEM 结果,
# 结果整理为 from、to、weight、p 数据框semCoefs <- coefs(pSem) |> dplyr::relocate( from = Predictor, to = Response, weight = Estimate, p = P.Value )semCoefs |> print(row.names = FALSE)from to weight p Std.Error DF Crit.Value Std.Estimate
distance age -0.2781 0.0079 0.1024 88 -2.7164 -0.2781 **
distance hetero 0.3460 0.0008 0.1000 88 3.4593 0.3460 ***
distance abiotic 0.4597 0.0000 0.0947 88 4.8562 0.4597 ***
age firesev 0.4539 0.0000 0.0950 88 4.7781 0.4539 ***
firesev cover -0.4371 0.0000 0.0959 88 -4.5594 -0.4371 ***
distance rich 0.2829 0.0021 0.0892 85 3.1710 0.2829 **
hetero rich 0.3383 0.0001 0.0825 85 4.1035 0.3383 ***
abiotic rich 0.2495 0.0037 0.0837 85 2.9815 0.2495 **
cover rich 0.2866 0.0005 0.0790 85 3.6297 0.2866 ***对于 lavaan、plspm 或者 brms 的结构方程模型的整理,见以下超链接,
👉 拟合 SEM:lavaan、piecewiseSEM 和 plspm 三个包(点击 👉 即可跳转。)4.0 用鼠标绘制 SEM 图示 APP(更新)
将结果整理为 from、to、weight、p 数据框,并存为 .csv 文件,还可以使用独家 APP,用鼠标点击,修改 SEM 图示的布局!
APP 有一个小更新,见本期另一篇推送。
# 保存数据write.csv(semCoefs, "sem_table.csv", row.names = FALSE)打赏并私信,可以获取此 APP!
详情见本期另一篇推送,
另见以下超链接,
👉 全球独家 SEM 工具:用鼠标点击,可视化结构方程模型(点击 👉 即可跳转。)4.1 第一步:将 SEM 结果变成网络
此处,对于 piecewiseSEM 的结果,可以这样转换为一个网络,
# igraph 和 tidygraph 包,将数据转换为一个图(graph)semGraph <- semCoefs |> select(from, to, weight, p) |> igraph::graph_from_data_frame(directed = TRUE) |> tidygraph::as_tbl_graph()semGraph# A tbl_graph: 7 nodes and 9 edges
#
# A directed acyclic simple graph with 1 component
#
# Node Data: 7 × 1 (active)
name
1 distance
2 age
3 firesev
4 hetero
5 abiotic
6 cover
7 rich
#
# Edge Data: 9 × 4
from to weight p
1 1 2 -0.278 0.0079
2 1 4 0.346 0.0008
3 1 5 0.460 0
# ℹ 6 more rows4.2 绘图:结构方程模型路径图
# 字体customizedFont <- "DejaVu Sans"# 简单的 SEM 可视化semGraph |> ggraph(layout = "tree") + # sugiyama, circle, tree geom_edge_arc( strength = -0.2, aes( edge_color = ifelse(weight > 0, "+", "-"), edge_width = abs(weight), label = round(weight, 3) ), start_cap = circle(0.5), end_cap = circle(0.5), label_size = 3, angle_calc = "along", vjust = 1.5, arrow = arrow(angle = 60, length = unit(5, "pt")), family = customizedFont ) + geom_node_label( aes(label = name), size = 4, family = customizedFont ) + scale_edge_color_brewer(palette = "Set2", direction = -1) + scale_edge_width_continuous(range = c(0.2, 2)) + coord_equal(clip = "off") + theme_void( base_family = customizedFont, base_size = 18 ) + theme( legend.position = "none", strip.text = element_text(size = 12) )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.4.3 绘图:直接和间接效应
使用代码,计算直接和间接效应,
好处是,适用性强,只要将结果整理为 from、to、weight、p 数据框(如下)就可以使用,
缺点是,无法检验间接路径显著性(但是 semEff 包可以,见后文),
getEffects 函数见文末,另存为 getEffects.R,然后 source("getEffects.R") 调用之。
# 使用函数source("getEffects.R")effectData <- getEffects(semGraph)# 绘制效应effectData |> ggplot() + aes( x = from, y = value, fill = name ) + geom_col( color = "grey0", lwd = 1.5 ) + geom_label( aes(label = round(value, 3)), size = 3, position = position_stack(vjust = 0.5), family = customizedFont ) + scale_fill_manual( name = NULL, values = c("Direct" = "green4", "Indirect" = "grey100"), na.translate = FALSE ) + facet_grid(. ~ paste("Effects to:", to)) + coord_cartesian(ylim = c(-0.2, 0.6)) + theme_minimal( base_family = customizedFont, base_size = 16, base_line_size = 2 ) + theme( panel.grid.major.x = element_blank(), panel.grid.minor = element_blank(), strip.text = element_text(size = 14) )Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_col()`).Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_label()`).6 SEM:重抽样的直接和间接效应
用了另一个现成的包,如下,
因为是重抽样获得,所以与上述根据结果计算略有不同,
而且,可以检验间接效应是否显著!
# https://murphymv.github.io/semEff/articles/semEff.htmllibrary(semEff)# 重抽样 1000 次(并计时)system.time( pSemBoot <- semEff::bootEff( mod = pSem, R = 100, # 可设置为 1000 或以上 seed = 1, parallel = "multicore" ))user system elapsed
13.331 7.083 2.675获取效应,
# 获取效应semEff <- semEff::semEff(pSemBoot)# summary(semEff) # 此代码,获取所有效应,去掉注释查看效果# 各类效应semEff::getTotEff(semEff, "rich") |> round(3) # 总效应(Intercept) distance age hetero abiotic firesev
0.000 0.453 -0.053 0.301 0.219 -0.117
cover
0.267semEff::getDirEff(semEff, "rich") |> round(3) # 总效应(Intercept) distance hetero abiotic cover
0.000 0.233 0.301 0.219 0.267semEff::getIndEff(semEff, "rich") |> round(3) # 总效应(Intercept) distance age firesev
0.000 0.220 -0.053 -0.1176.1 重抽样:直接和间接效应
因为是重抽样获得,所以与上述根据结果计算略有不同。
# 结果存储在 effRes 变量中effRes <- semEff::getEff(semEff, "rich")# 提取直接、间接、总效应dEf <- effRes$rich$DirectindEf <- effRes$rich$IndirecttotalEf <- effRes$rich$Total# 转换为数据框funGetEff <- function(Tp, Ef) { data.frame(effType = Tp, Effect = Ef, row.names = names(Ef)) |> tibble::rownames_to_column(var = "Predictor")}dEfDf <- funGetEff("Direct", dEf)indEfDf <- funGetEff("Indirect", indEf)totalEfDf <- funGetEff("Total", totalEf)# 合并数据框combinedDf <- rbind(dEfDf, indEfDf, totalEfDf) |> dplyr::filter(Predictor != "(Intercept)")# 绘图combinedDf |> dplyr::filter(effType != "Total") |> ggplot() + aes( x = Predictor, y = Effect, fill = effType ) + geom_col( color = "grey0", lwd = 1.5 ) + geom_label( aes(label = round(Effect, 3)), size = 3, position = position_stack(vjust = 0.5), family = customizedFont ) + scale_fill_manual( name = NULL, values = c("Direct" = "green4", "Indirect" = "grey100"), na.translate = FALSE ) + facet_grid(. ~ paste("Effects to: rich")) + coord_cartesian(ylim = c(-0.2, 0.6)) + theme_minimal( base_family = customizedFont, base_size = 16, base_line_size = 2 ) + theme( panel.grid.major.x = element_blank(), panel.grid.minor = element_blank(), strip.text = element_text(size = 14) )6.2 重抽样:直接和间接效应的置信区间
bootEff <- semEff$`Bootstrapped Effects`[["rich"]]lapply(bootEff, head, 3)$Direct
(Intercept) distance hetero abiotic cover
[1,] -6.001212e-17 0.2489151 0.3069260 0.04232212 0.32140204
[2,] 1.596619e-17 0.3025537 0.4142906 0.12890911 0.10502234
[3,] -1.513694e-17 0.2867629 0.1815622 0.29990303 0.09868543
$Indirect
(Intercept) distance age firesev
[1,] -6.001212e-17 0.1586614 -0.044291734 -0.17155155
[2,] 1.596619e-17 0.1707798 -0.014439474 -0.03006690
[3,] -1.513694e-17 0.1651902 -0.006246364 -0.01870544
$Total
(Intercept) distance age hetero abiotic firesev
[1,] -6.001212e-17 0.4075764 -0.044291734 0.3069260 0.04232212 -0.17155155
[2,] 1.596619e-17 0.4733334 -0.014439474 0.4142906 0.12890911 -0.03006690
[3,] -1.513694e-17 0.4519531 -0.006246364 0.1815622 0.29990303 -0.01870544
cover
[1,] 0.32140204
[2,] 0.10502234
[3,] 0.09868543
$Mediators
age hetero abiotic firesev cover
[1,] 0.006250031 0.13073737 0.02167396 -0.038041703 -0.20959325
[2,] 0.002024480 0.10452326 0.06423202 -0.012414993 -0.04248189
[3,] 0.001294425 0.03769445 0.12620135 -0.004951939 -0.02365738# 将 list 中的每个矩阵转换为数据框,保留其名字 names;并合并所有数据框 bind_rowsbootEffDf <- map(bootEff, as.data.frame) |> map2(names(bootEff), ~ mutate(.x, effType = .y)) |> bind_rows() |> dplyr::select(-`(Intercept)`)# 计算每个路径的均值和 95% 置信区间summaryStats <- bootEffDf |> # 将数据转换为长格式,以便进行计算 pivot_longer( cols = -effType, names_to = "Predictor", values_to = "tempEstimate" ) |> # 按分组 .by 计算均值和 95% 置信区间 summarise( Estimate = mean(tempEstimate, na.rm = TRUE), ci95Lower = quantile(tempEstimate, 0.025, na.rm = TRUE), ci95Upper = quantile(tempEstimate, 0.975, na.rm = TRUE), .by = c(effType, Predictor) )# 绘图summaryStats |> filter(!effType %in% c("Total", "Mediators")) |> ggplot() + aes( x = Predictor, y = Estimate, ymin = ci95Lower, ymax = ci95Upper, color = effType ) + geom_col( width = 0.9, lwd = 1.5, fill = "white", position = position_stack(vjust = 0.5) # position = position_dodge(1) ) + # geom_errorbar( # width = 0, # lwd = 0.5, # position = position_dodge(1), # show.legend = FALSE # ) + geom_text( label = "✱", size = 2, position = position_stack(vjust = 0.5), # position = position_dodge(1), family = customizedFont ) + scale_color_manual( name = NULL, values = c("Direct" = "green4", "Indirect" = "grey0"), na.translate = FALSE ) + facet_grid(. ~ paste("Effects to: rich")) + coord_cartesian(ylim = c(-0.2, 0.6)) + theme_minimal( base_family = customizedFont, base_size = 16, base_line_size = 2 ) + theme( panel.grid.major.x = element_blank(), panel.grid.minor = element_blank(), strip.text = element_text(size = 14) )Warning: Stacking not well defined when not anchored on the axisWarning: Removed 5 rows containing missing values or values outside the scale range
(`geom_col()`).Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text()`).