前回はトヨタ自動車の売上高を目的変数(被説明変数)に、トヨタ自動車の自動車販売台数およびドル円為替レートを説明変数として、重回帰を行い、決定係数、自由度調整済み決定係数、F統計量とそのp値について見てきました。
今回は、個別の係数についてゼロでないかを検定する\(t\)検定などの結果を見て、線形回帰モデルの構築結果の解釈をしていきたいと思います。
目次
(前回のおさらい) モデル構築結果
まず、前回のおさらいですが、予測結果はこのようになっていました。
縦軸がトヨタ自動車の売上高の前年比、横軸が自動車販売台数とドル円の為替レートを使った予測値です。
当たり前と言えば当たり前ですが、それだけでかなり当てられているのがわかります。
そして、statsmodelsのサマリは以下のようになっていました。
OLS Regression Results ======================================================================================= Dep. Variable: y R-squared (uncentered): 0.959 Model: OLS Adj. R-squared (uncentered): 0.950 Method: Least Squares F-statistic: 104.6 Date: Mon, 07 Feb 2022 Prob (F-statistic): 5.87e-07 Time: 23:17:27 Log-Likelihood: 28.862 No. Observations: 11 AIC: -53.72 Df Residuals: 9 BIC: -52.93 Df Model: 2 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ x1 1.0596 0.111 9.541 0.000 0.808 1.311 x2 0.5174 0.070 7.431 0.000 0.360 0.675 ============================================================================== Omnibus: 1.252 Durbin-Watson: 1.764 Prob(Omnibus): 0.535 Jarque-Bera (JB): 0.392 Skew: -0.461 Prob(JB): 0.822 Kurtosis: 2.927 Cond. No. 1.77 ==============================================================================
今回はこの中段にある係数に関する項目”coef”、”std err”、”t”、”P>|t|”、”[0.025 0.975]”についてい見ていきたいと思います。
係数の検定とは
まず、係数の検定とは何のためにするのか?を説明します。
検定の必要性
例えば、このような母集団のデータがあったとします。
母集団全体を使って回帰した結果が上図の点線になります。
ただ、母集団全体を観測できることはあまりなく、実際の観測できるのはその一部です。
ですので、観測されるデータ、つまり母集団のうちの一部の標本を使ってパラメータを推定します。
例えば以下の図のように、水色のマーカーのデータのみが観測されている場合に、観測された標本のみを使って回帰をしたとしまs。
母集団全体を使った回帰結果とは乖離がありますね。
仮に違うデータが観測できていた場合は、また違った回帰結果になってきます。
例えば、以下のような結果になります。
同じ母集団ですが違うデータが観測されているため、その前の回帰結果とは違う結果になっています。
ここからわかるように、母集団は未知のため、母集団の係数\(\beta\)を推定しようとした場合、推定したパラメータには標本の取り方による誤差が含まれていることがわかります。
実際に、取得できるデータが少なく、目的変数と選んだ説明変数がたまたま関係性をもっているように係数が推定される場合もあります。
そこで、標本から得られた係数について、どのぐらいの誤差があるかを求め、その誤差を使って係数の妥当性を判断することが重要になってきます。
モデルの仮定
少し理論的な話になりますが、まず、係数の検定をする上で必要なモデルの仮定について説明します。
今までのように単にモデルが線形であるという仮定と最小二乗法で妥当な係数を推定する、というだけでは検定はできません。
検定を行うために、母集団が以下で表されるという仮定をおきます。
母集団が以下のような関係が成り立っている。
$$\begin{align}
y=\beta_0+\beta_1x_1 + \cdots + \beta_px_p + \epsilon, \hspace{10pt}\epsilon \sim N(0, \sigma^2)
\end{align}$$
また、各観測値\(y^{(i)}\)は独立である(\(y^{(i)}\)は\(i\)番目のサンプルを意味する)。
つまり、観測される\(y^{(i)}\)は平均が\(\beta_0+\beta_1x_1^{(i)} + \cdots + \beta_px_p^{(i)}\)でその周りに分散が\(\sigma^2\)で正規分布しているというものです。
そして、この係数\(\beta\)を推定したいという場合を考えます。
また、\(x\)は確率変数ではなく固定値とみなして、\(x\)は\(y\)の平均を決定するための変数と考えます。
\(x\)は確率変数ではなく固定値である。
このような仮定を置くことで、係数の検定が可能になります。
係数の評価
では、ここから係数の評価をしていきたいと思います。
まず、最小二乗法の簡単な復習です。
最小二乗法の復習
最小二乗法は以下の行列を使った式に対し、誤差の二乗を最小化するものでした。
$$\begin{align}
{\bf{y}}&=X\beta
\end{align}$$
\({\bf{y}}\)および\({\bf{\epsilon}}\)は\(N\times 1\)、\(X\)は\(N\times(p+1)\)、係数ベクトル\(\beta\)は\((p+1)\times 1\)です。
そして、最小二乗法で係数は次のように求まります。
$$\hat{\beta}=\left(X^TX\right)^{-1}X^T{\bf{y}}$$
\(\hat{\beta}\)は\(y^{(i)}\)の線形結合であり、各\(y^{(i)}\)が正規分布に従うという仮定を置いているので、\(\hat{\beta}\)も正規分布になります。
ここから、この正規分布の平均と分散について調べていきます。
推定された係数の平均
まず、係数の平均について見てみましょう。
係数の推定値\(\hat{\beta}\)の平均を求めてみます。
\(X\)は固定値で、\(\mathbb{E}[{\bf{y}}]=X\beta\)なので
$$\begin{align}
\mathbb{E}[\hat{\beta}]&=\left(X^TX\right)^{-1}X^T\mathbb{E}[{\bf{y}}] \\
&=\left(X^TX\right)^{-1}X^T\left(X\beta\right)\\
&=\beta
\end{align}$$
となります。
\(\hat{\beta}\)の平均は真の係数\(\beta\)になります。
この関係は最小二乗法で求めた係数\(\hat{\beta}\)は偏りのない推定値ということで不偏推定量と呼ばれます。
不偏推定量は式からわかるように、平均が真の値になるという良い性質を持っています。
(よく標準偏差の計算はデータ数\(n\)で割るのではなく、\(n-1\)で割った標本標準偏差を使っていますが、これは標本標準偏差が不変推定量だからです)
推定された係数の標準偏差(標準誤差)
続いて、標準偏差について見てみます。
まず共分散は以下のように計算できます。
(係数の部分を二乗して、\({\bf{y}}\)の分散を取ります)
$$\begin{align}
\text{Var}(\hat{\beta})&=\left(X^TX\right)^{-1}\left(\left(X^TX\right)^{-1}X^T\right)^T\text{Var}({\bf{y}}) \\
&=\left(X^TX\right)^{-1}X^T \sigma^2
\end{align}$$
ここで、\(\sigma^2=\text{Var}({\bf{y}})\)は\({\bf{y}}\)の分散を表します。
推定したパラメータの標準誤差(標準偏差)を知りたいのですが、母集団の分布は未知なので、上の式で出てくる真の分散\(\sigma^2\)がわかりません。
そこで、\(\sigma^2\)の代わりに\(\sigma^2\)の推定値\(\hat{\sigma}^2\)を求めて、それを使います。
ここでは二乗和をデータ数で割るのではなくデータの自由度で割って分散を求めます。
$$\hat{\sigma}^2=\frac{1}{N-p-1}\sum_{i=1}^N \left(y^{(i)}-\bar{y}\right)^2$$
\(\bar{y}\)は\({\bf{y}}\)の平均を表し、自由度はデータ数\(N\)からパラメータ数\(p+1\)を引いた\(N-p-1\)です。
自由度で割ることにより不偏推定量になります。
(パラメータ数は切片項を入れているので\(p+1\)になっています)
そして、上の式で\(\sigma^2\)を\(\hat{\sigma}\)で置き換えた、
$$\text{std_var}(\hat{\beta}) = \left(X^TX\right)^{-1}X^T \hat{\sigma}^2$$
がパラメータ\(\hat{\beta}\)の共分散の推定値になります。
\(\text{std_var}(\hat{\beta})\)は共分散なので\((p+1)\times(p+1)\)の行列になっています。
注目するのは各係数の分散である対角成分で、それらを取り出し各々平方根を取ると、各係数の標準誤差が求まります。
$$\text{std_err}_i = \sqrt{\text{std_var}_i}$$
\(\text{std_var}_i\)は\(\text{std_var}(\hat{\beta})\)の\(i\)番目の対角成分になります。
以上より、\(i\)番目の係数\(\beta_i\)の推定値\(\beta_i\)の標準偏差の推定値が\(\text{std_err}_i\)であるということがわかりました。
statsmodelsの結果を確認
では、statsmodelsの結果を確認しましょう。
先ほどの式どおりに関数を作成します。
def calc_std_err(x, y, prediction): var_hat_y = np.sum((y - prediction) ** 2) / (len(y) - x.shape[1]) var = np.linalg.inv(np.matmul(x.T, x)) * var_hat_y return np.sqrt(var.diagonal())
2行目は\(\hat{\sigma}^2\)を計算しています。
3行目で\(\text{std_var}(\hat{\beta})\)を計算し、4行目でその対角成分の平方根を計算しています。
こちらを呼び出して計算してみましょう。
std_err = calc_std_err(x_ret, y_ret, prediction) print(std_err)
自動車販売台数の係数の標準誤差が0.111と、為替レートに対する誤差が0.0696となりました。
statsmodelsとも一致していますね。
推定された係数には、それぐらいの誤差が含まれているということで、これが小さければ小さいほど推定精度も高いと考えられます。
\(t\)値
では、上記で求めた標準誤差を使って、各係数の推定値が妥当か検定します。
妥当かと言いましたが、\(t\)検定では一般的に、各係数がゼロでないかどうかを検定します。
ゼロと変わらないのであれば、その説明変数は不要だと考えることができます。
検定の方法は、推定した\(i\)番目の係数\(\beta_i\)の推定値が\(\beta_i\)で、観測データによる誤差が\(\text{std_err}_i\)という情報を使います。
その際に、以下のように、推定された係数から帰無仮説で仮定する係数を引いて、それを標準誤差で割った\(t\)値を使います。
$$t_i=\frac{\hat{\beta}_i-\beta_i}{\hat{\sigma}_i}$$
この統計量が自由度\(N-p-1\)の\(t\)分布に従うことがわかっています。
ちなみに、\(t\)分布は以下のような分布で、正規分布と似ています。
自由度が大きくなると正規分布に近づくという性質があります。
この分布を使って以下の帰無仮説を検定します。
帰無仮説: 係数\(\beta_i\)はゼロである。
そして、有意水準を1%や5%と決めて、この帰無仮説が棄却されるとこの説明変数を入れることに意味があるという結果になります。
一方で棄却できない場合(受容される場合)は、ゼロではないという強い主張はできません。
(個人的にはまったく当てはまっていない説明変数ならともかく、スモールデータで分析している場合で棄却できなくても、機械的に説明変数から落とすのではなく、一定の合理性があれば説明変数として採用しても良いと考えています)
statsmodelsの結果を確認
ここではいったん、\(t\)値を確認しておきましょう。
t_values = results.params / std_err print(t_values)
\(t\)値は\(\beta_1\)、\(\beta_2\)それぞれ9.541, 7.431とstatsmodelsと一致します。
\(t\)値の解釈
係数\(\beta_1\)の\(t\)値は9.541でした。
これをどう見るかというと、\(t\)値が9.541というのは帰無仮説「係数がゼロである」が成立するとした場合に、この\(t\)値が9.541となったということです。
そしてこの値は自由度\(n-p-1\)の\(t\)分布に従います。
この\(t\)値は\(t\)分布のどこの値でしょうか。
例えば、\(t\)値9.541の場合です。
すごく端っこの方ですね。
帰無仮説が成立していた場合に\(t\)値がこんな値になるような係数\(\beta\)が推定されたということですが、その確率は非常に小さいようです。
つまり、帰無仮説は成立していないと言えそうです。
これが真ん中の方、つまりその\(t\)値を取る確率が高い場所に位置していたら、「もしかしたら係数はゼロかもしれない」、もう少し正確には「係数はゼロではないとは言い切れない」という話になります。
これを基準を決めて定量的に見ていくのが次の\(p\)値です。
\(p\)値 (P>|t|)
さきほどの帰無仮説のもとで、\(t\)値が9.541や7.431という値を取る確率はすごく低いことがわかりますが、実際にその確率はどれぐらいでしょうか?
それを\(p\)値と呼びます。
自由度\(n-p-1\)の\(t\)分布で\(t\)値の値を取る確率が\(p\)値です。
1つ目の係数\(\beta_1\)の\(p\)値は5.284e-06、2つ目の係数\(\beta_2\)の\(p\)値は3.971e-05となっています。
非常に小さいですが、この意味は帰無仮説が成立していると、この\(t\)値を取る確率は高々5.284e-06、3.971e-05ぐらいの確率である、つまり、帰無仮説が成立している確率は高々0.0005284%と0.0003971%程度ということです。
例えば、これが0.2(=20%)であれば、帰無仮説が成立している確率が20%もあるので、自信をもって係数はゼロではない!とは言えません。
ただ、係数がゼロである、とまで言えないことにご注意ください。
この検定の枠組みでは帰無仮説が棄却されることが重要で、その場合に限って、自信をもって係数がゼロではないという強い主張ができます。
また、何%で棄却か?というのは自動的には決まらず、一般的には有意水準を1%や5%などと決めておき、\(p\)値が有意水準以下であれば有意水準〇%で帰無仮説は棄却される、などと言います。
有意水準自体は強い自信を持って言える必要があるかどうかなどを鑑みて決定することが多いのかと思います。
statsmodelsの結果を確認
まず、scipyというパッケージからtというモジュールをインポートします。
NumPyをベースとした統計、最適化、積分など多数の科学・工学モジュールを集めたオープンソース・ライブラリ。
これを使えば\(t\)分布の分布関数などの値を取得できます。
from scipy.stats import t
そして、以下のようにデータ数-パラメータ数を自由度として(この場合9)、\(t\)分布の累積分布関数を当てはめまます。
p_values = (1 - t.cdf(t_values, df=results.df_resid)) * 2 print(p_values)
ここで2倍しているのは、両側検定を使うためです。
イメージとしては、\(t\)値が9.5より大きくなる確率と同様に、-9.5より小さい値を取る確率を計算し、その合計を計算しています。
# 左側と右側の確率を合計する t.cdf(-t_values, df=results.df_resid) + (1 - t.cdf(t_values, df=results.df_resid))
これで計算された\(p\)値は一つ目の係数については5.284e-06、2つ目の係数については3.971e-05となります。
信頼区間
では最後に信頼区間についてです。
statsmodelsのsummaryでは2.5%と97.5%の信頼区間が出ていますが、これは95%(=2.5%+(100%-97.5%))の確率で真の係数\(\beta\)はこの範囲に入るだろうというものです。
一つ目の係数を見ると、2.5%点が0.808で97.5%点が1.311なので、95%の確率で係数は0.808から1.311の間にあるということを意味します。

