Python实现PLS-PM：结构方程模型中小样本分析的新选择

简介：偏最小二乘路径建模算法的Python语言实现，定位为结构方程建模与因果预测分析的轻量级工具包，面向数据科学研究者、统计建模人员及需要处理中小样本或非正态数据的开发者。该程序源自R语言经典包的移植，并吸收其他扩展模块的功能，采用基于相关性的估计算法，允许以潜变量和显变量构建复杂的因果或预测模型，特别适合探索性研究场景。资源共六十三份文件，包括二十三份Python源码模块、三十三份CSV测试数据集，以及文档、配置文件、构建脚本等，压缩包仅九十九KB，代码与数据分离，便于理解算法原理并运行回归验证。目前已有两千六百三十六人学习下载，说明其在结构方程建模学习与实践中具有较高参考价值。读者可借此快速搭建分析流程，查看核心逻辑，也可基于测试数据复现文档示例，用于教学或二次开发。

1. 偏最小二乘路径建模（PLS-PM）的Python 3实现：中小样本结构方程终于不用死守R了

偏最小二乘路径建模（PLS-PM）是一种基于相关性的结构方程建模算法，它先用显变量估出潜变量得分，再在潜变量之间做路径回归，天然适合探索性研究、预测导向模型和中小样本场景，而且不要求数据满足多元正态假设。R的plspm包在这个领域几乎成了事实标准，但Python生态里一直缺一个等价的落地方案，直到这个移植包出现。这个项目把R的plspm完整迁移到Python 3，还吸收了R包seminr的附加功能，并且用回归测试锁定了和R结果的一致性。这篇笔记会从算法主干拆起，带你把配置、迭代、bootstrap和单维度检验整个流程跑通，最后把移植过程中最常翻车的地方一次性说清楚。

2. 拆开plspm包：从文件结构看懂PLS-PM的算法主干

2.1 潜变量、外模型与内模型：先建立共同词汇

用过AMOS、lavaan这类协方差型SEM工具的人，多少都经历过这样的场面：样本量只有一百多份，某个潜变量的指标稍多一点，模型就不收敛，或者修正指数改到最后面目全非。PLS-PM的思路完全不同，它不追求复现总体协方差矩阵，而是通过迭代的方式直接逼近潜变量得分，再回到潜变量之间做普通最小二乘回归。这个设计决定了它对样本量和数据分布都不敏感，才是真正适合业务探索场景的算法。

在理解代码前，先要建立三组词。显变量（indicators）是直接测到的题项，比如“这家店的服务让我满意”的打分；潜变量（latent variables）是测不到的抽象概念，比如满意度。外模型（outer model）描述一个潜变量和它名下那组显变量的关系，内模型（inner model）描述潜变量之间的影响关系。整个PLS-PM就是在“外模型迭代算权重”和“内模型回归算路径”之间来回切换。

对应到包的源码，这个分工非常清楚。outer_model.py负责外模型的载荷（loading）和权重输出，inner_model.py与inner_summary.py负责路径系数和R²的整理，weights.py是迭代中更新外权重的核心，estimator.py负责收敛后的最终估计。mode.py声明测量模式，scheme.py控制内模型加权方案，scale.py处理数据标准化与非度量数据转换。文件结构基本复刻了R plspm的对象体系，从R转过来的同学不会有陌生感。

把概念先立住是有实际意义的，因为排查翻车现场时，你得能快速定位问题在哪一层。路径系数方向不对，问题一定出在内模型或路径矩阵；载荷出现负值，优先怀疑测量模式和反向计分；迭代不收敛，则要回到权重更新和预处理环节。很多人一上来就盯着结果看，反而丢了排查的抓手。

2.2 模式A与模式B：反映型与形成型测量对应两套权重算法

测量模式的选择是PLS-PM里最容易被忽略、又最能改变结果的一个参数。mode.py里明确区分了A和B两种模式，它们的权重计算逻辑完全不同。

模式A是反映型（reflective）测量，潜变量是因，显变量是果。“忠诚度”高的顾客，大概率同时打高“再购意愿”“推荐意愿”和“价格容忍度”，这三道题的共性就是忠诚度。模式A下，外权重用指标与当前潜变量得分的协方差来计算：

