top of page

回帰係数のp値とt値が不明なときに、p値とt値を求める方法:python 2022/05/02


さて、回帰係数のp値とt値が不明なときに、p値とt値を求める方法を説明していきます。この問題は回帰係数の標準偏差を求める方法と言い換えることが出来、それについてまとまった文章がネット上に見つからなかったので書いておきます。



0:便利な関数を用いて、回帰係数とp値とt値を求める方法

1:便利な関数を用いずに、回帰係数とp値とt値を求める方法

2:便利な関数を用いて回帰係数を求め、便利な関数を用いずにp値とt値を求める方法

3:2のやりかたでエラー(nan)が出た場合の対処方法



https://stackoverflow.com/questions/27928275/find-p-value-significance-in-scikit-learn-linearregression/50077476?answertab=scoredesc#tab-top


0と1については、上のサイトのJARHの答えとほぼ同じ内容です。


2と3がこのブログで追加したポイントです。偏最小二乗(PLS)回帰を使用される場合、p値とt値が出てこなかったりしますし、エラーを吐いたりするので、その場合は2や3のコードが使えると思います。


0:便利な関数を用いて、回帰係数とp値とt値を求める方法


まず、デフォルトのバージョンとして、線形回帰の例を見てみましょう。



import numpy as np
import statsmodels.api as sm

from scipy.stats import t  # We only need the t class from scipy.stats
from sklearn import datasets, linear_model
from sklearn.linear_model import LinearRegression
from scipy import stats

#default_version 0

diabetes = datasets.load_diabetes()
X = diabetes.data
y = diabetes.target

X2 = sm.add_constant(X)
est = sm.OLS(y, X2)
est2 = est.fit()
print("default_version 0")
print(est2.summary())


簡単なコードですが、Xの変数として切片を挿入しているところに注意する必要があります。

こうすると、以下の結果が出ます。



default_version 0
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.518
Model:                            OLS   Adj. R-squared:                  0.507
Method:                 Least Squares   F-statistic:                     46.27
Date:                Mon, 02 May 2022   Prob (F-statistic):           3.83e-62
Time:                        15:02:13   Log-Likelihood:                -2386.0
No. Observations:                 442   AIC:                             4794.
Df Residuals:                     431   BIC:                             4839.
Df Model:                          10                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t     P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       152.1335      2.576     59.061     0.000     147.071     157.196
x1          -10.0122     59.749     -0.168     0.867    -127.448     107.424
x2         -239.8191     61.222     -3.917     0.000    -360.151    -119.488
x3           519.8398     66.534     7.813      0.000     389.069    650.610
x4          324.3904     65.422      4.958     0.000     195.805     452.976
x5         -792.1842    416.684     -1.901     0.058   -1611.169      26.801
x6          476.7458    339.035      1.406     0.160    -189.621    1143.113
x7          101.0446    212.533      0.475     0.635    -316.685     518.774
x8          177.0642    161.476      1.097     0.273    -140.313     494.442
x9          751.2793    171.902      4.370     0.000     413.409    1089.150
x10          67.6254     65.984      1.025     0.306     -62.065     197.316
==============================================================================
Omnibus:                        1.506   Durbin-Watson:                   2.029
Prob(Omnibus):                  0.471   Jarque-Bera (JB):                1.404
Skew:                           0.017   Prob(JB):                        0.496
Kurtosis:                       2.726   Cond. No.                         227.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.


1:便利な関数を用いずに、回帰係数とp値とt値を求める方法


次に、回帰係数自体がそもそも算出出来ない場合は、回帰係数を求めるところから始めます。


# Calculate OLS regression coefficients: version 1
n = X2.shape[0]
k = X2.shape[1]
beta_hat = np.matmul(np.matmul(np.linalg.inv(np.matmul(np.array(X2).transpose(), np.array(X2))), X2.transpose()), y)
# Calculate Variance of OLS estimate
residual = y - np.matmul(X2, beta_hat)  # calculate the residual
sigma_hat = sum(residual ** 2) / (n - k - 1)  # estimate of error term variance
variance_beta_hat = sigma_hat * np.linalg.inv(np.matmul(X2.transpose(), X2))  # Calculate variance of OLS estimate

print("Calculate OLS regression coefficients: version 1")
print("beta",beta_hat)
t_stat = beta_hat / np.sqrt(variance_beta_hat.diagonal())
print("t_values",t_stat)

p_value = 1 - 2 * np.abs(0.5 - np.vectorize(t.cdf)(t_stat, n - k - 1))
print("p_values",p_value)

こうすると以下の結果が出ます。betaは回帰係数、t_valueはt値、p_valueはp値を示しています。

デフォルトの結果と今回の結果が概ね一致しているので、それなりに正しく計算出来ていると言えます。



Calculate OLS regression coefficients: version 1

