Python(StatsModels) で重回帰分析を理解し、分析の精度を上げる方法

こんにちは、データサイエンティストのたぬ（@tanuhack）です！

重回帰分析は『数字の予測』や『優先順位付け』に強く、ビジネスシーンにおけるデータ分析の中で、最も多用されている分析手法です。

記号	意味
	目的変数(量的データ)
〜	複数の説明変数(量的データでも質的データでも可)
〜	偏回帰係数
	定数項(切片)

数式だけだとわかりにくいと思うので、具体例を踏まえ考えてみましょう。

もし、ある町の家賃の相場(目的変数)を予測したいとして、『最寄り駅までの距離』『専有面積』『築年数』を説明変数として、重回帰分析を行った結果、次のような回帰式が得られました。

架空の不動産の回帰式

1か月の家賃(円)
= 27000 + (-2500×最寄り駅までの時間(分)) + (3000×専有面積()) + (-1500×築年数(年))

気になった物件が、最寄り駅まで徒歩10分、専有面積が30、築年数10年だったとすると…。

上のように重回帰分析によって得られた回帰式によって、77,000円くらいと大まかな家賃を予測することができます。

Pythonで重回帰分析をする方法として、scikit-learnを用いる方法とStatsModelsを用いる方法の2つが存在しますが、前者の方法では解析の結果から得られた重回帰式の精度を表す各指標が見れないので使いません。

そこで今回は、PythonのStatsModelsモジュールを使って、重回帰分析をする方法を紹介し、得られた予測の精度を上げるために行う方法も紹介します！！

この記事の著者

データサイエンティスト / ブロガー

たぬ

tanu

詳しいプロフィール

プロフィール

福岡在住のデータサイエンティスト。Pythonでルーチンワークを瞬殺する仕組みを作ることやデータを視覚化することが得意です。

StatsModelsモジュールの使い方

まず、ターミナルからStatsModelsモジュールをインストールします。

以下のプログラムは、StatsModelsを使う上でのテンプレです。適宜、変数名を入れ替えて見て下さい。

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

# 何かしらのデータ
df = XXXXXXXXXXXX

# get_dummies()で質的データも対応可能
x = pd.get_dummies(df[['説明変数1', '説明変数2', '説明変数3']]) # ←もちろん増減可能 
y = df['目的変数']

# 定数項(y切片)を必要とする線形回帰のモデル式ならば必須
X = sm.add_constant(x)

# 最小二乗法によるモデリング
model = smf.OLS(y, X)
result = model.fit()

# 重回帰分析の結果を表示
result.summary()

それでは、実際にデータセットを用いてStatsModelsモジュールの挙動を確かめてみましょう

データセットの確認

データセットは、SIGNATE(シグネイト)と呼ばれる日本版kaggleの勉強用に提供されているお弁当の需要予測データを使用します。

また、このCSVファイル名をtrain.csvとして保存し、pandas.DataFrameのメソッドread_csv()でdf変数に取り込んだものとします。

import pandas as pd
import numpy as np

df = pd.read_csv('train.csv')

【Python】PandasでCSVファイルを読み込み/書き出しする実践テクニック集

時系列に沿った、お弁当の販売個数の推移を折れ線グラフで確認しておきましょう。

from matplotlib import pyplot as plt
import seaborn as sns; sns.set(font='IPAexGothic')
%config InlineBackend.figure_formats = {'png', 'retina'}
%matplotlib inline

df.plot(x='datetime', y='y', figsize=(12,5), title='時系列に沿ったお弁当の販売個数推移')

StatsModelsで重回帰分析

STEP1：複数の説明変数を選ぶ

まず手始めに、説明変数『week』と『temperature』が、目的変数の『お弁当の販売個数』に影響があるのか、重回帰分析してみます。

df[['temperature','week']]

STEP2：説明変数の中に質的データが含まれていないか確認する

重回帰分析を行う前に確認しないといけないのが、説明変数の中に質的データが含まれているかどうかです。

今回のパターンでは、説明変数『week』が質的データでした。

選んだ説明変数の中に質的データが含まれる場合、そのままでは重回帰分析できないので、ダミー変数化という処理を行います。

STEP3：質的データをダミー変数化させる

質的データをダミー化させる方法として、pandas.DataFrameのメソッドget_dummies()を用います。

pd.get_dummies(df[['temperature','week']])

get_dummies()の偉いところは、質的だった場合ダミー変数化して量的だった場合そのまま、自動でやってくれることですね。

では、整ったので重回帰分析を行ってみましょう＼(^o^)／

