Python3で線形モデルによる回帰分析とプロット

 2018/02/25    Pythonデータ解析入門    

この記事ではPython3で線形モデルによる回帰分析のやり方を分かりやすくご紹介します。サンプルcsvファイルを説明用に使いますので、記事を読みながら一緒に手を動かしたい方はぜひダウンロードして使って下さい。

データは以下のような形です。

    number  age  blood_pressure  lung_capacity sex  weight  disease
0        1   22             110           4300   M      79        1
1        2   23             128           4500   M      65        1
2        3   24             104           3900   F      53        0
3        4   25             112           3000   F      45        0
4        5   27             108           4800   M      80        0
5        6   28             126           3800   F      50        0

まずは単回帰分析

使うライブラリは、statsmodelsです。これを用いて最小二乗法を用いた線形モデルによる回帰分析を行います。

今回は、X血圧(blood_presssure)からY肺活量(lung_capacity)を予測することを考えましょう。それではいきなりですが、プログラムは以下のようになります。

import statsmodels.api as sm
import pandas as pd

# データを読み込む
data = pd.read_csv("sample.csv")

#回帰分析に使うデータの指定
x = data.iloc[:,[2]] #説明変数
y = data.iloc[:,3] #目的変数

#全要素が1の列を説明変数の先頭に追加(絶対必要!!)
X = sm.add_constant(x)

#モデルの設定
model = sm.OLS(y, X)

#回帰分析の実行
results = model.fit()

#結果の詳細を表示
print(results.summary())

このプログラムの説明を簡単にすると、まずdataにpadasのデータフレーム型でデータを読み込み、説明変数(血圧の列)と目的変数(肺活量の列)をそれぞれ変数x,yに代入しています。

また、 X = sm.add_constant(x) としている部分は、説明変数の一列目に新しく全要素が1.0の列を追加しています。これは、もし切片を必要とする線形回帰のモデル式ならば必ず必要な部分で、これを入れないと正しく回帰式が作成されません

model = sm.OLS(y, X) ではモデルを設定しています。今回は最小二乗法なのでOLSとしましたが、WLSなど他のモデルでも出来ます。

results = model.fit()で回帰分析を実行します。最後に結果 results.summary() をプリントしました。

結果は次のようになりました。

                            OLS Regression Results                            
==============================================================================
Dep. Variable:          lung_capacity   R-squared:                       0.235
Model:                            OLS   Adj. R-squared:                  0.208
Method:                 Least Squares   F-statistic:                     8.620
Date:                Sat, 03 Feb 2018   Prob (F-statistic):            0.00658
Time:                        17:40:20   Log-Likelihood:                -231.08
No. Observations:                  30   AIC:                             466.2
Df Residuals:                      28   BIC:                             469.0
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
const           6491.4949   1010.471      6.424      0.000    4421.640    8561.350
blood_pressure   -23.6714      8.062     -2.936      0.007     -40.187      -7.156
==============================================================================
Omnibus:                        2.242   Durbin-Watson:                   2.339
Prob(Omnibus):                  0.326   Jarque-Bera (JB):                1.287
Skew:                           0.166   Prob(JB):                        0.526
Kurtosis:                       2.042   Cond. No.                     1.25e+03
==============================================================================

表になっているので非常に見やすいのがメリットですね。ここでは色々な指標が出ましたが、中でも特に重要なものを以下の表にまとめました。

名称意味
R-squared決定係数。これが1に近いほど精度の高い分析と言える。
Adj. R-squared自由度調整済み決定係数。説明変数が多い時は決定係数の代わりに用いる。
AICモデルの当てはまり度を示す。小さいほど精度が高い。相対的な値である。
coef回帰係数
std err二乗誤差
tt値。
pp値。有意水準以下の値を取れば、回帰係数の有意性が言える。
[0.025 0.975]95%信頼区間。

p値が0.05よりもはるかに小さい値なので、有意水準5%で回帰係数に統計的な優位性が言えています。ただし、決定係数は非常に小さいのでモデルの当てはまり具合はそんなに良くはなそうです。重回帰分析についての説明の後、この単回帰分析の結果をプロットして、可視化してみます。

次に説明変数を増やして重回帰分析

重回帰分析を行う場合は、説明変数を増やすだけです。いかに重回帰分析のサンプルプログラムを示します。

import statsmodels.api as sm
import pandas as pd

# データを読み込む
data = pd.read_csv("sample.csv")

#回帰分析に使うデータの指定
x = data.iloc[:,[1,2,5,6]] #説明変数
y = data.iloc[:,3] #目的変数

#全要素が1の列を説明変数の先頭に追加,切片をつけるために必ず必要
X = sm.add_constant(x)

#モデルの設定
model = sm.OLS(y, X)

#回帰分析の実行
results = model.fit()

#結果の詳細を表示
print(results.summary())

先ほどのプログラムのxの変数に入れる列を増やすだけで重回帰分析に変わります。結果は以下のようになりました。

                           OLS Regression Results                            
==============================================================================
Dep. Variable:          lung_capacity   R-squared:                       0.838
Model:                            OLS   Adj. R-squared:                  0.812
Method:                 Least Squares   F-statistic:                     32.27
Date:                Sat, 24 Feb 2018   Prob (F-statistic):           1.54e-09
Time:                        23:44:46   Log-Likelihood:                -207.82
No. Observations:                  30   AIC:                             425.6
Df Residuals:                      25   BIC:                             432.7
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
const           2585.2366    648.246      3.988      0.001    1250.150    3920.324
age              -26.9601      5.461     -4.937      0.000     -38.207     -15.713
blood_pressure     1.8499      5.286      0.350      0.729      -9.037      12.736
weight            29.1946      4.096      7.127      0.000      20.758      37.631
disease          102.7966    110.368      0.931      0.361    -124.510     330.103
==============================================================================
Omnibus:                        5.566   Durbin-Watson:                   2.067
Prob(Omnibus):                  0.062   Jarque-Bera (JB):                1.903
Skew:                          -0.135   Prob(JB):                        0.386
Kurtosis:                       1.796   Cond. No.                     1.88e+03
==============================================================================

ただ、説明変数が増えているだけで、基本的な結果の見方は先ほどの単回帰分析の場合と同じです。決定係数は説明変数が増えれば増えるほど1に近づく性質があるので、説明変数が多い場合は、決定係数ではなく自由度調整済み決定係数の値をより重要視するようにしましょう。

結果をプロット

先ほどの単回帰分析の結果をプロットしたい場合、以下のようにプログラムを記述します。今回はライブラリとしてmatplotlib.pyplotを使いました。

import pandas as pd #データフレームを扱う用
import statsmodels.api as sm #回帰分析用
import matplotlib.pyplot as plt #プロット用

# データを読み込む
data = pd.read_csv("sample.csv")

#回帰分析に使うデータの指定
x = data.iloc[:,[2]] #説明変数
y = data.iloc[:,3] #目的変数

#全要素が1の列を説明変数の先頭に追加,切片をつけるために必ず必要
X = sm.add_constant(x)

#モデルの設定
model = sm.OLS(y, X)

#回帰分析の実行
results = model.fit()

#回帰係数と切片の推定値
a = results.params[0]
b = results.params[1]

#標本値をプロット
plt.plot(x, y,"o")

#回帰直線をプロット
plt.plot(x, a+b*x)

#プロットの表示
plt.show()

流れとしては、回帰分析→回帰係数と切片の値をa,bという変数に格納→標本値をプロット→回帰直線をプロット→プロットの表示という感じです。下図のようになりました。

  • スポンサーリンク

  • 関連コンテンツ

  • コメントを残す

    メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

    CAPTCHA