2018/02/25
2020/04/14
Python3で線形モデルによる回帰分析とプロット
この記事では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 | 二乗誤差 |
| t | t値。 |
| p | p値。有意水準以下の値を取れば、回帰係数の有意性が言える。 |
| [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という変数に格納→標本値をプロット→回帰直線をプロット→プロットの表示という感じです。下図のようになりました。
Recommended