beta [ 152.13348416  -10.01219782 -239.81908937  519.83978679 324.39042769
 -792.18416163 476.74583782  101.04457032  177.06417623 751.27932109
  67.62538639]
t_values [58.99287077 -0.16737594 -3.91263721  7.80412709  4.95267912 -1.89895643
 1.40455451  0.47487908  1.09526528 4.3653208   1.02368258]
p_values [0.00000000e+00 8.67152942e-01 1.06091213e-04 4.59632332e-14
 1.05375404e-06 5.82392749e-02 1.60875709e-01 6.35114256e-01
 2.74013518e-01 1.59181671e-05 3.06560773e-01]


2:便利な関数を用いて回帰係数を求め、便利な関数を用いずにp値とt値を求める方法


次に、便利な関数を用いて回帰係数は求められたが、p値とt値がわらかない場合のコードを示します。




#coefficients were known: version 2
model = LinearRegression()  #線形回帰モデルの呼び出し
model.fit(X, y)  #モデルの訓練

n = X2.shape[0]
k = X2.shape[1]

print("coefficients were known: version 2")

beta_hat2 = [model.intercept_] + list(model.coef_)
print("beta2",beta_hat2)
residual = y - np.matmul(X2, beta_hat2)  # calculate the residual
sigma_hat = sum(residual ** 2) / (n - k - 1)  # estimate of error term variance
variance_beta_hat2 = sigma_hat * np.linalg.inv(np.matmul(X2.transpose(), X2))  # Calculate variance of OLS estimate

t_stat = beta_hat2 / np.sqrt(variance_beta_hat2.diagonal())
print("t_values2",t_stat)

p_value = 1 - 2 * np.abs(0.5 - np.vectorize(t.cdf)(t_stat, n - k - 1))
print("p_values2",p_value)

このコードは、version1と基本的に一緒ですが、便利な関数を用いて見つけた回帰係数に切片が含まれているかどうかを確認することが大事です。この関数の場合は切片が別の名前(model.intercept_)で記録されていたので、それを用いましたが、もし切片が含まれていない場合は、X2 = sm.add_constant(X)のように、切片を含むようにして下さい。また、いくつかの関数ではmodel.coef_の中に切片が含まれている場合もあるので、その場合も注意が必要です。




coefficients were known: version 2
beta2 [152.1334841628965, -10.012197817470897, -239.81908936565566, 519.8397867901349, 324.3904276893769, -792.1841616283032, 476.7458378236605, 101.04457032134347, 177.0641762322499, 751.2793210873938, 67.62538639104415]
t_values2 [58.99287077 -0.16737594 -3.91263721  7.80412709  4.95267912 -1.89895643
 1.40455451  0.47487908  1.09526528 4.3653208   1.02368258]
p_values2 [0.00000000e+00 8.67152942e-01 1.06091213e-04 4.59632332e-14
 1.05375404e-06 5.82392749e-02 1.60875709e-01 6.35114256e-01
 2.74013518e-01 1.59181671e-05 3.06560773e-01]

このコードですと上のような結果が出ます。基本的に同じ結果になるので、正しく求められると言えそうです。


3:2のやりかたでエラー(nan)が出た場合の対処方法


さて、2のやり方で標準偏差が上手く求められないために、t値やp値が求められないことがあります。エラーの出る理由は逆行列のコードにあります。このコードでは、標準偏差を求める過程で逆行列を求めているのですが、その逆行列が存在しない場合にエラーが出されます。



variance_beta_hat2 = sigma_hat * np.linalg.inv(np.matmul(X2.transpose(), X2))

そのため、一般化(疑似)逆行列を求めるようにするために、上のコードを下のコードのように変更します。



variance_beta_hat2 = sigma_hat * np.linalg.pinv(np.matmul(X2.transpose(), X2))

少し分かり難いですが、np.linalg.invがnp.linalg.pinvに変更されています。





print("t_values were nan in version2: version 3")

beta_hat2 = [model.intercept_] + list(model.coef_)
print("beta2",beta_hat2)
residual = y - np.matmul(X2, beta_hat2)  # calculate the residual
sigma_hat = sum(residual ** 2) / (n - k - 1)  # estimate of error term variance
#variance_beta_hat2 = sigma_hat * np.linalg.inv(np.matmul(X2.transpose(), X2))  # Calculate variance of OLS estimate
variance_beta_hat2 = sigma_hat * np.linalg.pinv(np.matmul(X2.transpose(), X2))
t_stat = beta_hat2 / np.sqrt(variance_beta_hat2.diagonal())
print("t_values3",t_stat)

p_value = 1 - 2 * np.abs(0.5 - np.vectorize(t.cdf)(t_stat, n - k - 1))
print("p_values3",p_value)


上のようなコードを記入すると、下のような結果を得ます。



