R转Python多项式回归时变量显著性结果不一致问题咨询

问题描述 #

将R代码转换为Python时，多项式回归的显著性结果出现明显差异：

R代码（生成数据+拟合模型） #

set.seed(1)
x <- runif(n = 200, min = 0, max = 1)
y <- x^5 - 2 * x^4 + x^3 + rnorm(n = 200, sd=0.01)
df <- data.frame(x, y)

seventh_degree <- lm(formula=y~poly(x, degree=7), data=df)
summary(seventh_degree)

该代码输出的模型摘要显示1至5次多项式变量均显著。

Python原代码（复现流程） #

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.preprocessing import PolynomialFeatures

np.random.seed(1)
x = np.random.uniform(0, 1, 200)
y = x**5 - 2*x**4 + x**3 + np.random.normal(0, 0.01, 200)
df = pd.DataFrame({'x': x, 'y': y})

poly_features = PolynomialFeatures(degree=7,include_bias=False) 
x_poly = poly_features.fit_transform(df[['x']])  
x_poly = sm.add_constant(x_poly) 
lm_fit = sm.OLS(df['y'], x_poly).fit() 
print(lm_fit.summary())

但此代码输出的摘要显示无任何变量显著，与R的结果矛盾。

核心原因 #

两者的差异源于多项式特征的生成逻辑不同：

• R的poly()函数默认生成正交多项式，这类多项式的各阶项之间线性无关，彻底避免了多重共线性问题，回归系数的标准误不会被放大，因此能正确识别显著的项。
• Python的PolynomialFeatures默认生成原始幂次多项式（即$x, x^2, x^3,...,x7$），高次项与低次项之间高度相关，引发严重的多重共线性，导致回归系数的标准误急剧增大，最终表现为所有变量的显著性检验都不通过。

解决方法 #

在Python中生成正交多项式特征即可复现R的结果，有两种简便实现方式：

方式1：对原始多项式特征做正交化处理 #

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.preprocessing import PolynomialFeatures

# 生成模拟数据
np.random.seed(1)
x = np.random.uniform(0, 1, 200)
y = x**5 - 2*x**4 + x**3 + np.random.normal(0, 0.01, 200)
df = pd.DataFrame({'x': x, 'y': y})

# 生成原始多项式特征后做QR正交化
poly = PolynomialFeatures(degree=7, include_bias=False)
x_poly = poly.fit_transform(df[['x']])
# QR分解得到正交特征矩阵Q
Q, _ = np.linalg.qr(x_poly)
# 添加常数项
x_ortho = sm.add_constant(Q)

# 拟合模型并输出摘要
lm_fit = sm.OLS(df['y'], x_ortho).fit()
print(lm_fit.summary())

方式2：使用statsmodels公式接口直接生成正交多项式 #

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

# 生成模拟数据
np.random.seed(1)
x = np.random.uniform(0, 1, 200)
y = x**5 - 2*x**4 + x**3 + np.random.normal(0, 0.01, 200)
df = pd.DataFrame({'x': x, 'y': y})

# 用公式接口调用正交多项式拟合
lm_fit = smf.ols('y ~ poly(x, 7)', data=df).fit()
print(lm_fit.summary())

两种方法都能得到和R一致的结果：1至5次多项式变量显著，高次项不显著。

内容的提问来源于stack exchange，提问作者Boogie Ly