STEP4：重回帰分析を行う

import pandas as pd
import statsmodels.api as sm 
import statsmodels.formula.api as smf

df = pd.read_csv('train.csv')

x = pd.get_dummies(df[['temperature','week']]) # 説明変数
y = df['y'] # 目的変数

# 定数項(y切片)を必要とする線形回帰のモデル式ならば必須
X = sm.add_constant(x)

# 最小二乗法でモデル化
model = smf.OLS(y, X)
result = model.fit()

# 重回帰分析の結果を表示する
result.summary()

結果の見方

summary()によって、解析結果が表示されたら、まず最初に見るべき3つの指標があります。

1. Adj. R-squared：自由度調整済み決定係数
2. t：各説明変数ごとの統計量
3. p>|t|：各説明変数ごとのp値

STEP1：『Adj. R-squared：自由度調整済み決定係数』を確認する

この値は、一般的に『自由度調整済み決定係数』と呼ばれ、回帰式全体の精度、説明力を示しています。

説明変数が増えると、決定係数(R-squared)は大きくなると言われているので、重回帰分析では『自由度調整済み決定係数(Adj. R-squared)』の方を使いましょう。

今回の自由度調整済み決定係数は0.447で、なんとも言えない精度ですね。目標は0.85くらいだと言われています。

この数値の結果から、説明変数の選定が十分ではないことが分かります。

STEP2：『t：各説明変数ごとの統計量』を確認する

この数値の絶対値が大きければ大きいほど、目的変数に与える影響が大きいという指標になります。

constというカラム名は、X = sm.add_constant(x)で生成した『y切片』のことなので、無視して大丈夫です。

時点でt値の絶対値が高いのは、temperatureであるのが分かります。つまり『気温が上がればお弁当の販売個数は減る』ことに関係がありそうですね。

STEP3：『p>|t|：各説明変数ごとのp値』を確認する

このカラムには、1つ1つの説明変数毎のp値が表示されています。

p値の説明は不要だと思いますが、念のために…。p値は、0に近づけば近づくほど統計的な意味があり、逆に0.05よりも大きい値であればたまたまである可能性が高いとして、その説明変数は採用しないほうが良いとされています。

今回はすべて、0.05より低い値なので、全部説明変数として採択しても大丈夫そうですね。

結果から回帰式を作る

今回選んだ説明変数の重回帰分析で得られた、自由度調整済み決定係数の値は44%と精度は低いですが、とりあえず結果から回帰式を作る方法を紹介します。

『coef』に注目します。

気温20度で月曜日だった場合、93個くらい売れそうだなと予測を立てることができます！（精度44%ですが、笑）

ちなみに精度44.7%のグラフは下図のようになります。

pred = result.predict(X)
df['pred'] = pred

df.plot(y=['y','pred'], x='datetime', figsize=(12,5), title='精度44.7％のグラフ')

全然あってないですね！笑

次からは、重回帰分析の精度を上げるためにした方がいい方法を紹介します。

重回帰分析の精度を上げる

精度が低いモデルを作成しても何の価値もありません。

質の高いアウトカム(因果関係の推察、未来の予測)を実現するためには、モデリング(重回帰分析)にかける前のデータの質が大切になってきます。

『データの質を上げる』とは、モデリングの精度を上げることであり、その精度を上げるためにはデータの整形すなわち『前処理』が大切です。

データ分析の８割は、前処理だと言っても過言ではないでしょう。つまり、この泥臭い前処理を怠らず分析の精度の質を高めていくことこそが、データ分析する人の使命であり、力量が試されます。

Python(Pandas)でデータ分析するときに使う基本操作(前処理)まとめ

ステップワイズ法

説明変数の候補がたくさんあれば、『とりあえずすべての説明変数を重回帰分析にかけてp値が小さく、t値の絶対値が大きいものを探索する』というやり方をステップワイズ法と言います。

このときのp値は0.05以下ではなく、0.1以下で見ることが慣例のようです。

STEP1：カラムが『datetime』以外のすべての説明変数を入れて重回帰分析を行う

# すべての欠損値を0で置換する
df.fillna(0, inplace=True)
df.loc[df['precipitation'].str.contains('--'), 'precipitation'] = '0'

x = pd.get_dummies(df[['week', 'soldout', 'name', 'kcal', 'remarks', 'event', 'payday', 'weather', 'precipitation', 'temperature']]) # 説明変数
y = df['y'] # 目的変数

# 定数項(y切片)を必要とする線形回帰のモデル式ならば必須
X = sm.add_constant(x)

