如何找到每个系数的p值(显著性)?

lm = sklearn.linear_model.LinearRegression()
lm.fit(x,y)

当前回答

scikit-learn的线性回归不计算这个信息,但你可以很容易地扩展类来做:

from sklearn import linear_model
from scipy import stats
import numpy as np


class LinearRegression(linear_model.LinearRegression):
    """
    LinearRegression class after sklearn's, but calculate t-statistics
    and p-values for model coefficients (betas).
    Additional attributes available after .fit()
    are `t` and `p` which are of the shape (y.shape[1], X.shape[1])
    which is (n_features, n_coefs)
    This class sets the intercept to 0 by default, since usually we include it
    in X.
    """

    def __init__(self, *args, **kwargs):
        if not "fit_intercept" in kwargs:
            kwargs['fit_intercept'] = False
        super(LinearRegression, self)\
                .__init__(*args, **kwargs)

    def fit(self, X, y, n_jobs=1):
        self = super(LinearRegression, self).fit(X, y, n_jobs)

        sse = np.sum((self.predict(X) - y) ** 2, axis=0) / float(X.shape[0] - X.shape[1])
        se = np.array([
            np.sqrt(np.diagonal(sse[i] * np.linalg.inv(np.dot(X.T, X))))
                                                    for i in range(sse.shape[0])
                    ])

        self.t = self.coef_ / se
        self.p = 2 * (1 - stats.t.cdf(np.abs(self.t), y.shape[0] - X.shape[1]))
        return self

从这里偷来的。

您应该看一看Python中用于这种统计分析的统计模型。

其他回答

稍微了解一下线性回归理论,下面是我们计算系数估计器(随机变量)的p值所需的总结,以检查它们是否显著(通过拒绝相应的零假设):

现在,让我们用下面的代码段计算p值:

import numpy as np 
# generate some data 
np.random.seed(1)
n = 100
X = np.random.random((n,2))
beta = np.array([-1, 2])
noise = np.random.normal(loc=0, scale=2, size=n)
y = X@beta + noise

用scikit-learn从上面的公式计算p值:

# use scikit-learn's linear regression model to obtain the coefficient estimates
from sklearn.linear_model import LinearRegression
reg = LinearRegression().fit(X, y)
beta_hat = [reg.intercept_] + reg.coef_.tolist()
beta_hat
# [0.18444290873001834, -1.5879784718284842, 2.5252138207251904]

# compute the p-values
from scipy.stats import t
# add ones column
X1 = np.column_stack((np.ones(n), X))
# standard deviation of the noise.
sigma_hat = np.sqrt(np.sum(np.square(y - X1@beta_hat)) / (n - X1.shape[1]))
# estimate the covariance matrix for beta 
beta_cov = np.linalg.inv(X1.T@X1)
# the t-test statistic for each variable from the formula from above figure
t_vals = beta_hat / (sigma_hat * np.sqrt(np.diagonal(beta_cov)))
# compute 2-sided p-values.
p_vals = t.sf(np.abs(t_vals), n-X1.shape[1])*2 
t_vals
# array([ 0.37424023, -2.36373529,  3.57930174])
p_vals
# array([7.09042437e-01, 2.00854025e-02, 5.40073114e-04])

使用statmodels计算p值:

import statsmodels.api as sm
X1 = sm.add_constant(X)
model = sm.OLS(y, X2)
model = model.fit()
model.tvalues
# array([ 0.37424023, -2.36373529,  3.57930174])
# compute p-values
t.sf(np.abs(model.tvalues), n-X1.shape[1])*2 
# array([7.09042437e-01, 2.00854025e-02, 5.40073114e-04])  

model.summary()

从上面可以看出,两种情况下计算的p值完全相同。

获取p值的一个简单方法是使用statmodels回归:

import statsmodels.api as sm
mod = sm.OLS(Y,X)
fii = mod.fit()
p_values = fii.summary2().tables[1]['P>|t|']

你可以得到一系列你可以操作的p值(例如,通过计算每个p值来选择你想要保持的顺序):

你可以用scipy表示p值。此代码来自scipy文档。

>>> from scipy import stats >>>导入numpy为np x = np.random.random(10) y = np.random.random(10) >>>斜率,截距,r_value, p_value, std_err = stats. linreturn (x,y)

另外一个已经提出的选择是使用排列测试。将y的值洗牌后对模型进行N次拟合,计算拟合后模型的系数相对于原模型的系数值较大(单侧检验)或绝对值较大(双面检验)的比例。这些比例就是p值。

scikit-learn的线性回归不计算这个信息,但你可以很容易地扩展类来做:

from sklearn import linear_model
from scipy import stats
import numpy as np


class LinearRegression(linear_model.LinearRegression):
    """
    LinearRegression class after sklearn's, but calculate t-statistics
    and p-values for model coefficients (betas).
    Additional attributes available after .fit()
    are `t` and `p` which are of the shape (y.shape[1], X.shape[1])
    which is (n_features, n_coefs)
    This class sets the intercept to 0 by default, since usually we include it
    in X.
    """

    def __init__(self, *args, **kwargs):
        if not "fit_intercept" in kwargs:
            kwargs['fit_intercept'] = False
        super(LinearRegression, self)\
                .__init__(*args, **kwargs)

    def fit(self, X, y, n_jobs=1):
        self = super(LinearRegression, self).fit(X, y, n_jobs)

        sse = np.sum((self.predict(X) - y) ** 2, axis=0) / float(X.shape[0] - X.shape[1])
        se = np.array([
            np.sqrt(np.diagonal(sse[i] * np.linalg.inv(np.dot(X.T, X))))
                                                    for i in range(sse.shape[0])
                    ])

        self.t = self.coef_ / se
        self.p = 2 * (1 - stats.t.cdf(np.abs(self.t), y.shape[0] - X.shape[1]))
        return self

从这里偷来的。

您应该看一看Python中用于这种统计分析的统计模型。