回帰分析のt値の求め方:Pythonで実装

統計学




t検定とt値とは

t値は、モデルの説明変数のそれぞれについて、計算されます。
このt値が「2」以上なら、その説明変数は、「統計学的に、モデルに組み込むのは良い」または、「統計学的に支持される」ことを意味する。
逆に、t値が「1」未満なら、「統計学的に、支持できない」ことになります。つまり、説明変数として、用いるのは宜しくない、ということを意味する。

偏回帰係数について、「説明変数の被説明変数への効果が0である」=「偏回帰係数が0」という帰無仮説を検定する。

両側5%の基準での検定
t値の絶対値が2未満
→帰無仮説が誤っているとは考えにくい

t値の絶対値が2以上
→帰無仮説が誤っていると考え、効果のある変数だと見なせる

$$
t値 = \frac{推定係数}{推定係数の標準誤差}
$$

回帰分析のt値の求め方

回帰分析におけるt値の求め方について紹介します。

単回帰分析(説明変数の数が1つ)の時と、重回帰分析(説明変数の数が2つ以上)の時では、t値の求め方が違うみたい??なので別々に紹介します。

計算式を見れば分かると思いますが、単回帰分析も重回帰分析もt値の求め方は同じです。

ただ、単回帰分析と重回帰分析で分けて、t値の求め方を紹介しているので、数式が少し違っています。

単回帰分析におけるt値の求め方

回帰線の傾きを検定する。


\(\hat{ \beta }\)は、平均が\(\beta\)で、分散が\(\frac{\sigma^2}{\sum{(x_i-\bar{x })^2}}\)の正規分布に従う。
\(\sigma^2\)は、誤差分散を表す。

この誤差分散は、未知なので残差平方和\(SS_e\)を使って標本から推定する。

$$
\begin{eqnarray}
\sigma^2 \rightarrow \hat{ \sigma^2 } &=& \frac{SS_e}{n-2} \\
&=& \frac{\sum{(y_i – \hat{y_i})^2}}{n-2}
\end{eqnarray}
$$

\(\hat{\sigma^2}\)を用いて、\(\hat{\beta}\)のt値を求める。

$$
\begin{eqnarray}
t &=& \frac{\hat{\beta}-\beta}{\sqrt{ \frac{\hat{\sigma^2}}{\sum{(x_i-\bar{x})^2}}}} \\
&=& \frac{\hat{\beta}-\beta}{\sqrt{ \frac{SS_e}{(n-2)\sum{(x_i-\bar{x})^2}}}}
\end{eqnarray}
$$

帰無仮説を\(\beta=0\)を代入する。

$$
t = \frac{\hat{\beta}}{\sqrt{\frac{SS_e}{(n-2)\sum{(x_i-\bar{x})}}}}
$$

$$
|t| > t_{n-2,\frac{\alpha}{2}}
$$
上の条件の時、帰無仮説を棄却する。

おおよそ有意水準が0.05の場合、\(2 \leq |t|\)ならば、帰無仮説を棄却する。

重回帰分析におけるt値の求め方

重回帰分析におけるt値の求め方は下の計算式です。

$$
\begin{eqnarray}
t値 &=& \frac{\beta_i}{\sqrt{SS^{ii} \times V_e}} \\
&=& \frac{\beta_i}{\sqrt{SS^{ii} \times \frac{SS_e}{n-p-1}}}
\end{eqnarray}
$$

それぞれの文字の説明

誤差分散の推定

\(SS_e\)は、残差平方和です。
$$
\begin{eqnarray}
SS_e &=& \sum{(\hat{y_i} -y_i)^2} \\
\end{eqnarray}
$$

\(\frac{SS_e}{n-p-1}\)は、誤差分散を表しています。
$$
\sigma^2 \rightarrow \hat{ \sigma^2 } = \frac{SS_e}{n-p-1}
$$

次のようにVeで表すこともある。
$$Ve = \frac{SS_e}{n-p-1}$$

nはデータの数、pは説明変数の数を表している。

n-p-1は、誤差の自由度を表す。

\(SS^{ij}\)

\(SS^{ij}\)は、偏差平方和・偏差積和の行列\(SS_{ij}\)の逆行列を表しています。