# 最小二乗法でモデル化
model = smf.OLS(y, X)
result = model.fit()

# 重回帰分析の結果を表示する
result.summary()

説明変数を闇雲に増やすと重回帰式の精度は、44%→63%まで向上しました。

STEP2：p値が全て0.05(場合によって0.1)以下のもの以外を説明変数から取り除く

x = pd.get_dummies(df[['week', 'kcal', 'remarks', 'event', 'temperature']]) # 説明変数
x = x.drop(['remarks_0','remarks_手作りの味','remarks_料理長のこだわりメニュー','remarks_酢豚（28食）、カレー（85食）','remarks_鶏のレモンペッパー焼（50食）、カレー（42食）','week_木','week_金'], axis=1)
y = df['y'] # 目的変数

# 定数項(y切片)を必要とする線形回帰のモデル式ならば必須
X = sm.add_constant(x)

# 最小二乗法でモデル化
model = smf.OLS(y, X)
result = model.fit()

# 重回帰分析の結果を表示する
result.summary()

不要な説明変数を削除すると重回帰式の精度は、63%→71%まで向上しました。

STEP3：ステップワイズ法後の精度をグラフで確認する

pred = result.predict(X)
df['pred'] = pred

df.plot(y=['y','pred'], x='datetime', figsize=(12,5), title='精度71.6％のグラフ')

ステップワイズ法で重回帰分析の精度を上げていく様子が何となく掴めたのではないでしょうか。

ステップワイズ法+仮説【おすすめ】

説明変数が多いと、ステップワイズ法だけ使っても精度の高いモデルを作ることは中々できません。そこでおすすめしたいのが、ステップワイズ法に『仮説』を取り入れた方法です。

具体的にどうやって仮説を思いつくかと言うと、『データの可視化』を行い、傾向や偏り、変曲点、スパイクを目で捉えます。

STEP1：データの可視化により仮説を立てる

下の折れ線グラフからどういう仮説を立てられるでしょうか。少しだけ考えてみて下さい。

などなど…。

もちろん、これら以外にもいろんなことが考えられます。正解なんてないので、考えられる限りを書き出してみてくださいね。

他にもいろんなデータの可視化の手段があります。目的を考えてグラフを使い分けて、いろんな仮説を立ててみましょう！

STEP2：仮説から販売個数に影響がありそうなものを見つける

仮説を元に説明変数を探していると、どうやらメニュー名に『カレー』が含まれるか否かで販売個数に影響がありそうです。

# メニュー名に『カレー』を含む
df[df['name'].str.contains('カレー')]['y'].mean()

# メニュー名に『カレー』を含まない
df[df['name'].str.match('^(?!.*カレー).*$')]['y'].mean()

STEP3：説明変数に追加する

今回でいうと、name列にカレーという文字があれば1を無ければ0を返すcurry列を作ります。

pandas.Seriesのメソッドapply()を使って、1行ずつ追加させます。

def curry(x):
    if 'カレー' in x:
        return 1
    else :
        return 0
    
df['curry'] = df['name'].apply(lambda x : curry(x))

STEP4：同様に仮説を立て、新たな列として追加する

# remarks(特記事項)が『お楽しみメニュー』のときも売れてそう！
def fun(x):
    if x=='お楽しみメニュー':
        return 1
    else :
        return 0

df['fun'] = df['remarks'].apply(lambda x : fun(x))

# datetimeのカラムから年のカラムを作成
df['year'] = df['datetime'].apply(lambda x : x.split('-')[0])
df['year'] = df['year'].astype(np.int)

# datetimeのカラムから月のカラムを作成
df['month'] = df['datetime'].apply(lambda x : x.split('-')[1])
df['month'] = df['month'].astype(np.int)

STEP5：重回帰分析を行う

x = pd.get_dummies(df[['year', 'month', 'week', 'kcal', 'fun', 'curry', 'weather', 'temperature']]) # 説明変数
y = df['y'] # 目的変数

# 定数項(y切片)を必要とする線形回帰のモデル式ならば必須
X = sm.add_constant(x)

# 最小二乗法でモデル化
model = smf.OLS(y, X)
result = model.fit()

# 重回帰分析の結果を表示する
result.summary()

ステップワイズ法+仮説で、重回帰式の精度は71%→78%まで向上しました。

精度としては、85%までは上げたいところですが、残りはあなた自身でトライしてみてください！

STEP6：精度を目(グラフ)で確かめる

pred = result.predict(X)
df['pred'] = pred

