OLS 适合带有系数误差和转换目标的 python

Posted 2023-03-12

【中文标题】OLS 适合带有系数误差和转换目标的 python

【英文标题】：OLS fit for python with coefficient error and transformed target

【发布时间】：2021-10-09 02:14:21

【问题描述】：

似乎有两种方法可以让 OLS 适合 Python。 Sklearn 一和 Statsmodel 一。我偏爱 statsmodel 之一，因为它通过 summary() 函数给出了系数的误差。但是，我想使用 sklearn 中的 TransformedTargetRegressor 来记录我的目标。似乎我需要在 statsmodel 中的拟合系数误差和能够在 statsmodel 中转换我的目标之间做出选择。有没有一种好方法可以在任一系统中同时完成这两项工作？

在统计模型中会这样完成

import statsmodels.api as sm
X = sm.add_constant(X)
ols = sm.OLS(y, X)
ols_result = ols.fit()
print(ols_result.summary())

返回拟合系数及其误差

对于 Sklearn，您可以使用 TransformedTargetRegressor

from sklearn.compose import TransformedTargetRegressor
from sklearn.linear_model import LinearRegression
regr = TransformedTargetRegressor(regressor=LinearRegression(),func=np.log1p, inverse_func=np.expm1)
regr.fit(X, y)
print('Coefficients: \n', regr.coef_)

但是如果不自己计算系数，就无法获得系数的误差。有没有两全其美的好方法？

编辑

我在这里找到了我关心的特殊情况的一个很好的例子

https://web.archive.org/web/20160322085813/http://www.ats.ucla.edu/stat/mult_pkg/faq/general/log_transformed_regression.htm

【问题讨论】：

***.com/help/minimal-reproducible-example

@PaulH 我尽我所能帮助更好地说明问题。这是一个功能问题，所以我不确定如何最好地说明它。

您从TransformedTargetRegression 获得了哪些您无法通过自己转换结果获得的功能？

@coffeinjunky ols_result.summary() 如果我手动转换，则以转换后的目标单位给出结果。这就是我得到系数错误的方式。我需要它们在非转换单元中。如果有一种方法可以转换它，甚至只是转换系数及其误差，那么这将是一个值得赏金的解决方案。

澄清一下，您对系数及其标准误差感兴趣，因为它们在对转换结果进行回归拟合后与未转换结果相关。您说您通常会从TransformedRegressor 获得此信息，但标准错误除外。到目前为止这是正确的吗？如果是这样，为了清楚起见，您能否添加您将用来从TransformedTargetRegressor 获得所需系数的确切命令？

【参考方案1】：

只是在这里添加一个冗长的评论，我相信TransformedTargetRegressor 不会做你认为它做的事情。据我所知，逆变换功能仅在调用predict 方法时应用。它不以未转换结果为单位表示系数。

示例：

import pandas as pd
import statsmodels.api as sm

from sklearn.compose import TransformedTargetRegressor
from sklearn.linear_model import LinearRegression
import numpy as np
from sklearn import datasets

创建一些示例数据：

df = pd.DataFrame(datasets.load_iris().data)
df.columns = datasets.load_iris().feature_names

X = df.loc[:,['sepal length (cm)', 'sepal width (cm)']]
y = df.loc[:, 'petal width (cm)']

Sklearn 第一：

regr = TransformedTargetRegressor(regressor=LinearRegression(),func=np.log1p, inverse_func=np.expm1)
regr.fit(X, y)

print(regr.regressor_.intercept_)
for coef in regr.regressor_.coef_:
    print(coef)
#-0.45867804195769357
# 0.3567583897503805
# -0.2962942997303887

关于转换结果的统计模型：

X = sm.add_constant(X)
ols_trans = sm.OLS(np.log1p(y), X).fit()
print(ols_trans.params)

#const               -0.458678
#sepal length (cm)    0.356758
#sepal width (cm)    -0.296294
#dtype: float64

您会看到，在这两种情况下，系数是相同的。也就是说，使用带有TransformedTargetRegressor 的回归产生与带有转换结果的statsmodels.OLS 相同的系数。 TransformedTargetRegressor 不会将系数反向转换为原始未转换空间。请注意，除非变换本身是线性的，否则系数在原始空间中将是非线性的，在这种情况下，这是微不足道的（与常数相加和相乘）。 This discussion 这里指向一个类似的方向 - 在大多数/许多情况下，反向转换 beta 是不可行的。

