R与Python加权线性回归结果差异：Bug还是不同实现逻辑？

问题背景 #

学习加权线性回归时，对比R与Python（statsmodels库）的计算结果，发现二者在系数标准误、t值及F统计量上存在明显差异，以下是两个测试案例及核心疑问。

案例1：满秩矩阵的加权线性回归 #

数据集简单，系数估计值基本一致，但标准误、t值和F统计量存在差异。

Python代码及输出

import numpy as np
import statsmodels.api as sm
x = np.asarray ((18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29), dtype = float)
A = sm.add_constant (x)
y = np.asarray ((76.1, 77.0, 78.1, 78.2, 78.8, 79.7, 79.9, 81.1, 81.2, 81.8, 82.8, 83.5), dtype = float)
wt = np.asarray ((0.21, 0.15, 0.11, 0.0, -0.0, 0.11, 0.12, 0.15, 0.1, 0.2, 0.15, 0.4 ), dtype = float)
model = sm.WLS (y, A, weights = wt).fit ()
print (model.summary ())

输出：

WLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.992
Model:                            WLS   Adj. R-squared:                  0.991
Method:                 Least Squares   F-statistic:                     1241.
Date:                Fri, 19 May 2023   Prob (F-statistic):           8.07e-12
Time:                        11:55:00   Log-Likelihood:                   -inf
No. Observations:                  12   AIC:                               inf
Df Residuals:                      10   BIC:                               inf
Df Model:                           1
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         64.6927      0.456    141.786      0.000      63.676      65.709
x1             0.6452      0.018     35.223      0.000       0.604       0.686
==============================================================================
Omnibus:                        0.164   Durbin-Watson:                   2.140
Prob(Omnibus):                  0.921   Jarque-Bera (JB):                0.363
Skew:                           0.120   Prob(JB):                        0.834
Kurtosis:                       2.183   Cond. No.                         155.
==============================================================================

R代码及输出

A <- matrix(c( 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 ), nrow = 12, ncol = 1)
y <- c(76.1, 77.0, 78.1, 78.2, 78.8, 79.7, 79.9, 81.1, 81.2, 81.8, 82.8, 83.5)
wt <- c(0.21, 0.15, 0.11, 0.0, -0.0, 0.11, 0.12, 0.15, 0.1, 0.2, 0.15, 0.4 )
model <- lm (y ~ A, weights=wt)
summary (model)

输出：

Call:
lm(formula = y ~ A, weights = wt)

Weighted Residuals:
      Min        1Q    Median        3Q       Max 
-0.140456 -0.087469  0.007881  0.056603  0.166686 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 64.69272    0.51013   126.8 1.67e-14 ***
A            0.64524    0.02048    31.5 1.12e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1071 on 8 degrees of freedom
Multiple R-squared:  0.992,     Adjusted R-squared:  0.991 
F-statistic: 992.5 on 1 and 8 DF,  p-value: 1.121e-09

差异点：

• 单个系数的标准误和t值不同
• F统计量不同

案例2：秩亏矩阵的加权线性回归 #

R采用QR Householder方法处理秩亏矩阵，Python使用伪逆计算，导致重复变量的系数估计值不同，同时F统计量仍存在差异。

Python代码及输出

import numpy as np
import statsmodels.api as sm
x1 = np.asarray ([ -2, -1, 1, 2, 3, 4, 5 ], dtype=float)
x2 = np.asarray ([ 4, 1, 1, 4, 9, 16, 25 ], dtype=float)
x3 = np.asarray ([ 2, 4, 5, 7, 9, 12, 12 ], dtype=float)
y = np.asarray ([ 2, 3, 4, 5, 6, 7, 8 ], dtype=float)
A = np.asarray ([ np.ones (7, dtype=float), x1, x2, x2, x3 ]).T
wt = np.asarray ((0.11, 0.12, 0.0, 0.14, 0.15, 0.16, 0.15), dtype = float)
model = sm.WLS (y, A, weights = wt).fit ()
print (model.summary ())

输出：

WLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.998
Model:                            WLS   Adj. R-squared:                  0.996
Method:                 Least Squares   F-statistic:                     532.9
Date:                Fri, 19 May 2023   Prob (F-statistic):           0.000138
Time:                        11:53:14   Log-Likelihood:                   -inf
No. Observations:                   7   AIC:                               inf
Df Residuals:                       3   BIC:                               inf
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          3.0725      0.378      8.125      0.004       1.869       4.276
x1             0.6204      0.111      5.597      0.011       0.268       0.973
x2             0.0141      0.006      2.397      0.096      -0.005       0.033
x3             0.0141      0.006      2.397      0.096      -0.005       0.033
x4             0.0885      0.077      1.152      0.333      -0.156       0.333
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   2.656
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.452
Skew:                           0.542   Prob(JB):                        0.798
Kurtosis:                       2.389   Cond. No.                     1.46e+17
==============================================================================

R代码及输出

A <- matrix(c( -2, -1, 1, 2, 3, 4, 5,
              4, 1, 1, 4, 9, 16, 25,
              4, 1, 1, 4, 9, 16, 25,
              2, 4, 5, 7, 9, 12, 12), nrow = 7, ncol = 4)
y <- c(2, 3, 4, 5, 6, 7, 8)
wt <- c(0.11, 0.12, 0.0, 0.14, 0.15, 0.16, 0.15)
model <- lm (y ~ A, weights=wt)
summary (model)

输出：

Call:
lm(formula = y ~ A, weights = wt)

Weighted Residuals:
        1         2         3         4         5         6         7 
-0.040387  0.057390  0.000000 -0.017136  0.005993 -0.027315  0.022026 

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  3.07250    0.46316   6.634   0.0220 *
A1           0.62040    0.13577   4.570   0.0447 *
A2           0.02827    0.01445   1.957   0.1895  
A3                NA         NA      NA       NA  
A4           0.08849    0.09404   0.941   0.4461  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05695 on 2 degrees of freedom
Multiple R-squared:  0.9981,    Adjusted R-squared:  0.9953 
F-statistic: 355.2 on 3 and 2 DF,  p-value: 0.002808

核心疑问 #

问题1：案例1中R与Python的系数标准误和t值为何存在差异？ #

已分析R的标准误计算逻辑：

rse = np.sqrt(rss / residualDF)  # residualDF为8
standardError = rse * np.sqrt(np.diag(np.linalg.inv(A2.T @ A2)))

但不清楚Python的计算逻辑，是公式不同还是统计量定义有区别？

问题2：两个案例中R与Python的F统计量为何不同？哪种计算方式符合统计规范？ #

二者计算的RSS（残差平方和）和ESS（解释平方和）一致，但自由度使用不同：

• R的计算公式：(ess / DF_ess) / (rss / DF_rss)
• Python（statsmodels）的计算公式：self.mse_model/self.mse_resid

案例2的F统计量差异也源于类似的自由度处理差异，哪种方式是正确的统计规范？

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