0 引言
1 需要的包
2 需要的数据
3 使用 piecewiseSEM 拟合
4 结果的对比引言
在生态学统计分析中，结构方程模型（SEM）是一个强大的工具。怎么正确使用它，非常重要。
在做连续变量的交互作用时，
先分别把变量标准化，然后再相乘；
把变量相乘，然后再标准化。
那么，到底哪种是正确的变量标准化方式？
以下是我对一些基本概念的理解，也是提供给初学者的一些背景信息，
连续变量的交互作用，就是连续变量的乘积的效应；
标准化，就是使得变量缩放到均值是 0，标准差是 1；
piecewiseSEM 是 R 语言的一个 SEM 工具包，默认提供了标准化路径系数。
1 包
library(dplyr)
library(piecewiseSEM)2 数据
2.1 不预先标准化数据
# keeley 数据
# 来自 piecewiseSEM 包
dat1 <- keeley
# colnames(dat1)
# 假设
# age 与 firesev 有交互作用
# 添加交互作用项 age 乘以 firesev
dat1$age_firesev <- dat1$age * dat1$firesev
colnames(dat1)[1] "distance" "elev" "abiotic" "age" "hetero"
[6] "firesev" "cover" "rich" "age_firesev"2.2 预先手动标准化数据
注意，此处不标准化外生变量 rich。
2.2.1 分别 scale 再相乘
dat2 <- keeley
# scale 函数要修改一下
# 使得返回的结果是数值型
scale_num <- function(x) {
x |>
scale() |>
as.numeric()
}
# 注
# 先运行计算 age_firesev_s2 的代码
# 因为 age 和 firesev 变量没有重命名
# 因此 age 和 firesev 会被 scale 计算结果覆盖
dat2$age_firesev_s1 <- scale_num(dat2$age) * scale_num(dat2$age) # 分别 scale 再相乘
dat2$age <- scale_num(dat2$age)
dat2$firesev <- scale_num(dat2$firesev)
dat2$distance <- scale_num(dat2$distance)
# dat2$rich <- scale_num(dat2$rich)
c(
mean(dat2$age) |> round(5),
sd(dat2$age) |> round(5)
)[1] 0 1c(
mean(dat2$firesev) |> round(5),
sd(dat2$firesev) |> round(5)
)[1] 0 1c(
mean(dat2$age_firesev_s1) |> round(5),
sd(dat2$age_firesev_s1) |> round(5)
)[1] 0.98889 1.38173注意，
age_firesev_s1的均值不是 0，标准差也不是 1。因此，与标准化的目的（和结果）不相符。
2.2.2 相乘再 scale
dat3 <- keeley
# 注
# 先运行计算 age_firesev_s2 的代码
# 因为 age 和 firesev 变量没有重命名
# 因此 age 和 firesev 会被 scale 计算结果覆盖
dat3$age_firesev_s2 <- scale_num(dat3$age * dat3$firesev) # 相乘再 scale
dat3$age <- scale_num(dat3$age)
dat3$firesev <- scale_num(dat3$firesev)
dat3$distance <- scale_num(dat3$distance)
# dat3$rich <- scale_num(dat3$rich)
c(
mean(dat3$age_firesev_s2) |> round(5),
sd(dat3$age_firesev_s2) |> round(5)
)[1] 0 1此处，
age_firesev_s2的均值是 0，标准差是 1。因此，是正确的标准化。
3 使用 piecewiseSEM 拟合
3.1 用原始数据拟合
keeley_psem_1 <- psem(
lm(age ~ distance, data = dat1),
lm(firesev ~ distance + age, data = dat1),
lm(rich ~ distance + firesev + age_firesev, data = dat1),
age_firesev %~~% age,
age_firesev %~~% firesev
)
dSep(keeley_psem_1, .progressBar = FALSE) Independ.Claim Test.Type DF Crit.Value P.Value
1 rich ~ age + ... coef 85 1.764489 0.08124205fisherC(keeley_psem_1) Fisher.C df P.Value
1 5.021 2 0.0813.2 用预先手动标准化的数据拟合
3.2.1 分别 scale 再相乘的拟合
keeley_psem_s1 <- psem(
lm(age ~ distance, data = dat2),
lm(firesev ~ distance + age, data = dat2),
lm(rich ~ distance + firesev + age_firesev_s1, data = dat2),
age_firesev_s1 %~~% age,
age_firesev_s1 %~~% firesev
)
dSep(keeley_psem_s1, .progressBar = FALSE) Independ.Claim Test.Type DF Crit.Value P.Value
1 rich ~ age + ... coef 85 -0.07810522 0.9379279fisherC(keeley_psem_s1) Fisher.C df P.Value
1 0.128 2 0.938可以看到，这种标准化，与默认方法给出的 dSep 和 ficherC 不一致！
如果，这种标准化方法是错的，那么 SEM 将给出完全误导的结论。
3.2.2 相乘再 scale 的拟合
keeley_psem_s2 <- psem(
lm(age ~ distance, data = dat3),
lm(firesev ~ distance + age, data = dat3),
lm(rich ~ distance + firesev + age_firesev_s2, data = dat3),
age_firesev_s2 %~~% age,
age_firesev_s2 %~~% firesev
)
dSep(keeley_psem_s2, .progressBar = FALSE) Independ.Claim Test.Type DF Crit.Value P.Value
1 rich ~ age + ... coef 85 1.764489 0.08124205fisherC(keeley_psem_s2) Fisher.C df P.Value
1 5.021 2 0.081可以看到，这种标准化，与默认方法给出的 dSep 和 ficherC 的结果一致！
4 结果的对比
4.1 标准化后数据的对比
png("s1-s2.png", width = 5, height = 3, units = "in", res = 600)
par(
family = "Prompt",
mfrow = c(1, 1),
mar = c(4, 4, 1, 1), # bottom, left, top, right
cex = 1.4
)
plot(
dat2$age_firesev_s1,
dat3$age_firesev_s2,
pch = 16,
cex = 0.6,
xlab = "s1: scale bef. mult.",
ylab = "s2: scale aft. mult."
)
abline(a = 0, b = 1, lty = 2)
dev.off()quartz_off_screen
2
可以看到，分别标准化再相乘（X 轴），与相乘再标准化（Y 轴）数据的性质已经完全不同。
4.2 SEM 拟合结果的对比
rbind(
coefs(keeley_psem_1) |>
select(
Response,
Predictor,
Std.Estimate, # 标准化的系数
Estimate, # 未标准化的系数
P.Value
) |>
mutate(Method = "default"),
coefs(keeley_psem_s1) |>
select(
Response,
Predictor,
Std.Estimate, # 标准化的系数
Estimate, # 未标准化的系数
P.Value
) |>
mutate(Method = "method 1"),
coefs(keeley_psem_s2) |>
select(
Response,
Predictor,
Std.Estimate, # 标准化的系数
Estimate, # 未标准化的系数
P.Value
) |>
mutate(Method = "method 2")
) |>
arrange(Predictor, Response, Method) |>
relocate(Method, Std.Estimate, Predictor, Response) |>
filter(stringr::str_detect(Predictor, "~", negate = TRUE)) |>
print() Method Std.Estimate Predictor Response Estimate P.Value
1 default 0.4081 age firesev 0.0537 0.0001
2 method 1 0.4081 age firesev 0.4081 0.0001
3 method 2 0.4081 age firesev 0.4081 0.0001
4 default -0.2164 age_firesev rich -0.0378 0.1190
5 method 1 -0.1423 age_firesev_s1 rich -1.5560 0.1005
6 method 2 -0.2164 age_firesev_s2 rich -3.2683 0.1190
7 default -0.2781 distance age -0.3959 0.0079
8 method 1 -0.2781 distance age -0.2781 0.0079
9 method 2 -0.2781 distance age -0.2781 0.0079
10 default -0.1644 distance firesev -0.0308 0.0967
11 method 1 -0.1644 distance firesev -0.1644 0.0967
12 method 2 -0.1644 distance firesev -0.1644 0.0967
13 default 0.4813 distance rich 0.8235 0.0000
14 method 1 0.4869 distance rich 7.3547 0.0000
15 method 2 0.4813 distance rich 7.2710 0.0000
16 default -0.0665 firesev rich -0.6080 0.6147
17 method 1 -0.2399 firesev rich -3.6232 0.0072
18 method 2 -0.0665 firesev rich -1.0046 0.6147可以看到，
age_firesev_s1得到的标准化系数不同！
结语
先标准化，再相乘，得到的变量，
均值不为 0，标准差也不为 1，与标准化的目的相悖；
它代入 SEM 后获得的结果与 piecewiseSEM 默认给出的
Std.Estimate不一致；很可能是错的？
本文拟合的只是一个玩具 SEM，这个 SEM 内部的变量关系是臆造的，仅仅作为演示。此外，本公众号的管理员不是统计专家，本文的观点仅供参考。请大家讨论。
附注
转发到票圈，或转发到科研群，关注本公众号，并留下邮箱，可获得本文对应的 Quarto 文档。
提供可复现的数据和代码，对于科学研究的质量可能非常重要，
提供一份统计分析报告，也许，在审稿阶段就可以阻止很多错误的发生。
使用 Quarto 文档（https://quarto.org/）以确保统计分析结果可以彻底复现，
Quarto 可以理解为一个写代码的记事本，可以一气呵成地，把代码和结果导出为 PDF 或者 HTML 等格式；
如果代码报错，Quarto 的导出步骤就会停止。所以，能被 Quarto 导出的，一般是完全可复现的代码（不过，可人为设置跳过某些代码）；
尽管，可复现，不代表分析一定正确；
本文即是使用 Quarto 生成。