偏差平方和・偏差積和の行列\(SS_{ij}\)について、
説明変数の数は、pなので、\(SS_{ij}\)はp次正方行列になる。

$$
\begin{eqnarray}
SS_{ij} = \left(
\begin{array}{cccc}
SS_{ 11 } & SS_{ 12 } & \ldots & SS_{ 1p } \\
SS_{ 21 } & SS_{ 22 } & \ldots & SS_{ 2p } \\
\vdots & \vdots & \ddots & \vdots \\
SS_{ p1 } & SS_{ p2 } & \ldots & SS_{ pp }
\end{array}
\right)
\end{eqnarray}
$$

\(SS^{ij}\)は、\(SS_{ij}\)の逆行列を表すので、下のようになる。
$$
\begin{eqnarray}
SS^{ij} &=& \left(
\begin{array}{cccc}
SS_{ 11 } & SS_{ 12 } & \ldots & SS_{ 1p } \\
SS_{ 21 } & SS_{ 22 } & \ldots & SS_{ 2p } \\
\vdots & \vdots & \ddots & \vdots \\
SS_{ p1 } & SS_{ p2 } & \ldots & SS_{ pp }
\end{array}
\right)^{-1} \\
\\
&=&\left(
\begin{array}{cccc}
SS^{ 11 } & SS^{ 12 } & \ldots & SS^{ 1p } \\
SS^{ 21 } & SS^{ 22 } & \ldots & SS^{ 2p } \\
\vdots & \vdots & \ddots & \vdots \\
SS^{ p1 } & SS^{ p2 } & \ldots & SS^{ pp }
\end{array}
\right)
\end{eqnarray}
$$

偏差平方和・偏差積和は、次のような式で、表される。
偏差平方和は、分散にデータ数を掛けたもの
偏差積和は、共分散にデータ数を掛けたものである。

$$
SS_{ij} = \sum{(x_{ik} – \bar{x_i})(x_{jk} – \bar{x_j})}
$$

\(SS^{ii}\)は、対角成分を表す。

Pythonで重回帰分析のt値を求める

Pythonで、重回帰分析のt値を求めていきます。
今回は、numpyのみでプログラムしていきます。

今回は、Scikit-Learnにあるデータセット=ボストン住宅価格を使用します。

numpy

#import
from sklearn.datasets import load_boston
import numpy as np

#data_loading and data_setting
boston = load_boston()
data_X = boston.data
data_y = boston.target

#multiple regression
seg = np.ones((506, 1))
X = np.c_[seg, data_X]
X_t = X.T
beta_hat = np.dot(np.dot(np.linalg.inv(np.dot(X_t, X)), X_t), data_y)

#segmentation and coefficient
seg_value = beta_hat.T[0]
coef_value = beta_hat.T[1:]

#predict
y_hat = seg_value + np.dot(data_X, coef_value)

#calcurate t value
 
sse = (data_y - y_hat) **2
sse = np.sum(sse, axis=0)
sse = sse / (data_X.shape[0] - data_X.shape[1] -1 )

X_tmp = (X - X.mean(axis=0))
s = np.dot(X_tmp.T, X_tmp)
s = np.linalg.inv(s)

std_err = np.sqrt(np.diagonal(sse * s))

print("std_err:" ,std_err)

t = coef_value / std_err
print(t)

statsmodels

#chech t_value by statsmodels OLS
import statsmodels.api as sm
 
model = sm.OLS(data_y, sm.add_constant(data_X))
result = model.fit()
result.summary() 

sklearn-LinearRegression

lr = LinearRegression()
lr.fit(X, y)
coef_ = lr.coef_
y_hat = lr.predict(X)

sse =  np.sum((y - y_hat) **2, axis=0)
sse = sse / (X.shape[0] - X.shape[1] -1 )

s = np.linalg.inv(np.dot(X.T, X))

std_err = np.sqrt(np.diagonal(sse * s))

print("std_err : ", std_err)

t = coef_ / std_err

print("t : ", t)

参考

http://www.soumu.go.jp/ict_skill/pptx/ict_skill_3_4.pptx
http://lbm.ab.a.u-tokyo.ac.jp/~omori/kokusai/kokusai08_1218.html
https://www.kwansei.ac.jp/hs/z90010/sugakuc/toukei/rp12/rp12.htm

タイトルとURLをコピーしました