如何找到每个系数的p值(显著性)?
lm = sklearn.linear_model.LinearRegression()
lm.fit(x,y)
如何找到每个系数的p值(显著性)?
lm = sklearn.linear_model.LinearRegression()
lm.fit(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中用于这种统计分析的统计模型。
P_value是f个统计值之一。如果你想要得到这个值,只需使用这几行代码:
import statsmodels.api as sm
from scipy import stats
diabetes = datasets.load_diabetes()
X = diabetes.data
y = diabetes.target
X2 = sm.add_constant(X)
est = sm.OLS(y, X2)
print(est.fit().f_pvalue)
你可以用pingouin来写一行字。线性回归函数(免责声明:我是Pingouin的创建者),它使用NumPy数组或Pandas DataFrame与单/多变量回归一起工作,例如:
import pingouin as pg
# Using a Pandas DataFrame `df`:
lm = pg.linear_regression(df[['x', 'z']], df['y'])
# Using a NumPy array:
lm = pg.linear_regression(X, y)
输出是一个数据框架,其中包含每个预测器的beta系数、标准误差、t值、p值和置信区间,以及拟合的R^2和调整后的R^2。
稍微了解一下线性回归理论,下面是我们计算系数估计器(随机变量)的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值完全相同。
在多变量回归的情况下,@JARH的答案可能有错误。 (我没有足够的声誉来评论。)
在下面一行:
p_values = [2 * (1-stats.t.cdf (np.abs(我),(len (newX) 1)))我在ts_b),
t值遵循degree len(newX)-1的卡方分布,而不是遵循degree len(newX)-len(newX.columns)-1的卡方分布。
所以这应该是:
p_values = [2 * (1-stats.t.cdf (np.abs(我),(len (newX) len (newX.columns) 1)))我在ts_b)
(详见OLS回归的t值)