回帰係数の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)
