如何在 Python 中计算学生化残差？

Posted 2023-03-12

【中文标题】如何在 Python 中计算学生化残差？

【英文标题】：How to compute Studentized Residuals in Python?

【发布时间】：2018-01-11 02:44:02

【问题描述】：

我已经尝试寻找这个问题的答案，但到目前为止我还没有找到任何答案。我使用 statsmodel 在均值估算数据集上实现了普通最小二乘回归模型。我可以访问 OLS 结果中的残差列表，但不能访问学生化残差。如何计算/获取学生化残差？我知道计算学生化残差的公式，但我不确定如何在 Python 中编写这个公式。

提前致谢。

更新：我找到了答案。我可以从 OLS reults 的 outlier_test() 函数中获得一个包含学生化残差的数据框。

【问题讨论】：

欢迎来到 Stack Overflow！看起来您希望我们为您编写一些代码。虽然许多用户愿意为陷入困境的程序员编写代码，但他们通常只会在发布者已经尝试自己解决问题时提供帮助。展示这项工作的一个好方法是包含您迄今为止编写的代码、示例输入（如果有的话）、预期输出和您实际获得的输出（控制台输出、堆栈跟踪、编译器错误 - 不管是什么适用的）。您提供的详细信息越多，您可能收到的答案就越多。

查看 numpy 和 scipy 以及统计包是一个很好的起点

OLS 之后：Results 实例中有一个方法，更多信息请参见“另见”链接。 statsmodels.org/dev/generated/…

对于将来登陆这里的任何人，statsmodels 都会计算 Studentized Pearson Residuals。

【参考方案1】：

我正在处理同样的问题。解决方案是使用statsmodels 库：

from statsmodels.stats.outliers_influence import OLSInfluence

它包含一个resid_studentized_internal 方法。

【参考方案2】：

Nodar 的实现是不正确的，这里是来自 https://newonlinecourses.science.psu.edu/stat501/node/339/ 的更正公式以及删除的学生化残差，以防人们不想使用 statsmodels 包。两个公式返回的结果与上面链接中的示例相同

def internally_studentized_residual(X,Y):
    X = np.array(X, dtype=float)
    Y = np.array(Y, dtype=float)
    mean_X = np.mean(X)
    mean_Y = np.mean(Y)
    n = len(X)
    diff_mean_sqr = np.dot((X - mean_X), (X - mean_X))
    beta1 = np.dot((X - mean_X), (Y - mean_Y)) / diff_mean_sqr
    beta0 = mean_Y - beta1 * mean_X
    y_hat = beta0 + beta1 * X
    residuals = Y - y_hat
    h_ii = (X - mean_X) ** 2 / diff_mean_sqr + (1 / n)
    Var_e = math.sqrt(sum((Y - y_hat) ** 2)/(n-2))
    SE_regression = Var_e*((1-h_ii) ** 0.5)
    studentized_residuals = residuals/SE_regression
    return studentized_residuals

def deleted_studentized_residual(X,Y):
    #formula from https://newonlinecourses.science.psu.edu/stat501/node/401/
    r = internally_studentized_residual(X,Y)
    n = len(r)
    return [r_i*math.sqrt((n-2-1)/(n-2-r_i**2)) for r_i in r]

【参考方案3】：

使用OLSRresults.outlier_test() 函数生成包含每个观察的学生化残差的数据集。

例如：

#import necessary packages and functions
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

#create dataset
df = pd.DataFrame('rating': [90, 85, 82, 88, 94, 90, 76, 75, 87, 86],
                   'points': [25, 20, 14, 16, 27, 20, 12, 15, 14, 19])

#fit simple linear regression model
model = ols('rating ~ points', data=df).fit()

#calculate studentized residuals
stud_res = model.outlier_test()

#display studentized residuals
print(stud_res)

student_resid    unadj_p     bonf(p)
0   -0.486471   0.641494    1.000000
1   -0.491937   0.637814    1.000000
2    0.172006   0.868300    1.000000
3    1.287711   0.238781    1.000000
4    0.106923   0.917850    1.000000
5    0.748842   0.478355    1.000000
6   -0.968124   0.365234    1.000000
7   -2.409911   0.046780    0.467801
8    1.688046   0.135258    1.000000
9   -0.014163   0.989095    1.000000

本教程提供了完整的解释：https://www.statology.org/studentized-residuals-in-python/

【参考方案4】：

对于简单的线性回归，您可以使用以下方法计算学生化残差

将 X 和 Y 的平均值定义为：

mean_X = sum(X) / len(X) 
mean_Y = sum(Y) / len(Y) 

现在您必须估计系数 beta_0 和 beta_1

beta1 = sum([(X[i] - mean_X)*(Y[i] - mean_Y) for i in range(len(X))]) / sum([(X[i] - mean_X)**2 for i in range(len(X))]) 
beta0 = mean_Y - beta1 * mean_X

现在你需要找到拟合值，通过使用这个

y_hat = [beta0 + beta1*X[i] for i in range(len(X))]

现在计算残差，即 Y - Y_hat

residuals = [Y[i] - y_hat[i] for i in range(len(Y))]

我们需要找到H 矩阵，即 其中X 是我们的自变量矩阵。

要找到leverage，我们必须取H矩阵的对角元素，方法如下：

leverage = numpy.diagonal(H)

如果回归为找到标准误差

Var_e = sum([(Y[i] - y_hat[i])**2 for i in range(len(Y)) ]) / (len(Y) -2)
SE_regression = math.sqrt(Var_e*[(1-leverage[i]) for i in range len(leverage)])

现在您可以计算学生化残差

studentized_residuals = [residuals[i]/SE_regression for i in range(len(residuals))] 

请注意，我们有两种类型的学生化残差。一是Internally Studentized Residuals，二是Externally Studentized Residuals

我的解决方案找到内部学生化残差。

我在计算中进行了更正。对于外部学生化残差，请参阅@kkawabat 的答案

【讨论】：

我相信你的实现中有一个错误，因为这个结果与这里的玩具示例不匹配newonlinecourses.science.psu.edu/stat501/node/339

您缺少 SE_regression 的杠杆术语 sqrt(1-h_ii)

看来，我错过了计算中的杠杆项。谢谢@kkawabat 的更正。