R语言：结构方程模型简单遍历所有间接路径(借助 igraph包)

此文中的结构方程模型是使用 piecewiseSEM包构建的，如果使Lavaan、SEMinR或者plspm包构建的结构方程模型也可以尝试，只需要得到每条路径的系数即可。因为重点是igraph包计算节点到节点之间的所有路径。如果对象可以转化为igraph对象就都可以实现，本文使用的例子是个结构方程模型变量SWE.spring.psem ，长这样：

如果你 summary 这个变量就会得到他的所有信息。

> (new.summary <- summary(SWE.spring.psem, .progressBar = F))

Structural Equation Model of SWE.spring.psem 

Call:
  climate ~ geography
  soil ~ climate + geography + biophysical
  biophysical ~ geography + climate
  total_effects_spring ~ climate + biophysical + soil

    AIC
 -414667.352

---
Tests of directed separation:

                          Independ.Claim Test.Type    DF Crit.Value P.Value 
  total_effects_spring ~ geography + ...      coef 46566    -0.1314  0.8954 

--
Global goodness-of-fit:

Chi-Squared = NA with P-value = NA and on 1 degrees of freedom
Fisher's C = 0.221 with P-value = 0.895 and on 2 degrees of freedom

---
Coefficients:

              Response   Predictor Estimate Std.Error    DF Crit.Value P.Value Std.Estimate    
               climate   geography   0.9126    0.0096 47184    95.4825       0       0.4098 ***
                  soil     climate  -8.9855    0.1013 46567   -88.6675       0      -0.3270 ***
                  soil   geography  45.7612    0.2426 46567   188.6293       0       0.7480 ***
                  soil biophysical -17.1485    0.1646 46567  -104.1809       0      -0.3886 ***
           biophysical   geography   0.5216    0.0064 47183    81.9163       0       0.3763 ***
           biophysical     climate   0.0922    0.0028 47183    32.8412       0       0.1481 ***
  total_effects_spring     climate   0.8882    0.0130 46567    68.4619       0       0.3168 ***
  total_effects_spring biophysical   0.5661    0.0209 46567    27.1009       0       0.1257 ***
  total_effects_spring        soil  -0.0052    0.0005 46567   -10.9955       0      -0.0505 ***

  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05

---
Individual R-squared:

              Response method Marginal Conditional
               climate   none     0.18        0.27
                  soil   none     0.33        0.65
           biophysical   none     0.17        0.46
  total_effects_spring   none     0.16        0.16

我需要计算Geography到total_effects_spring的所有路径，再把路径系数相乘，计算这些路径的效应的加和什么方法最快速？一个个数效率太低下了。使用igraph包可以根据图论直接给你所有路径。首先先得到每个简单连接路径的系数，这个可以根据 summary 的结果进行整理得到，代码如下：

library(igraph)
kk1<-summary(SWE.spring.psem,.progressBar = F)
kk2 <- kk1$coefficients[,-(3:6)]
kk3 <- as.data.frame(kk2)
kk4 <- kk3[!grepl("~~", kk3$Predictor), ] %>% dplyr::select(2,1,4)
kk4
    Predictor             Response Std.Estimate
1   geography              climate       0.4098
2     climate                 soil      -0.3270
3   geography                 soil       0.7480
4 biophysical                 soil      -0.3886
5   geography          biophysical       0.3763
6     climate          biophysical       0.1481
7     climate total_effects_spring       0.3168
8 biophysical total_effects_spring       0.1257
9        soil total_effects_spring      -0.0505

我们就可以根据这个系数数据狂建立一个 igraph 对象，下面这个函数就是根据数据狂创建图对象的函数。

> g = graph_from_data_frame(kk4, directed=TRUE)
> plot(g)

这个对象建立完毕，可以对比一下刚刚那张图看看有不有丢下什么路径，接下来是计算 geography 到 total_effects_spring的所有路径，可以使用igraph::all_simple_path函数

> geography_paths <- all_simple_paths(g, "geography","total_effects_spring")
> geography_paths
[[1]]
+ 5/5 vertices, named, from 6a18791:
[1] geography            climate              biophysical          soil                
[5] total_effects_spring

[[2]]
+ 4/5 vertices, named, from 6a18791:
[1] geography            climate              biophysical          total_effects_spring

[[3]]
+ 4/5 vertices, named, from 6a18791:
[1] geography            climate              soil                 total_effects_spring

[[4]]
+ 3/5 vertices, named, from 6a18791:
[1] geography            climate              total_effects_spring

[[5]]
+ 4/5 vertices, named, from 6a18791:
[1] geography            biophysical          soil                 total_effects_spring

[[6]]
+ 3/5 vertices, named, from 6a18791:
[1] geography            biophysical          total_effects_spring

[[7]]
+ 3/5 vertices, named, from 6a18791:
[1] geography            soil                 total_effects_spring

一共有 7 条路径，由图知道这七条都是间接影响路径（因为 geography 没有直接连接到 total_spring_effect的路径），所以只需要遍历每条7 条路径计算其中简单路径系数相乘就行了。代码如下：

Vname <- vertex_attr(g)
geography_total_effects = 0
for(paths in geography_paths){
  basis = 1
  for(i in 1:(length(paths)-1)){
    from = Vname$name[paths[i]]
    to = Vname$name[paths[i+1]]
    basis = basis * kk4$Std.Estimate[kk4$Predictor==from&kk4$Response==to]
  }
  geography_total_effects = geography_total_effects + basis
}
geography_indirect_effects = geography_total_effects

其中vertex_attr函数的目的是为了得到每个节点的名字。

> Vname <- vertex_attr(g)
> Vname
$name
[1] "geography"            "climate"              "biophysical"         
[4] "soil"                 "total_effects_spring"

得到对象geography对 total_effects_spring的总效应是 0.1623233

> geography_indirect_effects
[1] 0.1623233

学会了你就动手试试吧。
因为本人的博客全都是很细节的东西，可能有这方面需求的人和看的人也不多。麻烦点个赞，谢谢。完了还可以出更多 R 语言的细节技术。下一期我出一个如何《如何传递一个字符串作为函数参数，然后还作为列名调用》下下期我出一个《如何R 语言并行 for 循环，并且把出错的循环自动跳过》。

👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