Python与R三阶多项式回归系数差异问题排查求助

Ah, this is a classic gotcha between Python's PolynomialFeatures and R's poly() function! The core issue boils down to different default polynomial types being used in each language:

Key Difference in Feature Generation #

• Python's PolynomialFeatures(degree=3): Generates raw polynomial features: 1, x, x², x³ (unscaled, non-orthogonal, just the raw powers of your input variable).
• R's poly(zn, 3): By default generates orthogonal polynomial features, not the raw x, x², x³ terms. These orthogonal features are created via the Gram-Schmidt process to eliminate multicollinearity between polynomial terms, which changes the coefficient values entirely (even though the fitted prediction curve will be nearly identical).

Fix 1: Make R Use Raw Polynomials (Match Python) #

To get R to produce the same coefficients as your Python code, force it to use raw polynomial terms with the raw=TRUE argument in poly():

# Use raw polynomial features (matches Python's PolynomialFeatures behavior)
fit.zn2 <- lm(crim ~ poly(zn, 3, raw=TRUE))
summary(fit.zn2)

Running this should give you coefficients that align closely with your Python results: 4.8461; -0.3322; 0.0065; -3.776e-05.

Fix 2: Make Python Use Orthogonal Polynomials (Match R) #

If you want Python to replicate R's default orthogonal polynomial coefficients, you can manually orthogonalize the raw polynomial features using Gram-Schmidt QR decomposition (which is what R's poly() does under the hood):

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

def poly_three_orthogonal(n):
    x = np.array(boston[n]).reshape(-1, 1)
    # Generate raw polynomial features (without bias term)
    raw_features = PolynomialFeatures(degree=3, include_bias=False).fit_transform(x)
    # Apply Gram-Schmidt orthogonalization to match R's poly()
    orthogonal_features, _ = np.linalg.qr(raw_features)
    # Add intercept term (like R's lm() includes by default)
    data = sm.add_constant(orthogonal_features)
    # Fit OLS model
    model = sm.OLS(boston['CRIM'], data).fit()
    return model.summary()

# Call the function for 'ZN'
poly_three_orthogonal('ZN')

This will output coefficients matching R's default result: 3.6135; -38.7498 ; 23.9398 ; -10.0719.

Quick Verification: Predictions Are Consistent #

Even though the coefficients look wildly different, the actual prediction curves from both models will be nearly identical. Try predicting a value (e.g., zn=20) in both languages—you'll see the predicted crim value is the same (or extremely close, accounting for floating-point precision).

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