さて、この信頼区間の求め方を見ていきます。
先ほど、\(t\)値が自由度\(n-p-1\)の\(t\)分布をするということを見ました。
$$t=\frac{\hat{\beta}-\beta}{\hat{\sigma}}\sim t(n-p-1)$$
\(t(n-p-1)\)は自由度が\(n-p-1\)の\(t\)分布の分布関数を表しています。
これを使って、\(t\)値が2.5%から97.5%の点に入る範囲は、
$$t_{0.025}(n-p-1)<t<t_{0.975}(n-p-1)$$
なので、
$$t_{0.025}(n-p-1)<\frac{\hat{\beta}-\beta}{\hat{\sigma}}<t_{0.975}(n-p-1)$$
となります。
これを変形すると、
$$\hat{\beta}+t_{0.025}(n-p-1)\hat{\sigma}<\beta<\hat{\beta}+t_{0.975}(n-p-1)\hat{\sigma}$$
と表されます。
つまり、真の係数\(\beta\)の2.5%の点\(\hat{\beta}+t_{0.025}(n-p-1)\)と97.5%の点\(\hat{\beta}+t_{0.975}(n-p-1)\)がこれで求まります。
statsmodelsの結果を確認
では、まずは2.5%タイル値を計算してみましょう。
print(results.params + t.ppf(0.025, df=9) * std_err)
t.ppfは\(t\)分布の逆関数を表しています。
1つ目の引数はパーセンタイル点、2つ目のdfは自由度です。
1つ目の係数\(\beta_1\)については0.808、2つ目の係数\(\beta_2\)については0.360となり、statsmodelsと一致しました。
97.5%タイル点についても計算して見ると、
print(results.params + t.ppf(0.975, df=9) * std_err)
\(\beta_1\)については1.311、\(\beta_2\)については0.675となり、こちらもstatsmodelsと一致していますね。
この幅が小さいほど、推定された係数の値が信頼できるということですね。
まとめ
今回は重回帰の係数に関する解釈・評価方法を見ていきました。
実務では、\(p\)値が5%以下の説明変数は除くといったこともされていると思います。
一方で、例えば〇〇をしたかどうか?を表す0, 1の変数で1となっているサンプルデータが少ない場合などでは、p値が少し大きく出ているけどこれは上式に照らし合わせて有効なはずだと考え、採用しているケースもありました。
個人的には、線形回帰を使って評価するような場合はデータが多くない場合も多く、より経験則をモデルに入れることが多いと思いますので、完全にデータだけで決めるのではなく、分析者・実務家の考えを入れるというのも良いと考えます。
最後に、もっと詳しく理解したいという方向けに、参考書を載せておきますので、参考にしてください。
改めて思いましたが、名著と呼ばれるような参考書は非常に厳密だし、体系的によくまとまっているので、非常に勉強になります。
このサイトでは、極力わかりやすく解説し、イメージを掴んでいただくことを目的としていますので、書き切れていないことも多いです。
ですので、このサイトでイメージを掴んでいただいたら、あとは参考書などでしっかりと学んでいただければ幸いです。
ということで、次はLassoやLidge回帰といった縮小推定、もしくはロジスティック回帰を見ていきたいと思います。
では!!
基本統計学
学問としてしっかりと解説されていますが、非常にわかりやすいです。
さらに平易に書かれている参考書もあると思いますが、この内容の濃さでは非常にわかりやすく、統計の基礎をしっかりと理解したい人にオススメです。
統計的学習の基礎
簡単ではないので、入門者にはあまりオススメしませんが、いつか統計・機械学習の専門家になりたい!という気合いの入っている人は持っていて良い一冊だと思います。
ボリュームが多く、個々の理論について非常に詳しく説明されています。
今回の記事で参考にしたので載せていますが、値段も高いのでご参考程度です(笑)。
ちなみに、この本では信頼区間の推定に\(t\)分布ではなく正規分布を使っていて、ちょっと混乱してしまいました…(そらデータが多ければ正規分布でいいんでしょうけど…)
Udemy講座
最後に、これからデータサイエンスを仕事にしたい!という人向けに、以下の記事でオススメのUdemy講座を紹介しています。
耳で聞いて理解する+本を読んで理解を深める、というのが効率の良い勉強法かと思いますのでご参考にしていただければと思います(+アウトプットがいいですね)。