df.plot(y=['y','pred'], x='datetime', figsize=(12,5), title='精度77.9％のグラフ')

一番重要な説明変数の見つけ方

『ステップワイズ法+仮説』を使って、モデルの精度が高く目的変数に影響がある複数の説明変数を見つけることができました。

今回得られた説明変数の中で、目的変数に最も影響を与えるのは、その中でも統計量の絶対値が大きいものだと紹介しましたね。

しかし「統計量の絶対値が大きいから」といって、それが一番重要な説明変数かと言われると、そうではありません。具体例を見てみましょう。

例えば、上の表のような結果が得られたとします。

そして、『どれぐらい説明変数を動かせる余地があり、また実際に説明変数をどれくらい動かせる手段があるのか』を考えなければなりません。

「temperatureとかweekとかweatherが統計量が高いので、これらが重要な説明変数です。」と報告されても、人間である以上天候や曜日を自在に操ることはできないので、果たしてこの分析結果は正しいと言えるでしょうか？

どう考えても、言えないですよね。

この結果の場合の重要な説明変数は、name(お弁当の名前)です。メニューを分析して、明らかに好まれているメニューとそうではないメニューを見つけ出し、変えていくことがお弁当の販売個数に対して、『最大効果』は低いかもしれませんが、確実にインパクトを与えられますね。

重回帰分析で注意すること

ステップワイズ法と仮説を組み合わせて、重回帰分析の精度を上げる方法を紹介しました。

ですが、精度だけ上げればそれでいいか、と言われればそうでもありません！！！

重回帰分析には、これから紹介する落とし穴が2つ存在します。それぞれ見てみましょう。

1. 過学習(オーバーフィッティング)
2. 多重共線性(マルチコ)

過学習(オーバーフィッティング)

基本的に説明変数が増えれば増えるほど、重回帰式の精度は高くなると紹介しましたが、それだけがいいことばかりとは限りません。

当てはめの精度は高いのに、予測精度が低くなることを過学習(オーバーフィッティング)と言います。

過学習になる原因は、『手持ちデータ』に過剰に適合しすぎたモデルを構築してしまったことです。こうなると、いまある『検証用データには当てはまりが良い』が『予測したい新しいデータに回帰式を当てはめると、当てはまりが悪くなる』といった減少が起きてしまいます。

過学習を回避するためには一般的に次に紹介する『クロスバリデーション法』を用います。

クロスバリデーション法とは

『回帰式を求める分析用のデータ』と、『その当てはまりの良さを確認するためのデータ』の2パターンを用意します。

今回、重回帰分析用に使用したデータセットには、回帰式を求める『train.csv』と当てはまりの良さを確認する『test.csv』の2つが用意されているので、test.csvを使います。

テスト用のデータセットがない場合

1つのデータを完全ランダムに半分抽出し、片方のデータで回帰式を求め、もう片方のデータで当てはまりの良さを検証するやり方を用います。(回帰式60%、テスト40%とかでも大丈夫です。)

# 回帰式を求める分析用のデータ
train = df.sample(frac=0.5, random_state=0)

# 回帰式の当てはまりの良さを確認するためのテストデータ
test = df.drop(train.index)

ただ、完全ランダムだと『時間の経過』などを説明変数に入れたい場合、不向きです。

時系列データの場合は、時間順にソートされているのを確認して半分にデータを分けた方が良いかもしれませんね。

# 時系列データの場合

