时间戳

首页 > 解决方案 > 为 OLS 的预测计算学生化删除残差和 outlier_test()

python - 为 OLS 的预测计算学生化删除残差和 outlier_test()

问题描述

我正在使用 statmodel OLS 进行迭代异常值消除。我已经安装了模型使用。

ols_result = sm.OLS(y,X).fit()

然后我可以得到学生化删除残差外部和bonferroni

ols_result.outlier_test(method="bonf")

我正在删除超过 bonferroni p > %10 的样本，其中厨师距离在每次迭代中也是最高的。直到我没有得到 bonf(p)>%10 的样本，然后我得到了原始样本的子集。

假设我有 400 个样本，删除异常值后我得到了 380 个样本。现在我想用 400 个样本再次找到学生化删除残差和 bonferroni，以拟合 380 个样本的回归。查看删除的异常值是否真的是异常值。

这就是问题开始的地方。我一直在寻找简单的方法来使用 statmodels OLS 模型来获取拟合值的残差和厨师距离，而不是自己编写这些函数。但是.outlier_test()，.get_influence()似乎适用于 OLS Result 对象。

你们有什么简单的方法可以在没有太多代码等的情况下实现这些测试吗？

提前致谢

标签： pythonstatisticslinear-regression

解决方案

经过对这种实用性的深入研究。我找不到使用 statmodels 库的简单方法。相反，我编写了自己的代码，也从 statmodels 源代码中获得了帮助。如果有人需要我在下面分享的代码

def internally_studentized_residual(X,Y,y_hat):
    """

    Calculate studentized residuals internal

    Parameters
    ______________________________________________________________________

    X:            List
                  Variable x-axis

    Y:            List
                  Response variable of the X 


    y_hat:        List
                  Predictions of the response variable with the given X

    Returns
    ______________________________________________________________________

    Dataframe :   List
                  Studentized Residuals   


    """

#     print(len(Y))
    X = np.array(X, dtype=float)
    Y = np.array(Y, dtype=float)
    y_hat = np.array(y_hat,dtype=float)
    mean_X = np.mean(X)
    mean_Y = np.mean(Y)
    n = len(X)

#     print(X.shape,Y.shape,y_hat.shape)

    residuals = Y - y_hat

    X_inverse = np.linalg.pinv(X.reshape(-1,1))[0]
    h_ii = X_inverse.T * X

    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,y_hat):
    """

    Calculate studentized residuals external

    Parameters
    ______________________________________________________________________

    X:            List
                  Variable x-axis

    Y:            List
                  Response variable of the X 


    y_hat:        List
                  Predictions of the response variable with the given X

    Returns
    ______________________________________________________________________

    Dataframe :   List
                  Studentized Residuals External  


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


def outlier_test(X,Y,y_hat,alpha=0.1):
    """

    outlier test for the points

    Parameters
    ______________________________________________________________________

    X:            List
                  Variable x-axis

    Y:            List
                  Response variable of the X 


    y_hat:        List
                  Predictions of the response variable with the given X

    alpha:        float
                  alpha value for multiple test

    Returns
    ______________________________________________________________________

    Dataframe :   studentized, unadjusted p values and benferroni multiple 
                  test dataframe  


    """
    resid = deleted_studentized_residual(X,Y,y_hat)
    df = len(X) - 1 
    p_vals = stats.t.sf(np.abs(resid),df) * 2
    bonf_test = multipletests(p_vals,alpha,method="bonf")
    df_result = pd.DataFrame()
    df_result.loc[:,"student_resid"] = resid
    df_result.loc[:,"unadj_p"] = p_vals
    df_result.loc[:,"bonf(p)"] = bonf_test[1]
    df_result.index = X.index

    return df_result

这可能是一个有点脏的代码。我没有时间清洁并提高效率。它使用numpy模块和scipy.stats

推荐阅读