t_values were nan in version2: version 3
beta2 [152.1334841628965, -10.012197817470897, -239.81908936565566, 519.8397867901349, 324.3904276893769, -792.1841616283032, 476.7458378236605, 101.04457032134347, 177.0641762322499, 751.2793210873938, 67.62538639104415]
t_values3 [58.99287077 -0.16737594 -3.91263721  7.80412709  4.95267912 -1.89895643
 1.40455451  0.47487908  1.09526528 4.3653208   1.02368258]
p_values3 [0.00000000e+00 8.67152942e-01 1.06091213e-04 4.59632332e-14
 1.05375404e-06 5.82392749e-02 1.60875709e-01 6.35114256e-01
 2.74013518e-01 1.59181671e-05 3.06560773e-01]


今回は逆行列が存在するデータだったので、結果は変わりませんが、逆行列が存在しないデータに適用する場合は、標準偏差の値が変わりますので、t値やp値が変わります。


そのため、論文で使用する場合は、「本解析では逆行列ではなく、疑似逆行列を用いたため、ここでのt値やp値は推定値である」というように注釈を入れた方が丁寧です。



なお、このヴァージョン3を用いた解析は下記の論文の表2で使用しましたので、興味のある方はご覧下さい。


Yokotani, K. Spread of gambling abstinence through peers and comments in online self-help chat forums to quit gambling. Sci Rep 12, 3675 (2022). https://doi.org/10.1038/s41598-022-07714-2


以下にコードの全文を掲載しておきますので、適宜ご利用ください。


import numpy as np
import statsmodels.api as sm

from scipy.stats import t  # We only need the t class from scipy.stats
from sklearn import datasets, linear_model
from sklearn.linear_model import LinearRegression
from scipy import stats

#default_version 0

diabetes = datasets.load_diabetes()
X = diabetes.data
y = diabetes.target

X2 = sm.add_constant(X)
est = sm.OLS(y, X2)
est2 = est.fit()
print("default_version 0")
print(est2.summary())


# Calculate OLS regression coefficients: version 1
n = X2.shape[0]
k = X2.shape[1]
beta_hat = np.matmul(np.matmul(np.linalg.inv(np.matmul(np.array(X2).transpose(), np.array(X2))), X2.transpose()), y)
# Calculate Variance of OLS estimate
residual = y - np.matmul(X2, beta_hat)  # calculate the residual
sigma_hat = sum(residual ** 2) / (n - k - 1)  # estimate of error term variance
variance_beta_hat = sigma_hat * np.linalg.inv(np.matmul(X2.transpose(), X2))  # Calculate variance of OLS estimate

print("Calculate OLS regression coefficients: version 1")
print("beta",beta_hat)
t_stat = beta_hat / np.sqrt(variance_beta_hat.diagonal())
print("t_values",t_stat)

p_value = 1 - 2 * np.abs(0.5 - np.vectorize(t.cdf)(t_stat, n - k - 1))
print("p_values",p_value)



#coefficients were known: version 2
model = LinearRegression()  #線形回帰モデルの呼び出し
model.fit(X, y)  #モデルの訓練

n = X2.shape[0]
k = X2.shape[1]

print("coefficients were known: version 2")

beta_hat2 = [model.intercept_] + list(model.coef_)
print("beta2",beta_hat2)
residual = y - np.matmul(X2, beta_hat2)  # calculate the residual
sigma_hat = sum(residual ** 2) / (n - k - 1)  # estimate of error term variance
variance_beta_hat2 = sigma_hat * np.linalg.inv(np.matmul(X2.transpose(), X2))  # Calculate variance of OLS estimate

t_stat = beta_hat2 / np.sqrt(variance_beta_hat2.diagonal())
print("t_values2",t_stat)

p_value = 1 - 2 * np.abs(0.5 - np.vectorize(t.cdf)(t_stat, n - k - 1))
print("p_values2",p_value)

#t_values were nan in version2: version 3

model = LinearRegression()  #線形回帰モデルの呼び出し
model.fit(X, y)  #モデルの訓練

n = X2.shape[0]
k = X2.shape[1]

print("t_values were nan in version2: version 3")

beta_hat2 = [model.intercept_] + list(model.coef_)
print("beta2",beta_hat2)
residual = y - np.matmul(X2, beta_hat2)  # calculate the residual
sigma_hat = sum(residual ** 2) / (n - k - 1)  # estimate of error term variance
#variance_beta_hat2 = sigma_hat * np.linalg.inv(np.matmul(X2.transpose(), X2))  # Calculate variance of OLS estimate
variance_beta_hat2 = sigma_hat * np.linalg.pinv(np.matmul(X2.transpose(), X2))
t_stat = beta_hat2 / np.sqrt(variance_beta_hat2.diagonal())
print("t_values3",t_stat)

p_value = 1 - 2 * np.abs(0.5 - np.vectorize(t.cdf)(t_stat, n - k - 1))
print("p_values3",p_value)