重回帰分析

単回帰分析はy = b0 + b1*xであり、求めたい係数はb0とb1だけであった。

また独立変数も1つだけである。

 

一方、重回帰分析はy = b0 + b1*x1 + b2*x2 + .... + bn*xn と求めたい係数も増え、また独立変数の数も二個以上である。

 

線形回帰で注意したいこと

以下の5つの仮定が成り立つ時にしか線形回帰ができない

 

・Linearity(線形性)

独立変数と従属変数の間に線形性がある

・homoscedasticity(等分散性)

独立変数と従属変数が同じ散らばり具合

・multivariate nomarity(多変量正規性)

多重線形回帰は、残差(従属変数yの観測値と予測値yˆの差が正規分布している。

・independence of errors(誤差の独立)

残差(従属変数yの観測値と予測値yˆの差は独立している。

・lack of multicollinearity(多重共線性の欠如)

独立変数間に強い相関があったり、一次従属な変数がある場合、解析が不可能であったり、たとえできたとしても信頼性が非常に低くなる。このような場合に多重共線性があるという。よって多重共線性がないことが大事。この仮定はVIFを使ってテストされる。

 

今回使うデータ

50社のスタートアップの資金投入や利益などのデータ(4社を抜き出している)

R&D Spend

Administration

Marketing Spend

State

Profit

165349.2

136897.8

471784.1

New York

192261.83

162597.7

151377.59

443898.53

California

191792.06

153441.51

101145.55

407934.54

Florida

191050.39

144372.41

118671.85

383199.62

New York

182901.99

 

R&D Spend:研究と開発にどれだけの費用を費やしたか

Administration:管理費

Marketing Spend:宣伝費

State:企業の所在地

Profit:利益

 

今回の従属変数はProfitであり、それ以外は独立変数である。

f:id:minmin_std:20191229173024p:plain

するとyがProfit、b1がR&D Spend、b2がAdmin、B3がMarketingであることは容易に想像できる。

 

しかしStateはどうだろうか??これはカテゴリー変数である。

そのためdummy variableを追加する必要がある。

 

StateにはNew YorkとCaliforniaの2つの変数がある。そのため1つの変数がわかればもう1つの変数も自動的にわかる。

 

具体的に新たな列を追加して、New Yorkの列が1ならその行の会社はNew Yorkにあるということになる。

 

そしてNew Yorkの列が1の時D1が1であり、0の時(Californiaの時)はD1が0になる。

つまり、係数が0になる。

 

しかし、この時Caloforniaの要素はb0に含まれている。

 

f:id:minmin_std:20191229174205p:plain

 

 

それでは、2つ目のdummy variable D2を加えたらどうなるか??

f:id:minmin_std:20191229183003p:plain

これは変数を複製していることになる。

どういうことか??

D2 = 1-D1という式が成り立つからである。

D2(California)が1の時、当然D1は0になり、逆もまた成り立つ。

 

これは上記の仮定のmulticollinearity(多重共線性)にあたり、正常に動作しなくなる。

これをdummy variable trap(ダミー変数トラップ)とよぶ。

 

 

p-valueとは??

Backward Elimination(変数減少法)に入る前に、p-valueについて、その仕組みを基本的に理解する必要がある。

https://haru-reha.com/p-value-mean/

このブログがわかりやすい。

 

 

モデルの作成

単回帰分析では1つの独立変数と1つの従属変数であったので簡単だった。しかし重回帰分析はいくつもの独立変数がある。

 

なのでそれらの中からどれを使い、どれを捨てるか決める必要がある。

 

なぜか??

 

理由は主に2つある。

 

1つ目は「garbage in ,gatbage out」である。これはコンピュータサイエンスにおける概念で、無意味なデータを入れると無意味なデータが返されるという意味である。つまり、無意味に独立変数を増やしすぎても、返されるデータは必ずしも良くならないというわけである。

 

2つ目は、ある変数が従属変数の振る舞いを予測しても、それを上司や仲間などに説明する必要がある(働いてたらそういうシチュエーションになるよねって話かな)。もし変数が1000もあったら、それを説明するのは大変。

 

これらに気をつけてモデルを作成する。

 

モデルを作る方法は5つある

1 : All-in

2 : Backward Elimination

3 : Forward Selection

4 : Bidirectional Elimination

5 : Score Comparison

 

Stepwise Regression(ステップワイズ回帰)はこの中の2、3、4である。が4のBidirectional Eliminationだけをステップワイズ回帰という人もいる。

 

 

1All-in

これはテクニカルな話ではなく、基本的に全ての変数を投げ入れることをいう。

これを使う時は、事前に知識がある場合である。

もしこれらの正確な変数が真の予測因子のものと知っていれば、事実であると既に知っているものを構築する必要はない。

専門知識から知っているかもしれないし、誰かの前にこのモデルを構築したので知っているかもしれない。

 

そして使わなければいけない時(仕事の都合などで)

 

最後に変数減少法(Backward Elimination)を用いて、準備する時

 

Backward Eliminationとは??

ステップ1:

モデルにとどまるための重要レベル(Sginificant Level)を選択する必要がある。したがって、デフォルトでは5パーセント、つまりSL(重要レベル)=0.05を使用し、次のステップで使用する。

 

ステップ2:

すべての可能な予測因子で完全なモデルを適合させるため、先ほど説明したアプローチですべてを実行し、すべての変数をモデルに入れて、それらを削除し始める。

 

ステップ3:

最も高いP値を持つ予測変数を検討する。

p>SLの時、ステップ4に進み、それ以外の時はモデルの準備が完了したことを意味する。

 

ステップ4:

基本的に最高のpを持つ変数を削除するためにその予測因子を削除する

 

 正直何を言っているかよくわからない。。。。。

 

が5つの方法の中で最も早く、また実践的なのがこれらしい。

 

 

pythonでモデルを作る

線形単回帰と全く同じである。

#Fitting Multiple Linear Regression to the Training set
from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(X_train,y_train)

#Predictiong the Testting set result
y_pred = regressor.predict(X_test)

 

実行結果

y_testの値は

103282
144259
146122
77798.8
191050
105008
81229.1
97483.6
110352
166188

y_predの値は

103015
132582
132448
71976.1
178537
116161
67851.7
98791.7
113969
167921

 

となり、しっかりと予測できていることがわかる。

 

とはいえ、なぜ線形単回帰と全く同じような実装で重回帰ができるか謎。、、

おそらく重回帰分析のやり方がLinearRegressionの中に入っていて、変数がたくさんある時は、重回帰分析になって1つの時は線形単回帰分析になるのではないかという予測。

まあ、中身で何をしているのかよくわかっていないけど。。

 

 

ここでしたやり方は全ての独立変数を使ったやり方であり、どの独立変数が予測の上で効果的かとか効果的でないのかとかは考えてない。

 

つまりAll-inのやり方だと思う。

 

backward eliminationを使って最適なモデルを作る

次はbackward eliminationを使って最適なモデルを作っていく。

 

最初にライブラリをインポートする

import statsmodels.formula.api as sm

 

次にbackward eliminationがちゃんと動くように少しトリックを加える必要がある。

具体的にいうと、特徴行列に1行加える必要がある。

 

なぜこれをする必要があるかは、重回帰分析の式を見るとわかりやすい。

 

y = b0 + b1*x1 + ..... + bn*xn

 

b0にx0がかかっていないが、実はx0が1だからである。よって本当の形は

 

y = b0*x0 + b1*x1 + .... + bn*xnとなる。

 

そして、今から使おうとしているstatsmodelライブラリはこのb0を考慮していない。

(ほとんどのライブラリ(例えばLinearRegression)ではb0は考慮されているのでx0を加える必要がない。)

だから、特徴行列に全ての要素が1のx0の列を加えてあげる必要がある。

 

こうすることで、statsmodelライブラリは重回帰分析の式を

 

y = b0 + b1*x1 + ..... + bn*xn

 

と認識することができる。

追加の仕方

np,appendを使う。

第一引数arrはarrの末尾に値を追加した配列を生成する。この時arr自体は変化しない。

第二引数valuesは追加したいもの。arrと同じ型にする。

第三引数axisは列に追加したいか行に追加したいか。0なら行に追加。1なら列に追加。

 

これでBackward eliminationをする準備が完了。

 

Backwardelimination

最初にすることは最終的に最適な特徴行列になる新しい行列を作ることだ。この行列には最終的に、従属変数の予測によくきく独立変数のみが入ることになる。

 

そして、backwardeliminationは最初に全ての独立変数を含み、そこから1つずつ重要でないものを消していく作業をする。

よって

X_opt = X[:,[0,1,2,3,4]]というふうに書く。

 

 

最初のステップ1は、モデルにとどまる有意水準(SL)を選択することです。

 

つまり、独立変数のp値が有意水準を下回る場合、独立変数がモデル内に留まるように、有意水準を選択する必要があります。ただし、独立変数のp値が有意水準のバーである場合、それはモデル内に留まらず、削除します。

 

単純にステップ1がすることはSLを選択すること。そしてその選択は0.05、つまり5%とする。

 

ステップ2ではFit the full model with  all possible predictorsを行う。

このall possible preditorsにあたるのが上で作ったX_optである。そしてFitはまだできていなかったので、新しいライブラリを使って行う。

import statsmodels.api as sm

このライブラリを使って重回帰モデルを機能XとYの将来の最適な行列に修正する。実際にフィットさせるにはyが必要になるからです。

 

regressor_OLS = sm.OLS(endog=y,exog=X_opt).fit()

このOLSはoridinary least squares (普通の最小二乗法)の略である。

第一引数endogには従属変数を入れる。

第二引数exogには観測の数とリグレッサーの数を持った配列を入れる。つまりX_opt。

 

 

ステップ3では最も高いp値を持つ独立変数を探す。そしてもしそのp値がSLの5%を超えていたらステップ4に進み、もし下回っていたらモデルの準備は完成ということになる。

 

そして高いp値を探すにはsummary関数を使う必要がある。これは、モデルを作るときに、重要となる全ての統計行列を含む優れたテーブルを返してくれます。そして、そのテーブルにはそれぞれの独立変数に対応するp値が含まれており、それをSLと比較することで、どれをモデルから削除するか決めれる。

regressor_OLS.summary()

出力

f:id:minmin_std:20191230184604p:plain



それぞれの独立変数に対してのp値が見れる。constはb0を表す。

重要なのはp値が低いほど、独立変数が従属変数に対してより重要であるということだ。

 

次にステップ4。独立変数の削除だ。

表を見て最も高いp値はx2の0.990である。これはSLの0.05を超えているため、この独立変数をモデルから削除する。インデックス番号で言うと、X_optの2番目の列にあたる。

X_opt = X[:,[0,1,3,4,5]]

これで二番目の列を削除した新たなX_optを作成

 

ステップ5

そして後は同じ手順でモデルを適用させテーブルを出す。そしてステップ3に戻る。
regressor_OLS = sm.OLS(endog=y,exog=X_opt).fit()
regressor_OLS.summary()

f:id:minmin_std:20191230184756p:plain

 

そして今度は新たなX_optのインデックス番号1番目の列が最もp値が高いのでこれを削除する。

 

f:id:minmin_std:20191230185234p:plain

60%なのでまた5%を超えているのでこの作業を続ける。

 

f:id:minmin_std:20191231000444p:plain

x1はR&D Spendだがこれは0.000となっていてprofitを予測するのに非常に強い独立変数であることがわかる。またp値は0になることはないので、これは0に非常に近いと言うことである。

x2は0.05よりわずかに大きい。SLを10%などと違う値を設定することでこの独立変数を維持することはできるが今回はこのルールに従う。

 

f:id:minmin_std:20191231001138p:plain

 

よって今回の最適な独立変数の組み合わせはR&D Spend1つのみと言う結果になった。

 

時間と共に指数関数的に増える従属関数を予測する時は、重回帰分析は向かない。