# 模式A（反映型）：外权重 = 块内指标与潜变量得分的协方差 w_blk = X_blk.T @ Z_blk / n_obs

这里X_blk是标准化后的指标矩阵，Z_blk是内模型加权得到的潜变量中间估计。这个公式本质上是用相关性加权，一个题项与潜变量得分相关性越高，它在构造潜变量时的话语权就越大。模式A的权重一般比较稳定，也很少出现大幅负值，前提是你没有把反向题忘了反转。

模式B是形成型（formative）测量，显变量是因，潜变量是果。比如“购买意愿”由“价格敏感度”“品牌认同”“功能需求”共同构成，而不是反过来。这时权重就不能用简单协方差了，因为多个指标之间可能有重叠信息，必须用多元回归来分离各自的贡献：

# 模式B（形成型）：外权重 = 最小二乘回归系数 w_blk = np.linalg.lstsq(X_blk, Z_blk, rcond=None)[0]

模式B的权重代表在控制其它指标后，该指标的增量贡献，所以它更敏感，也可能出现负号。如果这个块内指标共线性严重，回归系数会被放大到离谱，迭代还会震荡。碰到这种情况，先检查块内相关系数矩阵，再看要不要换成模式A，或者删掉一个冗余指标。

实际项目里怎么定？有一个经验可以抄：态度量表、感知量表这类“潜变量影响回答”的场景，用模式A；绩效指标、能力指标这类“几个维度拼出一个总分”的场景，用模式B。拿不准时，默认模式A，因为多数论文里的测量模型都是反映型，而且模式A迭代更稳健。如果你在同一模型里混用两种模式，建议在文档里把每个块的模式写清楚，不然后面解读bootstrap结果时会绕晕。

2.3 迭代主循环：权重更新与收敛判断是怎么串起来的

把整个PLS-PM压缩成一段逻辑，其实就是“内部估计 → 外部加权 → 更新潜变量得分”循环往复。主入口plspm.py干的就是这件事。我按这个包的结构整理了一份核心流程示意代码，可以直接对照阅读：

import numpy as np def plspm_core(X, blocks, path_matrix, modes, scheme='path', maxiter=100, tol=1e-6): """ 简化版 PLS-PM 迭代主循环 X 已标准化的显变量矩阵，行是样本，列是指标 blocks 每个潜变量对应的指标列索引列表 path_matrix 路径矩阵，行是被解释变量，列是解释变量 modes 每个潜变量的测量模式，'A' 或 'B' scheme 内部加权方案：'centroid' / 'path' / 'factor' """ n_obs, n_ind = X.shape n_lv = len(blocks) # 潜变量得分初始值：块内指标均值 Y = np.zeros((n_obs, n_lv)) for j, cols in enumerate(blocks): Y[:, j] = X[:, cols].mean(axis=1) W = np.ones(n_ind) / n_ind for it in range(maxiter): Y_old = Y.copy() # Step 1: 内部估计。对每个潜变量 j，找到邻居潜变量， # 按内部加权方案把邻居得分加权成 Z_j Z = np.zeros_like(Y) for j in range(n_lv): neighbors = [i for i in range(n_lv) if path_matrix[j, i] == 1 or path_matrix[i, j] == 1] if not neighbors: Z[:, j] = Y[:, j] continue for i in neighbors: r = np.corrcoef(Y[:, j], Y[:, i])[0, 1] if scheme == 'centroid': gamma = np.sign(r) elif scheme == 'factor': gamma = r elif scheme == 'path': gamma = 1.0 if path_matrix[j, i] == 1 else 0.0 Z[:, j] = Z[:, j] + gamma * Y[:, i] # Step 2: 外部加权。根据测量模式更新外权重 for j, cols in enumerate(blocks): X_blk = X[:, cols] if modes[j] == 'A': W[cols] = X_blk.T @ Z[:, j] / n_obs elif modes[j] == 'B': W[cols] = np.linalg.lstsq(X_blk, Z[:, j], rcond=None)[0] # Step 3: 更新潜变量得分并做标准化 for j, cols in enumerate(blocks): Y[:, j] = X[:, cols] @ W[cols] std = Y[:, j].std() if std > 1e-12: Y[:, j] = (Y[:, j] - Y[:, j].mean()) / std # 符号校正：防止某轮迭代权重翻转导致方向震荡 for j in range(n_lv): if np.corrcoef(Y[:, j], Y_old[:, j])[0, 1] < 0: Y[:, j] = -Y[:, j] if np.max(np.abs(Y - Y_old)) < tol: break return Y, W

