ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

Python 3实现PLS-PM偏最小二乘路径建模:从原理到代码

Python 3实现PLS-PM偏最小二乘路径建模:从原理到代码 简介偏最小二乘路径建模PLS-PM算法的 Python 3 实现包面向结构方程建模研究者、数据科学从业者及相关课程学习者。PLS-PM 是一种基于相关性的结构方程建模算法用于通过潜在变量和显变量估计复杂因果或预测模型它特别适合探索性研究能够处理中小样本乃至大型数据集且不要求满足多元正态假设。资源包共 63 个文件以 Python 脚本23 个与 CSV 数据文件33 个为主体另含 Markdown 文档、配置文件、Makefile 等压缩后约 99KB源码将参数配置、模式设定、权重计算、模型估计、Bootstrap 检验等功能拆分到独立模块便于学习与二次开发。目前已有 2636 人学习/下载关注度较高配套测试用例覆盖度量与非度量数据场景并参考 R 语言 plspm 与 seminr 的设计思路可结合 R 流程理解算法细节快速用 Python 完成 PLS-PM 建模与对照验证。1. 为什么营销分析里会有人找 PLS-PM 的 Python 3 实现做客户满意度研究时经常遇到一个尴尬模型里有“感知质量”“期望”“价值”“满意度”几个抽象概念每个概念又靠三五个问卷题在测。你可以用线性回归硬解但回归默认自变量没有测量误差而且潜变量之间的路径关系没法一次性估出来。偏最小二乘路径建模 (PLS-PM) 正是为这种场景生的它用主成分思想把多个指标压缩成潜变量得分再用普通最小二乘估计路径系数小样本、非正态、指标之间相关高都能扛这在社科、管理学和市场营销里几乎是标配。可惜最成熟实现是 R 的plspm包Python 3 生态里要么只有PLSRegression要么是零碎代码所以我干脆写了一个最小可复现的纯 NumPy 版本下文从原理讲到踩坑。2. PLS-PM 的算法结构先搞清楚潜变量得分是怎么迭代算出来的如果你以前跑过线性回归PLS-PM 给你的第一感觉是“这到底在优化什么”。回归是显式优化残差平方和PLS-PM 的优化目标不直接出现它是在一个迭代规则下收敛到一组权重。所以理解算法的最好方式不是看目标函数而是顺着每次迭代里数据怎么流动看一遍。2.1 结构模型与测量模型先给潜变量画一张路径图在动手写代码前每个潜变量必须知道三件事自己有哪些显变量作为测量指标结构模型里它是原因还是结果测量模式是反映型还是构成型。这三点在理论上决定收不收敛、收敛到哪。反映型测量潜变量Mode A的假设是所有显变量是同一概念的表现比如满意度的三个题项删掉任何一个剩下的仍然测的是满意度构成型测量潜变量Mode B的假设是显变量按一定逻辑合成概念比如“社会经济地位”由收入、学历、职业构成删掉“职业”概念本身就从“地位”变成了“收入水平”。在实际问卷里绝大多数李克特题项是反映型所以许多入门资料默认 Mode A但一旦遇到“组成型”指标仍用 Mode A路径系数会严重偏移。代码上我习惯用字典而不是配置类来定义模型。测量模型字典的 key 是潜变量名value 是它对应的显变量列名列表结构模型字典的 key 是某潜变量结果value 是它的直接原因潜变量列表。这个约定和 R 的 plspm 包类似但更简单换数据集时只需要改字典。2.2 迭代四步外估计、内估计、更新权重、收敛判定假设已经有了一批标准化后的显变量矩阵 X。对每个潜变量 j用 w_j 表示它的权重向量Y_j 表示它的得分向量。算法在下面四步之间循环第一步外估计。Y_j X_j w_j然后对 Y_j 做标准化。这一步把所有显变量按当前权重线性合成生成一个临时潜变量。第二步内估计。根据结构模型找出与 j 有直接路径关系的邻居计算 Y_j 与每个邻居 Y_k 的相关系数 r_{jk}再用这些相关系数作为权重组装邻居得分得到 Z_j sum(r_{jk} * Y_k)。如果没有邻居Z_j Y_j。这一步让结构模型的信息反哺测量模型。第三步更新权重。把 Z_j 当因变量X_j 当自变量做一次回归。Mode A 下用一元回归每个显变量单独与 Z_j 回归去掉截距后的系数就是新权重Mode B 下用多元回归把 X_j 整体对 Z_j 回归回归向量就是新权重。随后把权重归一化到模长 1。第四步收敛判定。比较新旧权重的差当平均变化小于阈值时停止。四个步骤里外估计和更新权重是测量模型的更新内估计是结构模型的约束。所以即使数据不变只要结构模型变了最后权重也会变这是 PLS-PM 和直接取主成分的关键差异。为了便于递归理解可以把核心循环写成下面这样# 注意这是伪代码只展示每次迭代的顺序 for epoch in range(max_iter): for lv in lvs: scores[lv] norm(X[lv] weights[lv]) for lv in lvs: inner[lv] sum(corr(scores[lv], scores[nb]) * scores[nb] for nb in neighbors(lv)) for lv in lvs: weights[lv] norm(X[lv].T inner[lv]) # Mode A if max(weights - old_weights) tol: break这段伪代码里norm()表示减均值除标准差X[lv]是当前潜变量下的显变量矩阵。很多初学者以为 PLS-PM 的收敛对象是潜变量得分其实不是收敛对象是权重。权重变了得分自然变得分变权重又跟着变。2.3 为什么这套办法能处理非正态和小样本与协方差结构方法的区别结构方程还有另一个大家族叫协方差结构CB-SEM代表工具是 AMOS、LISREL。CB-SEM 的解是从样本协方差矩阵拟合模型参数要求数据近似多元正态、样本量通常要 200 以上。PLS-PM 走的是成分提取路线它没有对噪声分布做强假设每一步都在做线性回归和相关系数估计这两个操作对尖峰厚尾数据有一定抵抗力。小样本时CB-SEM 经常出现不收敛或负方差而 PLS-PM 的迭代本质上是在做一系列普通最小二乘样本数大于某个潜变量下显变量个数就能跑。但这里有个边界要说清PLS-PM 对缺失值没有天然免疫对极端多重共线性也敏感。所谓“小样本也能跑”是指 50-100 个样本能给出数值解但不代表标准误可靠所以后文必须用 Bootstrap 补显著性检验。3. 用 Python 3 从零实现 PLS-PM一个最小可复现的类现在进入可直接抄作业的部分。我写类时避免依赖statsmodels、sklearn这类重量级库只用numpy和pandas。这样跑在干净环境里也不容易因为版本问题翻车。整个类的核心大约 150 行够做研究分析和业务场景验证。3.1 数据组织与模型定义用字典描述测量关系和路径关系先造一份合成数据方便对照真实路径系数验证结果。这里生成两个潜变量一个叫满意度一个叫忠诚度真实路径系数设为 0.6满意度用 3 个显变量测量忠诚度用 2 个显变量测量。import numpy as np import pandas as pd np.random.seed(42) n 100 # 潜变量满意度 - 忠诚度路径系数 0.6 lv1 np.random.normal(0, 1, n) lv2 0.6 * lv1 np.random.normal(0, 0.8, n) # 显变量每个潜变量由多个观测指标加上噪声生成 data pd.DataFrame({ sat1: 0.8 * lv1 np.random.normal(0, 0.3, n), sat2: 0.7 * lv1 np.random.normal(0, 0.4, n), sat3: 0.75 * lv1 np.random.normal(0, 0.35, n), loy1: 0.7 * lv2 np.random.normal(0, 0.4, n), loy2: 0.65 * lv2 np.random.normal(0, 0.45, n), }) # 测量模型潜变量 - 观测指标 measurement { Satisfaction: [sat1, sat2, sat3], Loyalty: [loy1, loy2], } # 结构模型结果变量 - 原因变量 structural { Loyalty: [Satisfaction] }这段代码里的关键参数是噪声标准差。噪声越小显变量和潜变量的相关越高路径系数恢复越准。实际数据中指标载荷一般在 0.6-0.8 之间所以我把生成系数设成 0.7 左右模拟真实问卷的表现。structural字典的 key 是结果变量value 是原因变量列表注意不要写反否则路径方向就错了。3.2 核心迭代实现 PLS-PM 的权重更新循环下面这个类实现了前面讲的迭代四步。为了减小篇幅我把方法写紧凑但保留了关键注释。class PLSPM: 极简偏最小二乘路径建模实现。 measurement: 潜变量 - 显变量列名列表 structural: 潜变量(结果) - 潜变量(原因)列表 mode: 默认 A反映型也可传字典指定某个潜变量为 B n_iter: 最大迭代次数 tol: 收敛阈值 def __init__(self, measurement, structural, modeA, n_iter100, tol1e-5): self.measurement measurement self.structural structural self.mode mode if isinstance(mode, dict) else {lv: mode for lv in measurement} self.n_iter n_iter self.tol tol self.lvs list(measurement.keys()) self.weights {} self.scores {} self.path_coeffs_ {} self.loadings_ {} def _init_weights(self, data): for lv in self.lvs: self.weights[lv] np.ones(len(self.measurement[lv])) self.weights[lv] / np.linalg.norm(self.weights[lv]) def _outer_estimation(self, data): scores {} for lv in self.lvs: cols self.measurement[lv] score data[cols].values self.weights[lv] scores[lv] (score - score.mean()) / score.std() return scores def _neighbors(self, lv): neighbors set() for target, sources in self.structural.items(): if target lv: neighbors.update(sources) if lv in sources: neighbors.add(target) return list(neighbors) def _inner_estimation(self, scores): inner {} for lv in self.lvs: neighbors self._neighbors(lv) if not neighbors: inner[lv] scores[lv] continue composite 0.0 for nb in neighbors: corr np.corrcoef(scores[lv], scores[nb])[0, 1] composite corr * scores[nb] inner[lv] (composite - composite.mean()) / composite.std() return inner def _update_weights(self, data, inner): for lv in self.lvs: X data[self.measurement[lv]].values z inner[lv] if self.mode[lv] A: new_w X.T z / len(z) else: new_w np.linalg.pinv(X.T X) (X.T z) self.weights[lv] new_w / np.linalg.norm(new_w) def fit(self, data): data data.copy() for col in data.columns: data[col] (data[col] - data[col].mean()) / data[col].std() self._init_weights(data) for i in range(self.n_iter): old_weights {lv: self.weights[lv].copy() for lv in self.lvs} scores self._outer_estimation(data) inner self._inner_estimation(scores) self._update_weights(data, inner) diff np.mean([np.linalg.norm(self.weights[lv] - old_weights[lv]) for lv in self.lvs]) if diff self.tol: break self.scores self._outer_estimation(data) # 方向对齐让每个潜变量的第一个显变量与得分正相关 for lv in self.lvs: first_col self.measurement[lv][0] if np.corrcoef(data[first_col].values, self.scores[lv])[0, 1] 0: self.weights[lv] -self.weights[lv] self.scores[lv] -self.scores[lv] # 路径系数标准化得分后结构方程截距为 0 for target, sources in self.structural.items(): X np.column_stack([self.scores[s] for s in sources]) y self.scores[target] beta np.linalg.pinv(X.T X) (X.T y) self.path_coeffs_[target] dict(zip(sources, beta)) # 载荷得分与显变量的相关系数 for lv in self.lvs: cols self.measurement[lv] self.loadings_[lv] [np.corrcoef(data[col].values, self.scores[lv])[0, 1] for col in cols] return self这个类里有两个参数容易被忽略。第一个是mode它决定了权重更新是走一元回归还是多元回归第二个是tol默认 1e-5数据量小的时候可能震荡可以调到 1e-4 或 1e-6。_inner_estimation里的相关性权重是因子方案如果你想用质心方案只需把corr替换成np.sign(corr)但质心方案在相关接近 0 时容易不收敛。3.3 跑通最小案例用合成数据验证路径系数把模型实例化并拟合model PLSPM(measurement, structural, modeA, n_iter100, tol1e-5) model.fit(data) print(路径系数, model.path_coeffs_) print(外部权重, model.weights) print(载荷, model.loadings_)预期输出里Loyalty对Satisfaction的路径系数应该在 0.6 上下因为真实值就是 0.6。载荷应该接近 0.7-0.8表示每个观测指标和潜变量的相关。外部权重没有再等权分配而是按迭代结果调整后的一组相对大小。这里的验证逻辑是数据由已知关系生成如果算法正确路径系数必然能还原到 0.6 附近。如果跑出来是 0.1 或 -0.2首要检查structural的方向有没有写反其次检查mode是否错误地设置成了 B。4. 参数设定与模型调优权重模式、收敛阈值和标准化这些细节很多人拿到代码第一反应是“直接默认参数跑”。默认参数在演示数据上没问题但换到真实问卷数据大概率会碰壁。下面三个细节是我在实际项目中调得最多的。4.1 Mode A 还是 Mode B反映型与构成型测量怎么选这是 PLS-PM 最绕的一个选择。统计上有建议用一致变量consistent PLS来判但业务上更常见的是直接问删掉一个指标概念本身变不变如果概念没变是反映型用 Mode A如果概念内涵变窄了是构成型用 Mode B。维度Mode A反映型Mode B构成型箭头方向潜变量 - 显变量显变量 - 潜变量权重公式一元回归系数多元回归系数典型场景满意度、信任度经济地位、数字化水平多重共线性问题不大可能使权重不稳定一个常见误区是“看 Cronbach‘s alpha 高就选 Mode A”。alpha 测的是内部一致性只能支持反映型但不能证明构建是对的。我通常做法是先按理论画出测量模型再跑两种模式看路径系数符号是否一致如果符号颠倒说明测量模型设定有问题。4.2 收敛阈值和最大迭代次数别让循环跑飞收敛阈值tol设 1e-5 在样本量 100 左右比较稳妥但显变量载荷很低时权重更新步长会变小可能看起来没收敛但实际已经稳定了。这时把阈值放松到 1e-4能少走一半迭代。反过来如果数据质量好1e-6 也不会有问题。最大迭代次数n_iter默认 100 基本够用但如果你的模型中存在一个潜变量只有一个显变量或者结构模型是个环可能会迭代到边界。建议在循环里打印每个 epoch 的权重变化量观察是持续下降还是来回震荡# 在 fit 里的 diff 计算后加一行调试输出 # print(i, diff)如果 diff 在某个值附近反复横跳说明权重在振荡常见原因是多个潜变量高度相关内估计让得分来回倒。这时可以把tol设成 1e-3接受一个次优解或者检查有没有指标被同时放进了多个潜变量。4.3 用 Bootstrap 评估路径系数的显著性给结论加个置信区间PLS-PM 本身不给出 p 值所以业内标准做法是用 Bootstrap 重采样出路径系数的分布。具体来说对原始样本做有放回抽样重跑模型收集同一路径的系数计算置信区间。下面是一个最小实现n_boot 100 path_samples [] for _ in range(n_boot): idx np.random.choice(len(data), len(data), replaceTrue) boot_data data.iloc[idx].reset_index(dropTrue) boot_model PLSPM(measurement, structural, modeA, n_iter100, tol1e-5) boot_model.fit(boot_data) path_samples.append(boot_model.path_coeffs_[Loyalty][Satisfaction]) path_samples np.array(path_samples) ci_low, ci_high np.percentile(path_samples, [2.5, 97.5]) print(f95% CI: [{ci_low:.3f}, {ci_high:.3f}])这个代码块里的关键参数是n_boot。业务报告建议至少 1000 次探索阶段 100 次够看趋势。置信区间用百分位法优点是简单缺点是偏差修正不足但作为常规验证足够。如果区间包含 0说明这个路径在统计上不显著建议不要写进结论。5. PLS-PM 避坑指南5 个让我翻过车的细节5.1 显变量标准化不一致导致权重尺度离谱现象跑出来的外部权重有几十几百路径系数却小得可怜甚至符号和理论相反。原因我在早期版本里只对显变量做了标准化但外估计算出来的潜变量得分没有标准化导致权重更新时回归系数被得分方差放大。得分方差越大权重数字越夸张。解决在每次外估计之后对得分强制标准化也就是(score - mean) / std。同时确保进入模型的显变量都做了 z-score不要一部分归一化到 [0,1]一部分用原始值。5.2 潜变量得分符号漂移路径系数差点变负现象同一份数据跑两次路径系数绝对值差不多但符号有时候正有时候负。原因PLS-PM 迭代不约束潜变量得分的方向。满意度得分可以整体取反忠诚度得分也可以整体取反只要二者相对关系不变路径系数就会出现正负翻倍。解决在迭代结束后加一个方向对齐步骤比如让每个潜变量第一个显变量与得分的相关为正。如果相关为负就把权重和得分同时取反。这样至少保证符号可复现也能和业务解释对上。5.3 构成型指标错用 Mode A直接算出个不可解释的路径现象某个潜变量由市场份额、门店数、营收增速三个完全不同的指标组成用 Mode A 跑路径系数超过 1而且 Bootstrap 置信区间极其宽。原因这三个指标之间只有弱相关甚至负相关反映型测量模型强假设它们都是同一概念的表现实际上它们应该线性组合成“规模”这个潜变量也就是构成型。解决把对应潜变量的 mode 设为 B。如果不想全部改成 B可以在实例化时传字典mode{Scale: B, Satisfaction: A}。改成 B 之后权重变成多元回归系数能显式处理指标间的竞争关系。5.4 数据有缺失值NaN 一路传染到崩溃现象模型在np.corrcoef报错或者路径系数全是 NaN但数据明明只缺了几个值。原因pandas 的mean()和std()会跳过 NaN但矩阵乘法不会一旦出现一个 NaN整个得分向量变成 NaN后面的相关和回归全是 NaN。解决在fit之前做缺失值处理。如果缺失比例低于 5%用均值填充如果高于 10%考虑删除对应的行或列。注意均值填充会压缩方差所以填充后要重新标准化。千万不要把含 NaN 的 DataFrame 直接丢进模型。5.5 只用一组路径系数下结论不带置信区间现象报告里写“满意度对忠诚度的路径系数为 0.62”但换了 Bootstrap 样本后路径系数在 0.1-0.9 之间横跳。原因PLS-PM 本质是偏最小二乘没有内生显式误差分布单次估计的路径系数只是点估计没有可信度信息。解决养成跑完点估计就跑 Bootstrap 的习惯。最少 100 次正式报告用 1000 次。如果置信区间太宽优先检查样本量不足或测量模型载荷太低。这个坑不是代码 bug是分析流程问题却是最容易让结论翻车的。6. 分类问卷题和多人群组两个让模型真正落地的技巧6.1 二分类显变量先分清有顺序还是没顺序问卷里经常混进“是否回购”“性别”“地区”这类分类变量。如果是二分类且只有 0/1直接纳入标准化即可因为标准化后只是改变尺度不影响方向。但如果是无序多分类比如地区分为华北、华东、华南不能直接编码成 1、2、3否则算法会把这几个数字之间的等距关系当成真实距离产生虚假相关。我一般做的处理是有顺序的类别如 1-5 分直接作为数值没有顺序的类别先转成哑变量再作为某个构成型潜变量的显变量。举个例子data[is_repeat] (data[purchase_history] yes).astype(int)然后把is_repeat放进 measurement 字典对应潜变量的指标列表。这样既保留了信息又不会让无意义的方向污染权重。6.2 多群组比较按人群拆开跑别让平均值骗了你做过一次零售项目全样本跑出来“满意度 - 忠诚度”路径系数是 0.45但拆开新客和老客后新客是 0.12老客是 0.71。合并回归的 0.45 对谁都指导不了。PLS-PM 的多群组比较并不复杂就是把样本按分组变量切开展分开拟合group_a data[data[segment] new] group_b data[data[segment] old] model_a PLSPM(measurement, structural, modeA) model_b PLSPM(measurement, structural, modeA) model_a.fit(group_a) model_b.fit(group_b) print(model_a.path_coeffs_[Loyalty][Satisfaction]) print(model_b.path_coeffs_[Loyalty][Satisfaction])这里关键的比较逻辑是先看两个路径系数的 Bootstrap 置信区间是否重叠如果完全不重叠说明群体间机制有显著差异如果区间重叠很多那 0.45 和 0.71 的差别可能只是抽样噪声。我习惯把两组路径系数和区间画在同一张坐标图上肉眼一看就知道该不该做分群策略。这个习惯帮我挡住过好几次“拿总模型给一线业务下指标”的冲动。希望帮到你。本文还有配套的精品资源点击获取
返回列表