線形回帰をstatsmodelsで実行。scikit-learnの結果と比べてみる。

statsmodelsとscikit-learn

Pythonで機械学習といえばscikit-learn。ですが、まずは統計学寄りのstatsmodelから触ってみる。statsmodelは予測モデルの表示に加えて、その名の通り、統計的な情報、例えば検定結果も計算して表示する。t値とかp値とか。

scikit-learnの実行例があったので、それと同じことをstatsmodelで書いて、結果（＝予測モデル）が一致するかどうかを確認する。

statsmodelsとscikit-learn
- scikit-learnの実行例
線形回帰とは
線形回帰をやってみる
今日の単語帳

scikit-learnの実行例

scikit-learn で線形回帰 (単回帰分析・重回帰分析) – Python でデータサイエンス

線形回帰とは

リンク先にそのままあるので詳細はそちらで。

y=b1x1+b2x2+b3x3+...+bnxn+b

x1, x2, x3, ... xnは、説明変数
yは、目的変数（説明変数を使って目的変数を予測する）
b1, b2, b3, .... bn は、係数
bはy切片

説明変数が１つだけの場合は「単回帰分析」、２つ以上ある場合は「重回帰分析」と呼ぶ。

線形回帰をやってみる

使用するデータ

ワインの品質（赤ワイン）のデータセットを使用する。

UCI Machine Learning Repository: Wine Quality Data Set

データはこんな感じ（先頭5行を表示）。

f:id:mkmkipy:20181013035937j:plain

quality列が品質を示す。定義より、qualityは0～10。
しかし、実際の値とそれぞれの件数は、print(df['quality'].value_counts()) より

5 681
6 638
7 199
4 53
8 18
3 10

つまり、quality列の値は、3, 4, 5, 6, 7, 8 の6種類である。

単回帰分析

ここでは、dencity(濃度)から、alcohol(アルコール度数)を予測する。

import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# データセットの読み込み
df = pd.read_csv("winequality-red.csv", sep=";")

### 単回帰
### ここでは、dencity(濃度)から、alcohol(アルコール度数)を予測する

# 回帰分析に使うデータの指定
x = df[['density']] #説明変数
y = df[['alcohol']] #目的変数
 
# 全要素が1.0の列を説明変数の先頭に追加
X = sm.add_constant(x)
 
# モデルの設定(OLS＝最小二乗法を指定)
model = sm.OLS(y, X)
 
# 回帰分析の実行
results = model.fit()
 
# 結果の詳細を表示
results.summary()

ポイント1：回帰分析に使うデータの指定

説明変数、目的変数はそれぞれ次のように指定。

x = df[['density']] #説明変数
y = df[['alcohol']] #目的変数

ポイント２：切片をモデル式に含める

statsmodelでは、切片を必要とする線形回帰のモデル式の場合（いや普通そうでしょ！）は、全要素が1.0の列を説明変数の先頭に追加する必要がある。

X = sm.add_constant(x)

結果の確認

参考文献にある単回帰分析の例にある回帰式の係数と切片が、今回の例と同じことを確認した。

[alcohol] = -280.16382307 × [density] + 289.675343383

重回帰分析

説明変数が2つ以上あるものを重回帰分析と呼ぶ。ここでは「quality」を目的変数に、「quality」以外を説明変数として、重回帰分析を行う。

import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# UCI Machine Leaning Repository 「Wine Quality Data Set (ワインの品質)」の赤ワインのデータセット
df = pd.read_csv("http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")

### 重回帰分析
### ここでは、quality（品質）以外の列から、qualityを予測する

# 回帰分析に使うデータの指定
x = df[['fixed acidity','volatile acidity','citric acid','residual sugar','chlorides','free sulfur dioxide','total sulfur dioxide','density','pH','sulphates','alcohol']] #説明変数
y = df[['quality']] #目的変数
 
# 全要素が1.0の列を説明変数の先頭に追加
# 切片を必要とする線形回帰のモデル式ならば必ず必要な部分で、これを入れないと正しく回帰式が作成されない
X = sm.add_constant(x)

# モデルの設定(OLS＝最小二乗法を指定)
model = sm.OLS(y, X)
 
# 回帰分析の実行
results = model.fit()
 
# 結果の詳細を表示
results.summary()

説明変数に複数列を指定する

重回帰分析の場合、データの指定で、x（説明変数）に複数列指定する（↓）。それ以外は単回帰の場合と同様でOK。

x = df[['fixed acidity','volatile acidity','citric acid','residual sugar','chlorides','free sulfur dioxide','total sulfur dioxide','density','pH','sulphates','alcohol']] #説明変数]]

結果の確認

参考文献にある重回帰分析の例にある回帰式の各係数と切片が、今回の例と同じことを確認した。

[quality] = -17.881164 × [density] + -1.874225 × [chlorides] +
-1.083590 × [volatile acidity] + -0.413653 × [pH] +
-0.182564 × [citric acid] + -0.003265 × [total sulfur dioxide] +
0.004361 × [free sulfur dioxide] + 0.016331 × [residual sugar] +
0.024991 × [fixed acidity] + 0.276198 × [alcohol] +
0.916334 × [sulphates] + 21.9652084495

もう少し簡単な書き方がわかったので書き直したのがこちら↓↓

mkmkipy.hatenablog.com