这段代码的每一行都能对应到包里的模块。scheme.py实现的是Step 1里的三种加权策略；weights.py实现的是Step 2的权重更新；scale.py管理Step 3的标准化逻辑。scheme参数的差异值得单独说一下：centroid只取相关系数的符号，最保守也最稳；factor直接用相关系数加权，信息量最足但对离群值敏感；path严格按路径矩阵方向加权，没有指向关系就不给权重。业务模型中有明确的中介或因果假设时，我一般选path，因为它给出的路径系数更贴近你设计模型时的理论意图。

有一个细节容易踩坑：潜变量得分的符号是不确定的，同一组权重乘出来可能是正版本也可能是负版本。上面代码里的符号校正，就是为了保证每一轮迭代的方向一致。包的实现里对这一点处理得比较仔细，但你自己写变体时要记得补上，否则就会看到迭代明明收敛了，路径系数却一轮正一轮负的诡异现象。

3. 跑通第一个模型：路径矩阵、blocks与参数设置

3.1 一个最小可运行的三潜变量模型

理论落回代码，先跑通一个最小模型。我常用感知质量 → 满意度 → 忠诚度这条链式结构来验证工具链，三个潜变量各配三个指标，数据量只有一百来行。这个样本量放到协方差型SEM里是灾难，但在PLS-PM里完全够用。

import pandas as pd from plspm import plspm # 1. 路径矩阵：行=被解释变量，列=解释变量，1 表示列影响行 path_matrix = pd.DataFrame( [[0, 0, 0], [1, 0, 0], [0, 1, 0]], index=['QUAL', 'SAT', 'LOY'], columns=['QUAL', 'SAT', 'LOY'] ) # 2. blocks：每个潜变量对应的指标列名 blocks = [ ['QUAL_1', 'QUAL_2', 'QUAL_3'], ['SAT_1', 'SAT_2', 'SAT_3'], ['LOY_1', 'LOY_2', 'LOY_3'] ] # 3. 测量模式：三个潜变量都用反映型 modes = ['A', 'A', 'A'] # 4. 执行 PLS-PM res = plspm(df, path_matrix, blocks, modes, scheme='path', scaled=True)

跑这段代码前，先确认df里没有缺失值。这个包和R plspm一样，没有内置缺失值插补，直接把含NaN的DataFrame传进去，某些计算会静默出错或者产出一堆NaN结果，定位起来非常痛苦。我的习惯是进入模型前先执行一遍df.isna().sum()，该删行删行，该填充填充。

scaled=True的作用是在迭代前把所有显变量做z-score标准化。绝大多数业务场景都该开着，因为指标刻度经常差几个量级，比如一个1到7的李克特题和一个0到100的百分比题放在同一个块里，不标准化，大刻度指标会变成幽灵权重，直接压过其它指标。如果你确实有理论理由不标准化，比如所有指标天生同量纲，那也要自己确认量纲一致再关。

3.2 关键参数逐一拆解：scheme、scaled、maxiter、tol

这个包的参数命名延续了R plspm的习惯，看一眼就知道对应关系。完整调用长这样：