# 回帰式を求める分析用のデータ
train = df[0:len(df)//2]

# 回帰式の当てはまりの良さを確認するためのテストデータ
test = df.drop(train.index)

重回帰式で得られた予測値をテストデータに代入させる

StatsModelsのメソッドformula.api.OLS().fit().predict()を使えば、予測値を求めることができます。

ステップワイズ法+仮説で作ったプログラムで試してみます。

import pandas as pd
import statsmodels.api as sm 
import statsmodels.formula.api as smf
from matplotlib import pyplot as plt
import seaborn as sns; sns.set(font='IPAexGothic')
%config InlineBackend.figure_formats = {'png', 'retina'}
%matplotlib inline

df = pd.read_csv('train.csv')

# すべての欠損値を0で置換する
df.fillna(0, inplace=True)
df.loc[df['precipitation'].str.contains('--'), 'precipitation'] = '0'

def curry(x):
    if 'カレー' in x:
        return 1
    else :
        return 0

def fun(x):
    if x=='お楽しみメニュー':
        return 1
    else :
        return 0

df['fun'] = df['remarks'].apply(lambda x : fun(x))
df['curry'] = df['name'].apply(lambda x : curry(x))

# datetimeのカラムから年のカラムを作成
df['year'] = df['datetime'].apply(lambda x : x.split('-')[0])
df['year'] = df['year'].astype(np.int)

# datetimeのカラムから月のカラムを作成
df['month'] = df['datetime'].apply(lambda x : x.split('-')[1])
df['month'] = df['month'].astype(np.int)

x = pd.get_dummies(df[['year', 'month', 'week', 'kcal', 'fun', 'curry', 'weather', 'temperature']]) # 説明変数
y = df['y'] # 目的変数

# 定数項(y切片)を必要とする線形回帰のモデル式ならば必須
X = sm.add_constant(x)

# 最小二乗法でモデル化
model = smf.OLS(y, X)
result = model.fit()

# 重回帰分析の結果を表示する
result.summary()

テストデータと回帰式を作るデータは、全く同じ説明変数を揃えなければいけません。今回で言うと、curry、year、month、funを新たに説明変数として追加したので、テストデータにもそれらを追加します

# すべての欠損値を0で置換する
test.fillna(0, inplace=True)
test.loc[test['precipitation'].str.contains('--'), 'precipitation'] = '0'

# 回帰式を作るデータと同じ説明変数を追加する
test['curry'] = test['name'].apply(lambda x : curry(x))
test['fun'] = test['remarks'].apply(lambda x : fun(x))
test['year'] = test['datetime'].apply(lambda x : x.split('-')[0])
test['year'] = test['year'].astype(np.int)
test['month'] = test['datetime'].apply(lambda x : x.split('-')[1])
test['month'] = test['month'].astype(np.int)

testX =  pd.get_dummies(test[['year', 'month', 'week', 'kcal', 'fun', 'curry', 'weather', 'temperature']]) # 説明変数
testX['const'] = 1.0 # 切片を追加（sm.add_constant()が使えないので）

# 回帰式を作るデータのカラム順と同じになるように並び替えする
# （ここでエラーが発生したら説明変数が足りていません。適宜追加してください。）
# testX['weather_雪'] = 0
# testX['weather_雷電'] = 0
testX = testX.ix[:,X.columns]

# 重回帰式で得られた予測値をテストデータに代入させる
pred = result.predict(testX)
test['pred'] = pred

この例で言うと、予測値がずれまくってることは無いと判断できます。（ひどいときはマイナスがついたり、桁がおかしかったりします）

こんな感じで、回帰式を求めるデータとテストデータ間で大きな乖離が無いか確認すれば、過学習に気づいて、より精度の高い重回帰式を作ることができるでしょう！

多重共線性(マルチコ)

多重共線性とは、説明変数間で非常に強い相関があることを指し、この値が大きいと回帰係数の分散が大きくなり、モデルの予測結果が悪くなることが知られています。

ただし、重回帰分析を行う目的が『因果関係の洞察』ではなく、『予測』であれば、気にしなくて大丈夫です。

summary()の結果でいう、Cond. No.が多重共線性をチェックする指標になります。

この例でいうと、と非常に大きい数字になっているので、多重共線性の可能性が大いにありえます。

多重共線性を避けるために、その説明変数を除外するかしないか考えたときに参考になる指標として、分散拡大係数(VIF：Variance Inflation Factor)というものがあります。

Pythonでは、StatsModelsのメソッドstatsmodels.stats.outliers_influence.variance_inflation_factorを使って、VIFを計算できます。

from statsmodels.stats.outliers_influence import variance_inflation_factor as vif

num_cols = model.exog.shape[1] # 説明変数の列数
vifs = [vif(model.exog, i) for i in range(0, num_cols)]

pd.DataFrame(vifs, index=model.exog_names, columns=['VIF'])

一般的にVIFの値が10(公式のリファレンスでは、5)を超えると、依存関係が強いため、適切な重回帰分析ができないと言われています。

今回でいうと、ダミー変数化で作成した『week』の列のVIF値がすべて『inf』となっており、依存関係が非常に強いです。

繰り返しになりますが、重回帰分析の目的が『因果関係の洞察』であれば説明変数から除外したほうが無難であり、『予測』が目的であれば除外しなくても大丈夫です。

さいごに

今回はPythonのStatsModelsモジュールで重回帰分析をして、その精度を上げていく方法を紹介しました。

説明変数の選択は、ステップワイズ法+仮説で行い、過学習と多重共線性に気を付けていれば、概ねうまくいきます。

それでは〜