此文中的结构方程模型是使用 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 循环,并且把出错的循环自动跳过》。
👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻👻