Python中基于QR分解处理共线特征的线性回归实现方法

现成可用方案 #

Python标准库没有自带和Rlm行为完全匹配的线性回归实现，不过第三方库statsmodels的OLS接口只要参数配置正确，可以直接复现你要的效果：存在完全共线性的冗余变量时，自动将对应系数、标准误、t值、p值填充为NaN，和你贴出的R输出逻辑完全对齐。

• 调用时需要手动添加截距项，标记矩阵包含常数项，同时指定用QR分解求解、设置和R一致的秩检测容差，模型会按输入列的先后顺序保留线性无关变量，把后出现的、可被其他变量线性表示的变量标记为未定义，变量选择规则和R完全一致。
• 要复现你贴的方差分析表，直接调用statsmodels的I型（序贯）平方和方差分析接口即可，R的anova.lm默认输出的就是I型平方和结果。

参考代码：

import numpy as np
import pandas as pd
import statsmodels.api as sm

# 假设dataset为pandas DataFrame，y列是因变量，其余列是自变量
X = dataset.drop(columns=["y"])
# 手动添加截距项，和R的y ~ . 公式逻辑一致
X = sm.add_constant(X)
# 初始化模型，指定QR分解求解
model = sm.OLS(endog=dataset["y"], exog=X, hasconst=True)
# 拟合时设置和R一致的秩检测容差，避免变量筛选结果不一致
fit_result = model.fit(method="qr", rcond=np.finfo(float).eps * max(X.shape))
# 输出和R格式类似的回归结果
print(fit_result.summary())
# 提取系数，共线性冗余变量会自动返回NaN
coefficients = fit_result.params
# 生成和R一致的I型方差分析表
anova_table = sm.stats.anova_lm(fit_result, typ=1)
print(anova_table)

手动实现方法（完全对齐R底层逻辑） #

如果不想依赖statsmodels，也可以手动复现Rlm的全部计算逻辑，核心步骤如下：

• 构建设计矩阵：在自变量矩阵最左侧拼接全1列作为截距项，和R的默认处理一致。
• 执行带列选主元的QR分解：R的lm底层调用LINPACK的dqrdc算法做选主元QR分解，会优先保留对因变量解释力强、输入顺序靠前的线性无关列，重排列顺序。
• 判定矩阵有效秩：用R默认的容差规则max(设计矩阵行数, 设计矩阵列数) * 双精度浮点数精度 * R矩阵对角线元素绝对值的最大值作为阈值，R矩阵对角线元素小于该阈值的列，就是存在完全共线性的冗余列。
• 求解系数：仅保留满秩部分的R子矩阵做回代求解非冗余变量的系数，冗余变量的所有统计量（系数、标准误、t值、p值）全部填充为NaN。
• 还原列顺序：将选主元过程中打乱的列顺序还原为原输入自变量的顺序，保证系数和输入变量一一对应。
• 计算后续统计量：残差、R方、调整R方、F统计量、方差分析表都仅基于保留的非冗余变量计算，结果和R的输出差异会在双精度浮点误差范围内。

注意：scikit-learn的LinearRegression不会返回NaN的核心原因是它默认用SVD求解，会对奇异值做软截断，给所有变量返回一个正则化后的数值解，不会显式标记共线性冗余变量，和Rlm的设计逻辑有本质区别。

结果对齐注意事项 #

• 秩检测的容差阈值必须和R保持一致，否则会出现保留/丢弃的变量数量、编号和R结果不匹配的问题，不要使用算法默认的宽松阈值。
• 如果要复现你贴出的方差分析结果，必须使用I型（序贯）平方和计算，不要用II型、III型平方和，R默认输出的方差分析表就是序贯检验结果。
• 截距项必须放在设计矩阵的第一列，否则选主元过程中可能误将截距项判定为冗余变量，导致结果和R不一致。

内容的提问来源于stack exchange，提问作者Ilyass Elamri