res = plspm(df, path_matrix, blocks, modes, scheme='path', # 'centroid' / 'path' / 'factor' scaled=True, # 是否标准化显变量 maxiter=300, # 最大迭代次数 tol=1e-6, # 收敛阈值 boot_val=False) # 是否执行 bootstrap

参数不多，但每个都有实际讲究，整理成表方便对照：

scheme是内模型加权方式，前面已经分析过三种方案的差异。补充一点：如果你模型里有中介变量，path方案通常比centroid更合适，因为它只让路径矩阵中明确写出的因果关系参与内部加权，不会引入理论上不该存在的间接关联噪声。

maxiter和tol共同控制循环何时停止。tol调到1e-7通常只会增加几十轮无效迭代，结果基本没变化；maxiter才是那个需要根据你数据情况调整的参数。特别是模式B模型，外权重靠回归估计，迭代起来比模式A慢得多，100轮经常不够。我第一次跑一个四潜变量模型时，警告信息显示迭代次数用满但收敛差一点点，把maxiter调到300后就安静了。这种小样本低成本调整，没必要死守默认值。

3.3 解读第一份输出：载荷、路径系数与R²

跑完后，res对象里包含外模型结果、内模型结果、潜变量得分和R²，结构上贴近R plspm的返回对象。我习惯先打印这三样：

print(res.outer_model) # 外模型：载荷与权重 print(res.inner_model) # 内模型：路径系数 print(res.r_squared) # 内生潜变量的 R²

典型输出大致是这样的结构：

OUTER MODEL weight loading QUAL_1 0.352 0.821 QUAL_2 0.338 0.774 QUAL_3 0.317 0.733 SAT_1 0.401 0.863 SAT_2 0.327 0.745 SAT_3 0.316 0.709 INNER MODEL Estimate QUAL -> SAT 0.782 SAT -> LOY 0.653 R² SAT 0.611 R² LOY 0.583

解读结果我分三层。第一层看载荷（loading），它是指标与潜变量的简单相关，反映“这个题项到底有没有测到目标概念”。载荷低于0.7要警惕，低于0.5基本可以淘汰。第二层看权重（weight），它代表指标在构造潜变量得分时的实际贡献。如果某指标载荷不低但权重异常低，说明它的信息被同块其它指标覆盖了，典型表现就是”这块指标测的东西高度重叠“。第三层才看路径系数和R²，路径系数是潜变量得分之间的标准回归系数，R²说明内生潜变量被前因变量解释的程度。本例中LOY被SAT解释了58.3%的方差，在管理研究里属于中等偏上的水平。

有一个细节经常被忽略：内模型路径系数是在潜变量得分上回归得到的，因此对潜变量的尺度敏感。全程保持scaled=True时，潜变量得分尺度的可比性有保障，路径系数可以直接比较相对大小；如果你关闭了标准化，路径系数的绝对数值就失去了解释意义，只能比较方向。

4. bootstrap与置信区间：给点估计补上稳健性证据

4.1 bootstrap.py 在做什么

PLS-PM本身不做显著性检验，它给的是点估计。路径系数是0.782，这个数字到底稳不稳，置信区间是多宽，都需要bootstrap来回答。bootstrap.py复刻了R plspm中的bootstrap逻辑：对原始数据做有放回抽样，每次抽取与原始样本等量的行，重新跑一遍完整的PLS-PM，重复若干次后，累积出每个参数的分布。

import numpy as np def bootstrap_plspm(df, path_matrix, blocks, modes, br=200, seed=None): """对原始数据行做有放回重抽样，逐次重估模型并收集路径系数。""" if seed is not None: np.random.seed(seed) n = len(df) boot_params = [] for b in range(br): idx = np.random.choice(n, size=n, replace=True) boot_df = df.iloc[idx].reset_index(drop=True) boot_res = plspm(boot_df, path_matrix, blocks, modes, scheme='path', scaled=True) boot_params.append(boot_res.path_coefs.copy()) return np.array(boot_params) # shape: (br, n_lv, n_lv)

这段代码的关键在于，重采样的对象是原始数据的行，不是潜变量得分。只有这样，重采样样本才能保留原始数据里显变量之间的相关结构，路径系数的分布才有统计意义。如果拿潜变量得分去重采样，等于把测量关系人为切断了，结果就是区间窄得离谱，误导性极强。

4.2 运行bootstrap并解读置信区间

在包里，bootstrap不是独立函数，而是主流程的一个开关，使用起来和R plspm的boot.val=TRUE完全对应：

res = plspm(df, path_matrix, blocks, modes, scheme='path', scaled=True, boot_val=True, br=500, seed=42) print(res.boot.path_coefs) print(res.boot.path_coefs_bounds) # 置信区间边界

500次重采样输出的路径系数置信区间大致长这样：

perc.025 perc.975 QUAL -> SAT 0.702 0.851 SAT -> LOY 0.548 0.714

判断规则非常直接：区间不含0，就认为该路径在统计意义上成立。需要特别强调的是，这是百分位法置信区间，不是正态近似，所以它允许分布不对称。比如QUAL → SAT的区间是0.702到0.851，右端压缩更紧，说明重采样分布有轻微左偏。这种不对称恰恰是非正态数据下bootstrap的价值所在。

至于br取多少，我的经验是按用途分档。快速探索用100到200次够用；正式写报告至少500次；如果模型里有中介效应检验，建议升到1000次，因为间接效应的分布通常更偏斜，需要更多重采样次数才能稳定住两端的百分位点。

4.3 非度量数据：scale.py与排序数据的处理

问卷里除了连续量表，还有大量顺序尺度和二分变量。R plspm对这类数据会在内部改用Spearman相关，以确保相关矩阵的稳健性。这个包的scale.py承担了对应的转换职责，而test_regression_nonmetric.py这个测试文件的存在，说明作者专门验证过非度量数据场景下与R结果的一致性。

触发条件是数据不是等距尺度。李克特5级量表通常可以直接当连续数据处理，但“从不/偶尔/经常/总是”这类纯顺序标签，就需要先编码成有序整数，并且留意迭代计算中把等级差异强行当作等距可能带来的偏差。包的实现一般会在标准化前先做秩转换，再进入标准化的流程。如果你手头有这类数据，跑模型前瞄一眼scale.py中走的处理分支，确认它选择的是Spearman分支而不是默认的Pearson分支。这一步能省掉不少“为什么载荷全部缩水到0.3”的困惑。

5. 避坑指南：从R移植到Python最常栽的几个跟头

实操下来，我在这套工具上栽过的跟头可以归纳成五条，每条都附上现象、原因和解决办法，对照排查比盲猜有效率得多。

5.1 路径矩阵方向搞反，路径系数全变号

现象：跑出来的路径系数方向和你理论假设完全相反，而且R²意外地高。

原因：R plspm的path_matrix是行代表被解释变量、列代表解释变量，1表示列影响行。很多从lavaan或AMOS切换过来的人，会把行和列的语义搞颠倒，传进去的其实是转置矩阵。

解决：打印path_matrix手工核对。三变量链式模型的非0元素应该在(2,1)、(3,2)两个位置，其余全是0；如果你看到非0元素跑到了(1,2)、(2,3)，直接path_matrix = path_matrix.T再跑一遍。这个问题在R里也经常发生，并不是Python移植版独有的毛病，但移植版没有R里那套兜底警告，静默出错的可能性更高。

5.2 blocks列顺序错位，潜变量张冠李戴

现象：外模型载荷看起来正常，但QUAL的指标跑到了LOY的名下，业务结论完全错位。

原因：blocks列表的顺序要和path_matrix的行列索引对应。比如path_matrix用了QUAL、SAT、LOY，blocks却被写成了SAT、LOY、QUAL的顺序。这个错在代码审查里极难用肉眼发现，因为所有数字都对得上，只有列名对不上。

解决：跑完后第一件事就是打印外模型载荷，核对每个潜变量名下的列名。我习惯提前建一个映射表：

block_map = {lv: blocks[i] for i, lv in enumerate(path_matrix.index)} print(block_map)

这个动作我几乎每次都做，因为bootstrap结果再漂亮，也救不了一个指标归属就错了的模型。

5.3 测量模式选错，载荷正负翻盘

现象：外模型载荷出现成片负值，或者某个块内权重正负交替得毫无规律。

原因：模式选错。比如某块实际是形成型测量，你按模式A跑，协方差加权就会把负向指标硬拉成负权重。另一种常见情况是反向计分题没有做反转，导致整个块的相关结构反了。

解决：先回到问卷设计确认测量性质，再确认是否处理过反向题。正向处理后如果还出现大面积负载荷，就测试切换模式B。辅助手段是用unidimensionality.py跑一遍单维度检验，块内出现负特征值往往是模式选择错了的信号，这个信号比看载荷更早也更客观。

5.4 不收敛时，先查数据预处理再调maxiter

现象：警告说迭代次数用完，或者同一份数据跑两次，结果会有微小差异。

原因：缺失值、零方差指标、模式B块内强共线性，这三个因素按概率排序。lstsq遇到含NaN的矩阵会静默产生垃圾结果，不报错，所以定位难度很大。

解决：严格按顺序排查。第一步df.isna().sum()查缺失；第二步逐列检查标准差，零方差列直接删；第三步看每个块内的相关矩阵，有没有超过0.95的强相关。全部排除后，再把maxiter调到300。这个顺序不要反，我见过有人上来就调maxiter到500，结果迭代仍然不稳定，最后发现是某个指标列里藏着NaN。

5.5 bootstrap不设种子，结果无法复现

现象：同一份数据两次跑bootstrap，点估计相同，置信区间边界却差0.03到0.05。

原因：bootstrap的随机性没有被固定。

解决：固定seed参数，并在分析记录里写清楚seed的值。还有一个小细节，种子一定要放在整个bootstrap循环之前，不能每轮重采样都重置；否则500次重采样会变成500遍一模一样的样本，置信区间窄得让评审老师一眼识破。

提示：seed只保证同一环境下的可复现性，不要跨平台跨版本强行追求完全相同的结果。

6. 进阶验证：单维度检验与跨语言交叉核对

6.1 用 unidimensionality.py 检验每个潜变量的单维度

拿到新数据，第一步要回答的问题不是“路径系数有多大”，而是“这一组指标到底能不能合成一个潜变量”。unidimensionality.py做的就是这件事。它针对每个块单独检验内部一致性，输出通常包括Cronbach alpha、Dillon-Goldstein rho和第一特征值。经验阈值是alpha大于0.7、rho大于0.7、第一特征值大于1。

from plspm.unidimensionality import compute_unidimensionality uni = compute_unidimensionality(df, blocks, modes) print(uni)

输出大致是这个结构：

block alpha DG.rho eig1 QUAL 0.83 0.88 3.21 SAT 0.81 0.86 2.98 LOY 0.79 0.84 2.77

如果某个块alpha在0.6以下，第一特征值逼近1，说明这一组题目内部一致性很弱，硬要合成一个潜变量，后面所有路径系数都是空中楼阁。正确的处理是先拆块，或者删掉贡献最小的指标，再重新评估。这一步务必放在bootstrap之前，因为置信区间再好看也救不了一个块本身不成立的结构。

6.2 与R plspm交叉验证：最后的保底

这个包自带test_regression_metric.py、test_regression_bootstrap.py和test_regression_seminr.py，作者在移植时明显做了和R端的一对一回归对齐。但我的态度是，测试文件通过不代表你的数据不出问题，跨语言交叉验证仍然是最省钱的保底手段。

具体做法很简单：在R里用相同数据跑一遍plspm，拿到路径系数和载荷，再到Python里跑一遍，逐项对比。两者差异通常应该小于千分之一量级。如果发现某几个系数在小数点后第二位开始分叉，不要急着归因于实现差异，先回到数据本身查一遍。这类差异往往是极端值、几乎重复的样本或者编码不一致引起的，不是算法问题。有一次我排查了一个小时，最后定位到数据里有一条重复记录被标准化放大成了极端值，删掉后两边结果完全对齐，路径系数也回到了合理区间。

6.3 收尾：我养成的最后一道工序

这套流程走了大半年后，我给自己定了一个强制收尾步骤：每次跑完正式模型，先做单维度检验确认块结构成立，再固定seed跑500次bootstrap留好置信区间，最后再花五分钟用R交叉核对一遍关键路径系数。最后一次迁移项目里，正是第三遍比对时救了我——Python和R的路径系数方向一致，但其中一个块的第一特征值只有1.1，回过头检查发现那个维度混入了两道几乎重复的题目。删掉冗余题后，模型重新迭代，R²和bootstrap区间都发生了实质性变化。

从那以后，哪怕赶进度，我也会把这三步强制走一遍，不会直接跳去画最终结论的汇报图。这个包的价值不仅在于把R的功能搬到了Python里，更在于它让整个PLS-PM流程有了可复现、可验证的完整链条。希望这篇拆解能帮你少走几趟弯路，把时间花在业务结论而不是怀疑工具上。

本文还有配套的精品资源，点击获取