如何在Statsmodels中进行鲁棒回归（RLM）并得到R平方值？

Question

谈到拟合优度的衡量 - R-Squared 似乎是“简单”线性模型常用（并且被接受）的一种度量方式。 但是在 statsmodels （以及其他统计软件）中，RLM 并没有包括 R-squared 与回归结果。 是否有一种“手动”计算的方法，可以类似于 Stata 中的方法？ 或者是否有另一种度量方法，可以从 sm.RLS 生成的结果中使用 / 计算？

这就是 Statsmodels 产生的结果：

import numpy as np
import statsmodels.api as sm

# Sample Data with outliers
nsample = 50
x = np.linspace(0, 20, nsample)
x = sm.add_constant(x)
sig = 0.3
beta = [5, 0.5]
y_true = np.dot(x, beta)
y = y_true + sig * 1. * np.random.normal(size=nsample)
y[[39,41,43,45,48]] -= 5   # add some outliers (10% of nsample)

# Regression with Robust Linear Model
res = sm.RLM(y, x).fit()
print(res.summary())

输出结果为：

                    Robust linear Model Regression Results                    
==============================================================================
Dep. Variable:                      y   No. Observations:                   50
Model:                            RLM   Df Residuals:                       48
Method:                          IRLS   Df Model:                            1
Norm:                          HuberT                                         
Scale Est.:                       mad                                         
Cov Type:                          H1                                         
Date:                 Mo, 27 Jul 2015                                         
Time:                        10:00:00                                         
No. Iterations:                    17                                         
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          5.0254      0.091     55.017      0.000         4.846     5.204
x1             0.4845      0.008     61.555      0.000         0.469     0.500
==============================================================================

- Primer

---

由于RLM是通过迭代加权最小二乘法进行估计的，因此您可以尝试复制WLS实例wls_results = WLS(mod.endog, mod.exog, weights=mod.weights).fit()，其中mod是拟合后的RLM模型。但无法保证其准确性。WLS结果的R平方值具有加权残差的R平方值，这将是降低异常值影响的度量标准。但是，如果它们的权重不同，我认为您不能通过R平方比较模型。 - Josef

合适的答案在这里：https://github.com/statsmodels/statsmodels/pull/1341，其中包括基于 SAS 定义的 R方。 - Josef

谢谢，mod = sm.RLS(y, x); r2_wls = sm.WLS(mod.endog, mod.exog, weights=mod.fit().weights).fit().rsquared 帮了我，使得 R2=0.948。与 OLS 的 R2=0.731 相比，看起来有点“过于完美了” :-) - Primer

谢谢提供链接 - 在Github搜索类似问题时没有看到它。补丁中的函数生成R2 = 0.721。略低于OLS的R2...但是BIC从181下降到177（这是一个显著的变化吗？）？是否有其他措施来证明RLS清楚地并且以数字方式显示“最佳拟合”？ - Primer

我刚刚发现了这个链接：https://stat.ethz.ch/pipermail/r-help/2008-November/179773.html。首先，PR 1341 也修复了一些在 robust 中没有使用但在当前的 RLM 中需要的 bug。1341 中的 rsquared 是基于似然函数（或 M-估计目标函数）而不是基于残差平方和的伪 R 平方，并且 OLS 的 AIC 是基于正态分布的。虽然我有一段时间没有看过这个内容了，但是展示“更好的拟合”有点含糊，因为 RLM 会降低那些“不拟合”的观测值的权重，并将它们视为异常值。 - Josef

谢谢提供的链接 - 我从那个讨论中得出的结论是，可以使用“残差的比例估计...来确定预测精度”，而不是R2。看起来，在RLS结果中，比例是通过.scale属性进行访问的。然而，我找不到任何关于如何解释这个参数以及它实际上意味着什么的说明。在搜索过程中，我发现了一些可能感兴趣的论文：R平方的鲁棒版本，鲁棒AIC。 - Primer

Answer 1

由于OLS返回R2，我建议使用简单线性回归将实际值与拟合值进行回归。无论拟合值来自何处，这种方法都会为您提供相应的R2指示。

- majeed simaan

Answer 2

R2不是RLM模型拟合优度的好指标。问题在于异常值对R2值有巨大影响，以至于它完全由异常值决定。之后使用加权回归是一个有吸引力的替代方案，但最好查看估计系数的p值、标准误差和置信区间。

将OLS摘要与RLM进行比较（由于随机化不同，结果略有不同）：

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.726
Model:                            OLS   Adj. R-squared:                  0.721
Method:                 Least Squares   F-statistic:                     127.4
Date:                Wed, 03 Nov 2021   Prob (F-statistic):           4.15e-15
Time:                        09:33:40   Log-Likelihood:                -87.455
No. Observations:                  50   AIC:                             178.9
Df Residuals:                      48   BIC:                             182.7
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.7071      0.396     14.425      0.000       4.912       6.503
x1             0.3848      0.034     11.288      0.000       0.316       0.453
==============================================================================
Omnibus:                       23.499   Durbin-Watson:                   2.752
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               33.906
Skew:                          -1.649   Prob(JB):                     4.34e-08
Kurtosis:                       5.324   Cond. No.                         23.0
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

                    Robust linear Model Regression Results                    
==============================================================================
Dep. Variable:                      y   No. Observations:                   50
Model:                            RLM   Df Residuals:                       48
Method:                          IRLS   Df Model:                            1
Norm:                          HuberT                                         
Scale Est.:                       mad                                         
Cov Type:                          H1                                         
Date:                Wed, 03 Nov 2021                                         
Time:                        09:34:24                                         
No. Iterations:                    17                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.1857      0.111     46.590      0.000       4.968       5.404
x1             0.4790      0.010     49.947      0.000       0.460       0.498
==============================================================================

If the model instance has been used for another fit with different fit parameters, then the fit options might not be the correct ones anymore .

您可以看到，在OLS到RLM的转换中，无论是截距项还是斜率项，标准误差和置信区间的大小都有所降低。这表明估计值更接近真实值。

- Rob

Answer 3

为什么不使用model.predict来获取r2？例如：

r2=1. - np.sum(np.abs(model.predict(X) - y) **2) / np.sum(np.abs(y - np.mean(y)) ** 2)

- Tom Petrillo

---

这将被离群值所主导。 - Josef

@Josef - 通常，我会使用WLS机制，并在外样本数据上比较R2值（或研究特定指标）。是否有更好的机制？ - Tom Petrillo