前回は線形回帰モデルの基本について見てきました。
特に、説明変数が一つの単回帰について見ましたが、今回は重回帰を行い、さらにもう少し統計的な結果の解釈をしていきたいと思います。
具体的には決定係数、自由度調整済み決定係数、モデルが有意かどうかを検定するF検定などについて見ていきたいと思います。
一部少し難しい話が入っていくるところもありますので、雰囲気だけわかれば良いという方は、そういったところは飛ばしていただいても大丈夫です。
また、皆様にも同じように手を動かして確認していただけると幸いです。
目次
statsmodels
前回はScikit-learnを使いましたが、今回はstatsmodelsという非常に便利なライブラリを紹介しつつ解説をしていきたいと思います。
statsmodelsも機械学習モデルを作成する上では非常に便利なツールですので、是非使い方と出力結果の意味を理解していただければと思います。
statsmodelsの公式ドキュメントはこちらです。
回帰分析をするにあたってScikit-learnとstatsmodelsのどちらを使うのが良いのかはなかなか難しいところですが、statsmodelは簡単に統計量を出してくれる点は非常に便利です。
またstatsmodelsは回帰分析だけでなく、ロジスティック回帰モデルなどの分類モデル、主成分分析、ARIMAモデルなどの時系列モデルなどなど様々なモデルを利用できます。
利用できるモデルは以下のドキュメントに載っていますのでご参考です。
重回帰モデル
前回はトヨタ自動車の売上高をトヨタ自動車の販売台数で回帰をしました。
今回は他にも売上高に影響を与えるであろうファクターのUSDJPY為替レートを取り入れてみたいと思います。
また、今回は売上の水準そのものではなく、前年度比何%増減するか？を予測するモデルにしたいと思います。
\(y\)はトヨタ自動車の売上の前年度比、\(x_1\)は自動車販売台数の前年度比、\(x_2\)はUSDJPY為替レートの前年度比として、以下のモデル式で表されます。
$$y=\beta_1 x_1 + \beta_2 x_2$$
\(\beta_1\)、\(\beta_2\)は最小二乗法で求める係数ですね。
まずはstatsmodels.OLSで単回帰
まず、statsmodelsとnumpyをインポートします。
import statsmodels.api as sm import numpy as np
USDJPY為替レートも用意し、データを作成します。
car_sales = np.array([8139584, 8423984, 8334485, 9692948, 10133095, 10167489, 10093781, 10250567, 10441126, 10602559, 10456617, 9919759]) sales = np.array([189509, 189936, 185836, 220641, 256919, 272345, 284031, 275971, 293795, 302256, 299299, 272145]) years = np.arange(2009, 2021) fx = np.array([92.46079167, 84.9741, 78.87079167, 83.30581667, 100.3575, 110.7883333, 120.0966667, 108.3425, 110.685, 111.0325, 108.3983333, 106.1825])
単回帰の結果と比較したいので、まずは単回帰用の自動車販売台数だけを説明変数にしたデータを作成します。
トヨタ自動車の自動車生産台数および売上高の前年比を計算し、それをそれぞれx_ret、y_retとします。
x = car_sales.reshape(len(car_sales), -1) x_ret = x[1:] / x[:-1] - 1 y_ret = sales[1:] / sales[:-1] - 1
では、ここから回帰分析をしていきますが、今回はstatsmodelsというライブラリのOLSというモジュールを使います。
OLSはOrdinary Least Squaresの略で最小二乗法のことです。
sm.OLSでデータをセットし、fitで回帰を実行しするだけです。
引数の順番がy、xなのでご注意ください。
model = sm.OLS(y_ret, x_ret) results = model.fit()
これで回帰分析が実行されましたが、その様々な結果がresultsという変数に入ります。
係数を見るにはresults.paramsを確認します。
print(results.params)
すると係数は1.31となっています。
つまり、自動車販売台数が1%増えるとトヨタ自動車の売上は1.31%増加することになります。
ここは概ね1に近くなっていればよいと考えられますが、ちょっと増えすぎという感じがしますね。
では次に予測結果が妥当かプロットをして確認しましょう。
まず予測を実行します。(model.predictではなく、results.predictであることにご注意ください)
prediction = results.predict(x_ret)
予測結果と実績を比較すると以下のようになります。
横軸が予測値で縦軸が実績値になっています。
45度の点線上に乗っていればいるほど予測が正しいということになりますが、ある程度うまく説明できているように見えますね。
もちろん数%の誤差はありますが、ある程度増減は当たっています。
（補足）statsmodelsで切片項を追加する
SklearnのLinearRegressionの場合は、fit_interceptにTrueまたはFalseを指定することで切片項を使用する/使用しないを設定しました。
statsmodelのOLSでは切片項が必要な場合は、add_constantで説明変数に1を追加します。
x_ret = sm.add_constant(x_ret)
ただし今回は、切片項は使わないのでこちらは必要ありません。
為替レートを追加した重回帰
では、トヨタ自動車の売上に影響すると考えられるUSDJPY為替レートを説明変数に追加してみましょう。
\({\bf{x}}\)に為替レートを加えて変化率を計算します。
x = car_sales.reshape(len(car_sales), -1) x = np.concatenate([x, fx.reshape(len(fx), -1)], axis=1) x_ret = x[1:] / x[:-1] - 1 y_ret = sales[1:] / sales[:-1] - 1
先ほどと同じように、最小二乗法でパラメータを求めて、結果を表示します。
model = sm.OLS(y_ret, x_ret) results = model.fit()
では係数を見てみましょう。
print(results.params)
\(\beta_1=1.060\)、\(\beta_2=0.517\)ということで、自動車販売台数が1%増加すると売上高は1.06%増加し、USDJPY為替レートが1%円安になると売上高は0.517%増加するという結果になっています。
直感には合ってますね。
(高級車かどうかなど)どういった車が売れるかによって売上の増加幅は違うと思いますが、平均的には1.06%となっており、おおむね1%に近くなっています。
同じように結果をプロットしてみましょう。
先ほどの単回帰の場合よりも点線(45度線)に乗っていることがわかりますね。
為替レートを追加することで説明力が上がり、より予測がうまくいっていることがわかります。
プロットして良くなっているということがわかりますが、この結果を数値を使って評価・確認していきましょう。
モデルの結果を確認・評価する
自動車販売台数に為替レートを追加することで予測が精緻になることがわかりました。
では結果の評価をしていきましょう。
resultsには様々な結果が入っており、例えば
print(results.resid)
とすると、残差(回帰で予測したあとの誤差)を表示してくれます。
モデルの結果をまとめてくれるsummary()
statsmodelsでは、results.summary()を呼び出すと結果をきれにサマリしてくれます。
print(results.summary())
このような結果が出ます。
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
==============================================================================
まず、上段の左側をサラッと見ると、ModelはOLS、MethodはLeast Squares(最小二乗法)になっています。
“No. Observation”はデータ数で11個になっています。
Df Modelはモデルの自由度を表し、ここではパラメータの個数です。(DfはDegree of Freedomを表します)
つまり、モデルがどれだけ自由に動かせる係数を持っているか、を表しています。
今回は切片項を使っていないので、パラメータは2つになります。
Df Residualsは残りの自由度なので、データ数からDf Modelを引いたものです。
データに比べDf Modelが大きい(=モデルが大きい)とDf Residualが小さくなり、極端な場合だとデータ数と説明変数の数が同じ場合はDf Residualがゼロになりますが、その場合は、方程式と変数の数が同じ連立方程式になっていますので、誤差を最小化するまでもなく、一意に係数が決まってしまいます(そうならない場合もありますがここでは割愛します)。
では、今回は上段右側の“R-squared”、“Adj. R-squared"、"F-statistic”、“Prob（F-statistic）”を少し詳しく見ていきたいと思います。
Log-Likelihood、AIC、BICはまた別の説明から必要になりますので、またの機会(次々回ぐらい？)にしたいと思います。
また、長くなってしまうので、係数の検定などの説明は次回にしたいと思います。
モデルの評価
決定係数(R-squared)
決定係数をざっくり言うと目的変数\(y\)の分散のうち説明変数\(X\)により説明できる分散の割合を表します。
決定係数は0から1の範囲を取り、1に近いほど予想値が実際のデータに近くなります。
決定係数が0の場合は、説明変数で説明できる割合がゼロということで、なにも説明力がないと考えます。
決定係数が1の場合は、説明変数で完全に目的変数の変動を説明できるということです。
また、相関係数の方が馴染みが多い人もいると思いますが、決定係数は相関係数の2乗になります。
先補dのトヨタ自動車の売上の結果だと、決定係数は、為替レートを入れない場合は0.706、為替レートを入れた場合は0.959となっています。
自動車販売台数だけだとトヨタ自動車の売上高変化率の分散の70.6%しか説明できなかったのですが、為替レートを追加することで全体の95.9%を説明できるようになりました。
為替レートを追加することによって大幅に改善していますね。
決定係数の計算方法
ここから決定係数はどうやって計算されているのかを説明しますが、若干話が難しくなるので、次のところまで読み飛ばしていただいても大丈夫です。
まず、目的変数\(y\)の分散を以下の式で求めます。
$$\sigma_y^2=\frac{1}{n}\sum^n_{i=1}\left(y^{(i)}-\bar{y}\right)^2$$
\(\bar{y}\)は\(y\)の平均を表しています。
先ほどの例だと、トヨタ自動車の売上高前年度比の分散に該当します。
これはトヨタ自動車の売上高が、平均に対してどれだけ散らばっているか？を表します。
そして、説明変数で条件付けした\(y\)の予測値を使った分散(残差の二乗和の平均)
$$\sigma_{y|x}^2=\frac{1}{n}\sum^n_{i=1}\left(y^{(i)}-\hat{y}_{x^{(i)}}\right)^2$$
を考えます。
\(\hat{y}\)は説明変数を使った\(y\)の予測値です。
先ほどの例で自動車販売台数のみのモデルを考えた場合では、自動車販売台数から予測されるトヨタ自動車の売上高から実際の売上高がどれぐらい散らばっているかを表します。
良いモデルであれば、単純に分散を取った\(\sigma_y^2\)よりも、予測を使った残差の分散\(\sigma_{y|x}^2\)の方が小さくなっていると考えられます。
そして、分散がどれだけ減ったかを表す
$$\sigma_r^2 = \sigma_y^2 - \sigma_{y|x}^2$$
を考えます。
これは無条件の分散\(\sigma_y^2\)が自動車販売台数を考慮したことによってどれだけ分散が減少するかを意味しています。
つまり、説明変数で説明できる分散になります。
そして決定係数はその比、
$$R^2=\frac{s_r^2}{s_y^2}=1-\frac{s_{y|x}^2}{s_y^2}$$
で表されます。
計算式からわかるように、決定係数は\(y\)の分散のうち\(x\)ｎよって説明される部分の割合を意味しています。
ちなみに、相関係数は1つの説明変数に対しては
$$\rho = \pm \sqrt{R^2}$$
で表されます。
回帰係数の符号がプラスであれば相関係数もプラス、回帰係数の符号がマイナスであれば相関係数もマイナスになります。
プラスであれば、説明変数の値が増加したときに目的変数の値も増加する、説明変数の値が減少したときに目的変数の値も減少し、これを正の相関や順相関と呼びます。
マイナスであれば、それぞれ逆に動きし、負の相関や逆相関と呼びます。。
なお説明変数が複数ある場合の相関係数は必ずプラスになります。
複数説明変数があると、増加する変数も減少する変数も混ざっているので、相関がプラスだとかマイナスだとかは言えないからです。
Statsmodelsの決定係数を確認
ではstatsmodelsの決定係数をpythonで計算して確認してみましょう。
注意点は、statsmodelsの決定係数はR-squared (uncentered)となっており、目的変数\(y\)の分散を求める際に平均を引いていません。
ですので、statsmodelsでは以下のように計算しています。
$$\begin{align}
\sigma_y^2&=\frac{1}{n}\sum^n_{i=1}{y^{(i)}}^2\\
\sigma_{y|x}^2&=\frac{1}{n}\sum^n_{i=1}\left(y^{(i)}-\hat{y}_{x_i}\right)^2\\
R^2 &= 1-\frac{\sigma_{y|x}^2}{\sigma_y^2}
\end{align}$$
この式に従って、以下のコードで計算することが可能です。
mse_y = np.sum(y_ret ** 2) / len(y_ret)
mse_y_x = np.sum((y_ret - prediction) ** 2) / len(y_ret)
rsquared = 1 - mse_y_x / mse_y
print(f'{rsquared:0.3f}')
計算結果は0.959となり、statsmodelsの結果と一致していることが確認できます。
自由度調整済み決定係数(Adjusted R-squared)
説明変数が複数ある場合や、説明変数の数が違うモデルを比較するときには自由度調整済み決定係数を使うことが一般的です。
自由度調整済み決定係数は「説明変数の数を考慮に入れた決定係数」と言えます。
通常、説明変数を増やせば増やすほど、決定係数は大きくなり、説明変数とデータ数が同じ数になると、完全に説明できることになります(説明変数間の相関が1や-1でない限り)。
ただ、それで良いモデルか？と言われるとそんなことはありません。
そこで、説明変数の数が複数ある場合は以下の式で表される自由度調整済み決定係数を用いることが一般的です。
$$\bar{R}^2=1-\left(1-R^2\right)\frac{n-1}{n-k}=R^2-\frac{k-1}{n-k}\left(1-R^2\right)$$
\(n\)はデータの数、\(k\)はパラメータの数、\(R^2\)は通常の決定係数を表します。
決定係数と自由度調整済み決定係数の関係
一番右の式で、\(R^2\)から正の値を引いていますので、自由度調整済み決定係数は通常の決定係数\(R^2\)よりも小さくなります。
パラメータが1つ\(k=1\)のときは決定係数と自由度調整済み決定係数は同じになります。
説明変数の数\(k\)が多くなると\(\bar{R}^2\)が小さくなり、\(R^2\)と\(\bar{R}^2\)の差は大きくなります。
逆に説明変数の数よりもデータ数がはるかに大きいときは、\(R^2\)と\(\bar{R}^2\)の差はほとんどなくなります。
トヨタ自動車の売上の例だと、為替レートを追加した場合、自由度調整済み決定係数でも0.950もあります。
データは少ないものの非常にうまく予測できています(当たり前と言えば当たり前ですが)。
自由度調整済み決定係数の導出
ここからは話が若干難しくなるので、読み飛ばしていただいても大丈夫です。
先ほどの決定係数を計算する際に分散を求めていましたが、その分散を不偏分散と呼ばれるものにします。
\(\sigma_y^2\)を求める際に\(n\)で割っていましたが、\(n\)で割るのではなく、\(n-1\)で割ることで不偏分散と呼ばれるものになります(不偏って？不偏分散って？という方は最後にご紹介する統計の本をオススメします)。
$$\bar{\sigma}^2_y=\frac{1}{n-1}\sum^n_{i=1}\left(y^{(i)}- \bar{y}^{(i)}\right)^2$$
また、\(\sigma_{y|x}\)計算時も\(n\)で割っていましたが、これをモデルの自由度である「データ数-パラメータ数」で割ってやります。
$$\bar{\sigma}_{y|x}^2=\frac{1}{n-k}\sum^n_{i=1}\left(y^{(i)}-\hat{y}_{x_i}\right)^2$$
最終的に自由度調整済み決定係数は決定係数のときと同じように計算して以下で表されます。
$$\begin{align}
\bar{\sigma}_y^2&=\frac{1}{n-1}\sum^n_{i=1}{y^{(i)}}^2\\
\bar{\sigma}_{y|x}^2&=\frac{1}{n-k}\sum^n_{i=1}\left(y^{(i)}-\hat{y}_{x_i}\right)^2\\
\bar{R}^2 &= 1-\frac{\hat{\sigma}_{y|x}^2}{\bar{\sigma}_y^2}
\end{align}$$
また、決定係数と自由度調整済み決定係数を見比べると、これらのの関係は以下のように計算できます。
$$\bar{R}^2=1-\left(1-R^2\right)\frac{n-1}{n-k}=R^2-\frac{k-1}{n-k}\left(1-R^2\right)$$
statsmodelsの結果を確認
statsmodelsでは、目的変数の分散は平均を引かないものを使います。
そして、予測値からの残差の二乗和については、自由度で割ります。
$$\begin{align}
\bar{\sigma}_y^2&=\frac{1}{n}\sum^n_{i=1}\left(y^{(i)}\right)^2\\
\bar{\sigma}_{y|x}^2&=\frac{1}{n-k}\sum^n_{i=1}\left(y^{(i)}-\hat{y}_{x_i}\right)^2
\end{align}$$
こちらのコードで計算できます。
mse_y = np.sum(y_ret ** 2) / len(y_ret)
mse_y_x = np.sum((y_ret - prediction) ** 2) / (len(y_ret) - 2)
rsquared_adj = 1 - mse_y_x / mse_y
print(f'{rsquared_adj:0.3f}')
結果は0.950とstatsmodelsの結果と一致しました。
念のため、決定係数と自由度調整済み決定係数の関係を使って、求めてみましょう。
statsmodelsでは、割り算の箇所の分子が\(n-1\)ではなく\(n\)になります。
$$\bar{R}^2=1-\left(1-R^2\right)\frac{n}{n-k}$$
rsquared_adj = 1 - (1 - rsquared) * (len(y_ret) ) / (len(y_ret) - n_params)
print(f'{rsquared_adj:0.3f}')
結果は、0.950と一致しましたね。
F統計量(F-statistic)
最後にF検定と呼ばれる検定の結果です。
ここでは検定について詳しくは解説しませんので、最後にご紹介する本などをご参考にしていただければと思います。
F統計量とは「係数がすべてゼロである」という仮説(帰無仮説と呼びます)が成立するかどうかを検定するものです。
帰無仮説\(H_0\)：係数\(\beta\)がすべてゼロである
帰無仮説が成立しないという結果になることを棄却されると言いますが、そのときは「係数がゼロではないのでモデルは有意である」と主張することができます。
つまり、ざっくり言うと「係数がすべてゼロ」と仮定した場合に、実際に推定された係数は、その仮説に当てはまるか、当てはまらないかを確認するものです。
この検定では帰無仮説を棄却することで強い主張ができるので、基本的には帰無仮説を棄却したいという思いがあります。
さて、検定の方法ですが、単純に最小二乗法で回帰しただけでは検定はできないので、モデルについて以下のような仮定を置きます。
$$y=\beta_0+\beta_1x_1 + \cdots + \beta_px_p + \epsilon$$
ここで\(epsilon\)は平均ゼロ、分散\(\sigma^2\)の正規分布に従います。
つまり、\(y\)の平均は\(\hat{y}=\beta_0+\beta_1x_1 + \cdots + \beta_px_p\)で、その周りに分散が\(\sigma^2\)でデータが分布することになります。
こういった仮定を置くことで検定が可能になります。
上記の仮定のもと、F統計量は以下で表されます。
$$\begin{align}
F&=\frac{SE_{y}-SE_{y|x}}{k} \left/ \frac{SE_{y|x}}{n-k} \right.\\
SE_y &= \sum^n_{i=1}(y^{(i)})^2 \\
SE_{y|x} &= \sum^n_{i=1}\left(y^{(i)}-\hat{y}^{(i)}\right)^2 \\
\end{align}$$
このF統計量は自由度\(\nu_1\)がパラメータ数\(k\)、\(\nu_2\)がデータ数\(n\)-パラメータ数\(k\)のF分布に従います。
ですので、これが大きいほど、帰無仮説「すべての係数がゼロである」が棄却されることになり、この回帰モデルが有意であるということができます。
ただ、これで検定しようとすると、この統計量とF分布の分布関数を照らし合わせないといけません。
そこで次で解説するp値を使うことで確認がしやすくなります。
statsmodelsのF statisticsを確認
statsmodelsのF値を確認しておきましょう。
F統計量はこちらのコードで計算できます。
se_y = np.sum((y_ret) ** 2)
se_y_x = np.sum((y_ret - prediction) ** 2)
f_statistic = ((se_y - se_y_x) / n_params) / (se_y_x / (n- n_params))
print(f'{f_statistic:0.3f}')
104.631となりstatsmodelの結果と一致しますね。
F検定のp値
続いてp値です。
詳しいことは最後に乗せておく参考書などを参照していただきたいのですが、ざっくり言うと、p値とは帰無仮説が成立する(この場合すべての係数がゼロである)と仮定した場合に、実際の観測値(この場合は推定されたパラメータ)の取る値は、どれぐらいの確率で起こるか？を表します。
ですので、このp値が小さいと帰無仮説は棄却される、つまり係数がゼロではないと言うことができます。
もう少し正確に言うと、p値が1%だとすると、本来の係数がすべてゼロだとすると、このような係数になる確率はたかだか1%しかないと言えます。
そこで、有意水準と呼ばれるものを5%とかに設定すると、この帰無仮説は棄却される(無に帰す)、つまりモデルが有意水準5%で有意だと言えることになります。
一方でp値が大きい、例えば40%だとすると、本来の係数がすべてゼロでも40%ぐらいの確率でこういった係数になるよね、という考え方です。
そして、検定が棄却されなかった場合に言えるのは、係数が有意(すべてゼロでない)とは言えない、ということになります。
「係数がすべてゼロである」とまでは言えないことにご注意ください。
あくまで、係数はもしかしたらゼロかもしれないけど「このp値だと係数がゼロではないとは言えない」、ということであって、「係数がゼロである」とまでは言えません。
この検定だと帰無仮説が棄却されないと強い主張はできないということですね。
トヨタ自動車の売上の例では、p値が5.87e-0.7となっていて非常に小さい値を取っています。
有意水準を5%としても十分帰無仮説を棄却でき、モデルの係数が有意であることがわかります。
statsmodelsでF検定のp値を確認
では、statsmodelsの結果を確認しましょう。
上記の F 統計量がさきほどの自由度を持つF分布に従うので、以下で計算することができます。
print(f'p値:{stats.f.sf(f_statistic, n_params, n-n_params)}')
結果は5.87e-0.7となりstatsmodelsの結果と一致しました。
（補足）F 検定について
F 検定は一般的に分散の比較・検定に使うものです。
それを応用して、係数をすべてゼロにした場合の分散と予測値を使った分散を比較し、予測値を使った場合に分散が減少しているかどうかを検定するものです。
まとめ
今回は、重回帰モデルを使って回帰し、評価結果を確認しました。
また、決定係数、自由度調整済み決定係数、F 統計量について詳しく見ました。
まだ、実務でも重要視されている係数の t 検定は説明していませんので、次回に説明したいと思います。
最後に今回の記事で参考にした以下の2冊の本をご紹介しておきます。
まずは超オススメの非常にわかりやすい一冊です。
最近よく見かける超わかりやすくした解説本ではありませんが(読んでいませんが)、統計学をこれから学び始めるという人で、しっかりと理解したい人には非常によい本です。
今回もいざ自分で説明すると非常に難しかったので、この本をすごく参考にしています。
あらためて説明うまいなぁと思いました。
もう一冊はこちらです。
1冊目と比べて少し難しいですが、何度も読めばよくわかります。
大学院時代の先生のオススメでよく読んでいました。
あと、この記事で作成しているグラフはPlotlyというPythonのライブラリを使っています。
こちらもオススメですので以下の記事をチェックしてみてください！
では！！