简介这份资源面向数据建模与分析人员、科研工作者及结构方程模型学习者提供偏最小二乘路径建模PLS-PM算法的 Python 3 实现。该算法属于结构方程建模方法可基于潜在或表现变量估计复杂因果与预测模型特别适合探索性研究、中小样本和非正态分布数据可视为 R 语言 plspm 包在 Python 生态中的替代方案。压缩包共 63 个文件以 23 个 Python 源码模块为主体覆盖自助重抽样、参数估计、内部与外部模型、权重与量表计算等核心环节另有 33 个 CSV 示例数据集用于回归测试与演示辅以说明文档、配置文件及许可证整体仅 99KB轻量且目录规整。目前已有 2636 人学习下载。借助完整源码、测试套件和文档读者既能从底层理解各步骤的实现逻辑也能直接复用或改进这些模块将算法应用到自己的因果推断、客户满意度及社会科学建模项目中。1. 偏最小二乘路径建模的 Python 3 实现为什么你需要一个自己的轮子拿到一份 120 条样本的问卷数据变量基本都是 1-5 分的评分非正态、还有少量缺失。用协方差结构方程模型SEM拟合要么算不收敛要么报出负方差。这种场景下偏最小二乘路径建模PLS-PM往往是最务实的选择。这篇文章要讲的是一份不依赖 SmartPLS 或 SPSS 插件的 Python 3 实现——从迭代权重到路径系数全部用 numpy 手写。适合做用户满意度、服务质量、顾客忠诚这类小样本因果模型的从业者。读完你可以照着自己的数据结构改一改跑通一个最小案例并且知道哪些坑会让你白忙一场。2. 从测量模型到结构模型PLS-PM 的迭代原理与收敛条件2.1 反映式还是形成式测量模型选错整个模型就翻车PLS-PM 与传统的协方差 SEM 最大的不同在于它不要求显变量服从多元正态分布也不需要动辄几百的样本量。但天下没有免费的午餐它换来的是对测量模型假设的高度敏感。测量模型描述“潜变量如何通过显变量被测量”最常见的是反映式reflective潜变量是原因显变量是结果指标之间高度相关。另一种是形成式formative显变量是原因潜变量是结果指标之间不必相关。在 Python 实现里这两种模式对应不同的权重更新算法。反映式用模式 A权重等于显变量与潜变量得分的协方差形成式用模式 B权重等于显变量对潜变量得分做多元回归的回归系数。选错模式的后果很直接权重不收敛或者收敛了但载荷出现负值、符号与理论相反。常见判断标准是如果删掉一个指标其余指标含义依然能被潜变量覆盖那就是反映式如果指标各自代表潜变量某一独立维度缺一个都会影响潜变量定义那就是形成式。我在实际项目里处理过最典型的翻车案例用户用 5 个题目测“品牌形象”其中两个题目是“这个品牌很时尚”三个是“这个品牌很可靠”。从语义上看这是两个子维度本应拆成两个形成式潜变量但他硬塞进一个反映式潜变量迭代后其中一个题目的载荷变成 -0.1。这不是算法的问题是模型定义错了。2.2 重心法、因子法与路径法内部权重的三种选择内部权重决定了一个潜变量在结构模型中“参考邻居”的方式。三种常见策略重心法centroid与邻居潜变量得分的相关系数只取其符号即 1 或 -1计算简单迭代稳定但忽略了相关程度。因子法factor直接用皮尔逊相关系数作为权重把握了相关强度但容易放大噪声。路径法path区分前因和后果对前因变量用回归系数对后果变量用相关系数最符合路径建模直觉但计算量稍大。三种方法的选择并不玄学。样本量小、模型结构简单时重心法和因子法结果几乎一致模型复杂、路径数多时路径法通常更稳。下面这个表格可以直接抄进方案文档里内部权重法权重公式适用场景潜在风险重心法sign(corr(y_i, y_j))小样本、探索性分析丢失强度信息系数容易偏饱和因子法corr(y_i, y_j)指标质量较高对噪声敏感偶发不收敛路径法前因用回归 β后果用相关系数 r有明确因果方向实现复杂需维护方向表2.3 迭代收敛条件权重变化小于 1e-6 还是用 GOF 兜底PLS-PM 的求解本质是迭代固定点问题。每一轮的外层估计和内层估计相互交替更新权重向量直到新旧权重的最大绝对差小于阈值。我一般把收敛阈值设成 1e-6最大迭代次数设成 100。阈值太小容易在高维数据上耗尽迭代次数太大则得到的载荷和路径系数不稳定。阈值通常不需要低于 1e-7因为样本本身噪声远大于这个噪声量级。很多人忽略的一点是迭代收敛的是权重不是路径系数。权重稳定后潜变量得分才可信路径系数是用稳定得分回归出来的。如果迭代提前退出潜变量得分还处于半成品状态后续路径系数、R²、Bootstrap 置信区间全都会失真。所以实现时要同时记录迭代次数和权重最大变化量方便判断“收敛了”是真的收敛还是撞上了 max_iter 上限。3. 用 Python 3 从零写一个 PLS-PM 核心类3.1 模型定义与数据标准化直接喂 DataFrame 给迭代器先定义潜变量与显变量的对应关系以及结构模型中的路径。我习惯用字典和元组列表因为这两个结构能直接映射到 SmartPLS 里的“测量模型”和“路径图”。import numpy as np import pandas as pd # 潜变量 - 该潜变量对应的显变量列名 manifest_sets { SQ: [sq1, sq2, sq3], # 服务质量 CS: [cs1, cs2, cs3], # 顾客满意度 CL: [cl1, cl2, cl3] # 顾客忠诚度 } # 结构路径左 - 右 structural_relations [(SQ, CS), (CS, CL)] # 模拟 120 条观测3 个潜变量各 3 个显变量 rng np.random.default_rng(42) n 120 def make_latent(alpha): return rng.normal(0, 1, n) def add_noise(latent, reliab0.7): err rng.normal(0, np.sqrt(1 - reliab), n) return latent * np.sqrt(reliab) err潜变量得分需要做标准化但显变量进入迭代前也必须标准化。这里使用 z-score即减去均值除以标准差。为什么要标准化因为不同问卷题目的量纲可能不同比如一个题是 1-5 打分另一个题是 1-10 打分若直接进迭代权重会被量纲大的题目绑架。DataFrame 的列名必须和manifest_sets里的键值完全一致否则后面的矩阵切片会直接空指针。3.2 核心迭代外部估计与内部估计的交替更新下面这个PlsPm类是完整的最小实现。它只依赖 numpy不调任何第三方结构方程库。重点是_outer_estimate和_inner_estimate两个私有方法前者计算外部得分后者计算内部得分权重更新是两者的桥梁。class PlsPm: def __init__(self, manifest_sets, structural_relations, modeA, internal_weightscentroid, tol1e-6, max_iter100): self.manifest_sets manifest_sets self.structural_relations structural_relations self.mode mode # 测量模型模式A 反映式B 形成式 self.internal_weights internal_weights # 内部权重算法centroid/factor self.tol tol # 收敛阈值 self.max_iter max_iter # 最大迭代次数 self.all_cols [col for cols in manifest_sets.values() for col in cols] def _standardize(self, data): X data[self.all_cols].copy() return (X - X.mean(0)) / X.std(0, ddof1) def _neighbors(self, lv): nbs [] for a, b in self.structural_relations: if a lv: nbs.append(b) elif b lv: nbs.append(a) return nbs def _outer_estimate(self, X, weights): scores {} for lv, cols in self.manifest_sets.items(): scores[lv] X[cols].values weights[lv] scores[lv] (scores[lv] - scores[lv].mean()) / scores[lv].std(ddof1) return scores def _inner_estimate(self, scores): inner_scores {} for lv in self.manifest_sets: nbs self._neighbors(lv) inner np.zeros_like(scores[lv]) for nb in nbs: corr np.corrcoef(scores[lv], scores[nb])[0, 1] if self.internal_weights centroid: w np.sign(corr) elif self.internal_weights factor: w corr else: raise ValueError(internal_weights 只支持 centroid/factor) inner w * scores[nb] inner_scores[lv] (inner - inner.mean()) / inner.std(ddof1) return inner_scores def _update_weights(self, X, inner_scores): new_weights {} for lv, cols in self.manifest_sets.items(): X_lv X[cols].values if self.mode A: w X_lv.T inner_scores[lv] / len(X) elif self.mode B: w np.linalg.pinv(X_lv.T X_lv) X_lv.T inner_scores[lv] else: raise ValueError(mode 只支持 A 或 B) w w / np.linalg.norm(w) new_weights[lv] w return new_weights def fit(self, data): X self._standardize(data) self.X X # 等权初始化 weights {lv: np.ones(len(cols)) / np.sqrt(len(cols)) for lv, cols in self.manifest_sets.items()} it 0 for it in range(1, self.max_iter 1): old {lv: w.copy() for lv, w in weights.items()} scores self._outer_estimate(X, weights) inner_scores self._inner_estimate(scores) weights self._update_weights(X, inner_scores) diff max(np.max(np.abs(weights[lv] - old[lv])) for lv in weights) if diff self.tol: break self.weights_ weights self.scores_ self._outer_estimate(X, weights) self.iterations_ it self._compute_loadings() self._compute_paths() self._compute_r2() return self这段代码的逻辑顺序是标准化X→ 初始化等权权重 → 进入迭代循环 → 每次循环先用当前权重算外层得分再用外层得分算内层得分然后用内层得分更新权重。权重更新的核心是模式 A 的协方差公式X_lv.T inner_scores[lv] / len(X)即每个显变量与内层得分的协方差。权重最后做 L2 归一化是因为权重本身只代表相对比例绝对值不影响潜变量得分的标准化。收敛判定用的度量是最大绝对权重差。如果设置了tol1e-6通常 20 轮以内就能收敛。若超过 100 轮仍未达到说明模型定义或数据有严重问题应优先检查模式是否选对、是否存在缺失值。3.3 载荷、路径系数与 R² 的一次性计算收敛后需要输出三个关键结果载荷、路径系数、内生潜变量 R²。载荷是显变量与潜变量得分的相关系数路径系数是内生潜变量对它前因潜变量得分的回归系数R² 是回归模型解释的方差比例。def _compute_loadings(self): self.loadings_ {} for lv, cols in self.manifest_sets.items(): self.loadings_[lv] {} score self.scores_[lv] for col in cols: self.loadings_[lv][col] np.corrcoef(self.X[col], score)[0, 1] def _compute_paths(self): self.path_coefs_ {} endogenous set(b for _, b in self.structural_relations) for target in endogenous: predictors [a for a, b in self.structural_relations if b target] y self.scores_[target] X_p np.column_stack([self.scores_[p] for p in predictors] [np.ones(len(y))]) beta, _, _, _ np.linalg.lstsq(X_p, y, rcondNone) self.path_coefs_[target] { predictors: predictors, coefs: beta[:-1], intercept: beta[-1] } def _compute_r2(self): self.r_squared_ {} for target, info in self.path_coefs_.items(): pred_mat np.column_stack([self.scores_[p] for p in info[predictors]]) y_hat pred_mat info[coefs] info[intercept] y self.scores_[target] ss_res np.sum((y - y_hat) ** 2) ss_tot np.sum((y - y.mean()) ** 2) self.r_squared_[target] 1 - ss_res / ss_tot路径系数矩阵里我统一加了截距项因为潜变量得分虽然均值为 0但结构模型里前因变量未必完全中心化。lstsq方法返回的beta[:-1]是各自变量的系数beta[-1]是截距。R² 的计算采用最传统的一减残差平方和比总平方和。这段代码不需要额外传递参数它直接从self.scores_读取数据。4. 跑通一个最小案例服务质量、满意度与忠诚度的完整分析4.1 模拟一份有真实感的问卷数据继续用上一章的模拟函数生成三个潜变量SQ服务质量、CS顾客满意度、CL顾客忠诚度。设定 SQ 对 CS 的真实路径系数 0.6CS 对 CL 的真实路径系数 0.5每个潜变量的显变量载荷都在 0.7 附近。SQ_true make_latent(0) CS_true 0.6 * SQ_true rng.normal(0, 0.8, n) CL_true 0.5 * CS_true rng.normal(0, 0.8, n) data pd.DataFrame({ sq1: add_noise(SQ_true, 0.75), sq2: add_noise(SQ_true, 0.72), sq3: add_noise(SQ_true, 0.70), cs1: add_noise(CS_true, 0.74), cs2: add_noise(CS_true, 0.68), cs3: add_noise(CS_true, 0.71), cl1: add_noise(CL_true, 0.73), cl2: add_noise(CL_true, 0.69), cl3: add_noise(CL_true, 0.72), })这里add_noise的第二个参数是信度reliability约等于指标对潜变量的解释方差。0.7 是问卷研究里常用的心理测量阈值。信度设置太高会让结果显得过于完美太低则需要更大样本量才能收敛。模拟数据的价值在于你知道真实系数能反向检验实现是否正确。4.2 用 PlsPm 类拟合模型并检查迭代日志实例化并调用fit然后打印权重、载荷和迭代次数model PlsPm( manifest_sets, structural_relations, modeA, internal_weightscentroid, tol1e-6 ) model.fit(data) print(迭代次数:, model.iterations_) print(\n权重:) for lv, w in model.weights_.items(): print(lv, np.round(w, 4)) print(\n载荷:) for lv, loading in model.loadings_.items(): print(lv, {k: round(v, 4) for k, v in loading.items()})预期输出权重矩阵中每个潜变量下三个显变量权重大致均衡且方向一致。如果某个权重出现负号说明该指标与潜变量方向相反需要检查题干是不是反向计分。中心化后正负号本身没有绝对意义但与理论方向的匹配是模型有效性的第一步。4.3 结果解读路径系数和 R² 具体怎么看路径系数可以直接从字典里读取for target, info in model.path_coefs_.items(): for pred, coef in zip(info[predictors], info[coefs]): print(f{pred} → {target}: {coef:.4f}) print(f{target} R²: {model.r_squared_[target]:.4f})模拟数据里 SQ→CS 的路径系数应该约 0.55-0.65CS→CL 约 0.45-0.55。R² 紧接着给出解释力CS 被 SQ 解释约 35%-40%CL 被 CS 解释约 25%-30%。这些数值和真实参数基本吻合说明迭代实现没有系统性偏差。这里的路径系数是标准化系数可以直接比较不同路径的相对强弱。如果出现大于 1 的路径系数要警惕潜变量得分之间存在严重共线性或者内部权重算法选择不当。5. PLS-PM Python 实现避坑指南常见错误、翻车现象与排查方法5.1 权重不收敛或震荡最常见的杂症与解法现象diff长时间不降到 1e-6 以下迭代在 50 轮后仍在 0.001 量级震荡或者直接lstsq报奇异矩阵错误。原因最常见的是测量模型模式选错。一个实际是形成式的潜变量被定义成反映式模式 A权重更新公式会反复追逐不存在的相关性。另一个可能是指标间存在严重多重共线性pinv虽然能跑出结果但每次更新方向都不稳定。第三种是潜变量只有两个指标信息量不足权重在正负符号之间来回弹。解决先检查测量模型定义把疑似形成式的潜变量切到模式 B 试试。如果切换后收敛了说明定义确实有问题。如果仍然震荡用相关矩阵检查同一潜变量内部的指标相关性找出相关系数超过 0.9 的指标考虑合并。还可以把tol放宽到 1e-5但只适合探索阶段正式报告不建议低于 1e-5。5.2 潜变量得分符号方向反了别慌这是 PLS-PM 的身份问题现象其他结果都正常但某个潜变量下所有显变量的载荷都是负值且该潜变量到内生变量的路径系数也是负的。理论上这个关系应该是正相关。原因PLS-PM 迭代过程不固定潜变量得分的符号。由于权重向量有正负两个可行解迭代可能收敛到与理论方向相反的镜像解。这不算算法错误但严重影响解释。在问卷数据里常见于“满意度”这类正向潜变量如果某道题是反向计分符号翻转更容易发生。解决手动强制方向。在fit结束后检查理论路径系数如果符号与假设相反将该潜变量的权重向量全部取反再重新计算得分和路径系数。更稳妥的做法是在迭代前就固定锚点指标比如选定理论上的正向指标要求它的权重为正否则在每一轮更新后反转整个权重向量。这个动作要在权重标准化之前做否则会影响收敛速度。5.3 数据缺失与量纲问题为什么我的标准化没生效现象数据里有 NaNnp.corrcoef直接返回 NaN迭代循环崩溃。或者指标量纲差异极大权重被量纲大的指标主导载荷出现异常。原因PLS-PM 的标准迭代不处理缺失值X.std()对含 NaN 的列返回 NaN导致标准化结果全是 NaN。另一个隐蔽点是std(ddof1)与ddof0的选择样本量小于 30 时差异明显但实际影响不大更重要的是列名拼写错误导致data[cols]返回全 NaN 列。解决迭代前显式检查data.isnull().sum().sum()如果有缺失用均值插补或删除样本。删样本前要评估缺失是否随机的如果某道题缺失超过 10%建议先做多重插补不要直接硬删。量纲问题可以统一用StandardScaler再验证一遍确保均值为 0、方差为 1。我的血泪经验是不要相信apply(zscore)手动写(X - X.mean(0)) / X.std(0)更可控。5.4 模式 B 遇上多重共线性结构方程里的隐形炸弹现象使用模式 B 的潜变量其显变量权重数值异常大正负交错但绝对值都超过 0.8明显不合理。原因模式 B 的权重更新公式(X_lv.T X_lv) ^ -1 X_lv.T inner_scores本质上是对显变量做多元回归。当显变量之间相关系数高于 0.8矩阵求逆不稳定回归系数会被放大到荒谬地起“抵消”作用。解决模式 B 的显变量本应代表潜变量的不同维度所以先检查显变量之间的 VIF。如果 VIF 大于 5考虑删除冗余指标或者改用 PLS Regression 里的稀疏权重方法。实在不行就退回模式 A并在文档里说明测量模型是反映式形成式在共线性面前太脆弱。6. 进阶用 Bootstrap 验算路径系数的显著性不再被黑匣子迷惑6.1 写一个自带重抽样的 fit_bootstrap 方法验证路径系数是否显著最不容易被审稿人怼的方法是 Bootstrap。原理很简单从原始数据有放回抽样得到 B 组样本每组都重新跑一遍PlsPm.fit记录路径系数最后计算百分位数置信区间。若区间不包含 0就判定显著。def bootstrap_plspm(model_class, data, n_boot200, alpha0.95): boot_coefs {target: [] for target in set(b for _, b in model_class.structural_relations)} boot_r2s {target: [] for target in boot_coefs} indices np.arange(len(data)) for _ in range(n_boot): sample_idx rng.choice(indices, sizelen(data), replaceTrue) sample data.iloc[sample_idx] m PlsPm( model_class.manifest_sets, model_class.structural_relations, modemodel_class.mode, internal_weightsmodel_class.internal_weights, tolmodel_class.tol ).fit(sample) for target, info in m.path_coefs_.items(): for pred, c in zip(info[predictors], info[coefs]): boot_coefs[target].append({f{pred}-{target}: c}) for target in boot_r2s: boot_r2s[target].append(m.r_squared_[target]) return boot_coefs, boot_r2s注意每次重抽样后必须重新实例化类不能复用上一次的权重。因为 Bootstrap 样本间的权重初始值如果继承会让每一次拟合都从不同起点出发得到的结果差异会混淆抽样误差和迭代噪声。6.2 输出置信区间与显著性标记计算置信区间可以写一个简单函数def ci_from_boot(boot_values, alpha0.95): lower (1 - alpha) / 2 upper 1 - lower return np.percentile(boot_values, lower * 100), np.percentile(boot_values, upper * 100)一般 B200 是探索的最低标准正式分析建议 1000。每组 Bootstrap 都要重新迭代所以耗时线性增加。我一般先跑 100 次验证代码没写错再放 1000 次过夜。这个过程中最大的坑是每次重抽样本里会出现重复样本导致某些显变量的方差变成 0fit抛异常。所以实际代码里要捕获异常跳过那轮 Bootstrap而不是整体崩溃。6.3 我踩过的最大一个坑把显著性当效应量看Bootstrap 置信区间不包含 0 只代表“统计上显著”不代表路径系数大。在小样本里哪怕只有 0.1 的路径也可能显著因为 Bootstrap 的置信区间受样本量影响很大。我最早做案例时把 0.08 的值刻在里面画星星结果报告被导师一眼看穿说这根本没有实际意义。后来我习惯输出三条信息点估计值、置信区间、以及基于重抽样的效应量分布中位数。效应量至少达到 0.2 才值得在结论里大写特写。这个手写实现虽然简单但让我彻底摆脱了对黑匣子的恐惧。每次看到 SmartPLS 里那些默认参数我能立刻拆出它背后在跑什么。如果你也要自己实现记住先用模拟数据验证正确性再上真实数据最后再做 Bootstrap顺序不能乱。希望帮到你。本文还有配套的精品资源点击获取