该怎么做？

如果解释是您的目标，我相信最接近您希望达到的目标是使用预测值，您可以改变回归量或系数。因此，让我举个例子：如果您的目标是说明sepal length 的一个标准误差较高对未转换结果的影响，您可以创建拟合的预测值以及 1- sigma 场景（通过系数回火，或通过 X 中的相应列回火）。

示例：

# Toy example to add one sigma to sepal length coefficient
coeffs = ols_trans.params.copy()
coeffs['sepal length (cm)'] +=  0.018 # this is one sigma


# function to predict and translate predictions back:
def get_predicted_backtransformed(coeffs, data, inv_func):
    return inv_func(data.dot(coeffs))

# get standard predicted values, backtransformed:
original = get_predicted_backtransformed(ols_trans.params, X, np.expm1)
# get counterfactual predicted values, backtransformed:
variant1 = get_predicted_backtransformed(coeffs, X, np.expm1)

然后你可以说例如关于未转换结果的均值偏移：

variant1.mean()-original.mean()
#0.2523083548367202

【讨论】：

对于这种情况，是的，您可以通过利用exp(sum of x_i) 是exp(x_i)s 的乘积的属性来做到这一点。本质上，exp(ln(y)) = exp(alpha + beta_1*x1 + beta_2*x2) = exp(alpha) * exp(beta_1*x1) * exp(beta_2*x2)。 .

对于连续变量，您需要为 x1 和 x2 选择适当的值（请参阅 Stata 的 GLM/LDV 模型的边际效应选项）。在你的情况下，这很容易，因为你有一个假人，所以它只会是exp(beta_2)。您应该通过将预测与所有的虚拟设置为 1 和所有的虚拟设置为 0 进行比较来获得相同的效果。在这两种情况下，我认为TransformedTargetRegressor 不会为你做任何事情，除非我忽略了一些事情。 @基思

只是为了非常清楚并确保我们在这里找到正确的树：如果您再次查看该公式，您会发现这都是乘法而不是加法。如果您仔细阅读 Stata 文章，他们会谈论几何平均值而不是算术平均值。你确定这是你最终想要的吗？如果不是，我仍然会使用预测值并简单地改变回归量的值。至于标准错误或 CI，我几乎可以肯定这并不那么简单。

你能描述一下你想如何使用标准错误吗？这是为了论文吗？转换后的模型中的推论是有效的，因此统计测试很好。 @基思

好的。因此，如果您的系数在统计上是显着的（使用适当的标准误差），那么您的推断就完成了。剩下的问题是它是否具有经济意义。在这里，上面描述的练习似乎给了你一个答案，不是吗？使用相关假人/交互 = 1 与交互 = 0 的预测值（以 $ 为单位）的差异。或者甚至熟悉几何平均故事并计算 exp(beta)。看起来您不需要这里的标准错误。

【参考方案2】：

简而言之，Scikit learn 无法帮助您计算系数标准误。但是，如果您选择使用它，您可以自己计算错误。在问题Python scikit learn Linear Model Parameter Standard Error@grisaitis 中提供了一个很好的答案，解释了它背后的主要概念。

如果您只想使用可与 sciait-learn 配合使用的即插即用功能，您可以使用：

def get_coef_std_errors(reg: 'sklearn.linear_model.LinearRegression',
                        y_true: 'np.ndarray', X: 'np.ndarray'):
    """Function that calculates the standard deviation of the coefficients of 
    a linear regression. 

    Parameters
    ----------
    reg : sklearn.linear_model.LinearRegression
        LinearRegression object which has been fitted 
    y_true : np.ndarray
        array containing the target variable
    X : np.ndarray
        array containing the features used in the regression

    Returns
    -------
    beta_std
        Standard deviation of the regression coefficients 
    """
    y_pred = reg.predict(X) # get predictions
    errors = y_true - y_pred # calculate residuals
    sigma_sq_hat = np.var(errors) # calculate residuals std error

    sigma_beta_hat = sigma_sq_hat * np.linalg.inv(X.T @ X)
    
    return np.sqrt(np.diagonal(sigma_beta_hat)) # diagonal to recover variances

【讨论】：

所以您的解决方案是使用 sklearn 并为错误编写代码，而不是使用 statsmodel 并手动转换？这就是我所倾向的方向。

是的，但这只是因为我更熟悉 Scikit-learn API

这似乎与@grisaitis 的方法不一样。我会用他的