はじめに
今回は文部科学省のページで公開されている情報Ⅰの教員研修用教材の「第3章 情報とデータサイエンス 前半」内の「重回帰分析とモデルの決定」についてみていきたいと思います。
「重回帰分析とモデルの決定」の前章ではExcelとPythonを使用して、データの整形・前処理・操作・欠損値と異常値の取扱いなどについて扱っていますが、この章からは「重回帰分析とモデルの決定」を行っていきます。
しかしRでの実装例しか書いていないので、pythonでの実装について考えていきたいと思います。
教材
高等学校情報科「情報Ⅱ」教員研修用教材（本編）：文部科学省
第3章 情報とデータサイエンス 前半 (PDF:8.9MB) PDF
環境
- ipython
- Colaboratory - Google Colab
概要
教材の重回帰分析とモデルの決定内のp131から説明されている「（2）統計ソフト R での重回帰分析の実行」
重回帰分析のコンピュータでの実行と出力
の箇所でRで書かれてる実装例をpythonに書き換えていきたいと思います
前提
今回、前章までで分析を行う対象のデータはすでに整形を行った以下のファイルを使用しました。
high_male_data.csv
重回帰分析とモデルの決定
重回帰分析の実行
pythonのソースコード
stasmodelには、様々な線形回帰モデルがあります。 基本的なもの（最小二乗法、OLS）からより複雑なもの（反復再重み付け最小二乗法、IRLS）まであります． stasmodelの線形モデルには、2つの主なインターフェースがあります。
配列ベースのものと、formula式ベースのものです。
(第2版 Pythonによるデータ分析入門 NumPy,pandasを使ったデータ処理 より)
教材のコードは、最小二乗法(OLS)を使用しているので、ここでもOLSを使います。
まずは、教材のソースコードに近いformula式ベースのコードと配列ベースのコード両方見ていきたいと思います。
formula式
import pandas as pd
import statsmodels.formula.api as smf
high_male = pd.read_csv('/content/high_male_data.csv')
model = smf.ols('X50m走 ~ 立ち幅跳び + ハンドボール投げ + 握力得点 + 上体起こし得点', data = high_male)
results = model.fit()
print(results.summary())
配列ベース
import pandas as pd
import statsmodels.api as sm
high_male = pd.read_csv('/content/high_male_data.csv')
x = high_male[['立ち幅跳び', 'ハンドボール投げ', '握力得点','上体起こし得点']]
y = high_male['X50m走']
x = sm.add_constant(x)
model = sm.OLS(high_male['X50m走'], x)
results = model.fit()
print(results.summary())
pythonによる出力結果
formula式
OLS Regression Results
==============================================================================
Dep. Variable: X50m走 R-squared: 0.525
Model: OLS Adj. R-squared: 0.510
Method: Least Squares F-statistic: 36.18
Date: Mon, 20 Jul 2020 Prob (F-statistic): 2.39e-20
Time: 04:16:26 Log-Likelihood: -41.639
No. Observations: 136 AIC: 93.28
Df Residuals: 131 BIC: 107.8
Df Model: 4
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 10.8194 0.325 33.331 0.000 10.177 11.462
立ち幅跳び -0.0120 0.002 -7.648 0.000 -0.015 -0.009
ハンドボール投げ -0.0144 0.006 -2.367 0.019 -0.026 -0.002
握力得点 -0.0402 0.024 -1.677 0.096 -0.088 0.007
上体起こし得点 -0.0255 0.020 -1.264 0.208 -0.065 0.014
==============================================================================
Omnibus: 3.064 Durbin-Watson: 1.813
Prob(Omnibus): 0.216 Jarque-Bera (JB): 2.955
Skew: 0.359 Prob(JB): 0.228
Kurtosis: 2.927 Cond. No. 2.61e+03
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.61e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
配列ベース
OLS Regression Results
==============================================================================
Dep. Variable: X50m走 R-squared: 0.525
Model: OLS Adj. R-squared: 0.510
Method: Least Squares F-statistic: 36.18
Date: Mon, 20 Jul 2020 Prob (F-statistic): 2.39e-20
Time: 04:15:11 Log-Likelihood: -41.639
No. Observations: 136 AIC: 93.28
Df Residuals: 131 BIC: 107.8
Df Model: 4
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 10.8194 0.325 33.331 0.000 10.177 11.462
立ち幅跳び -0.0120 0.002 -7.648 0.000 -0.015 -0.009
ハンドボール投げ -0.0144 0.006 -2.367 0.019 -0.026 -0.002
握力得点 -0.0402 0.024 -1.677 0.096 -0.088 0.007
上体起こし得点 -0.0255 0.020 -1.264 0.208 -0.065 0.014
==============================================================================
Omnibus: 3.064 Durbin-Watson: 1.813
Prob(Omnibus): 0.216 Jarque-Bera (JB): 2.955
Skew: 0.359 Prob(JB): 0.228
Kurtosis: 2.927 Cond. No. 2.61e+03
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.61e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
- coef:係数の推定値
- R-squared: 寄与率(決定係数)
- Adj. R-squared: 自由度修正済みR2
そのほか、LinearRegressionを使う方法などもあります。
今回、細かいパラメータについての説明は省きます。
[参考]Rのソースコード(教材より)
# データの読み込み
high_male <- read.csv("high_male_data.csv")
# 4 つの説明変数による重回帰分析
res <- lm(X50m走 ~ 立ち幅跳び + ハンドボール投げ + 握力得点 + 上体起こし得点, data = high_male)
summary(res)
lm()は引数パラメータweights = nullの場合、OLS使用(ordinary least squares is used.)
[参考]Rによる出力結果
Call:
lm(formula = X50m走 ~ 立ち幅跳び + ハンドボール投げ +
握力得点 + 上体起こし得点, data = high_male)
Residuals:
Min 1Q Median 3Q Max
-0.66229 -0.22824 -0.03057 0.22797 0.99194
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.819417 0.324609 33.331 < 2e-16 ***
立ち幅跳び -0.012005 0.001570 -7.648 3.89e-12 ***
ハンドボール投げ -0.014393 0.006081 -2.367 0.0194 *
握力得点 -0.040185 0.023969 -1.677 0.0960 .
上体起こし得点 -0.025489 0.020163 -1.264 0.2084
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3349 on 131 degrees of freedom
Multiple R-squared: 0.5249, Adjusted R-squared: 0.5104
F-statistic: 36.18 on 4 and 131 DF, p-value: < 2.2e-16
- Estimate:係数の推定値
- Residual standard error:（残差）標準誤差
- degrees of freedom（残差）自由度
- Multiple R-squared:寄与率 R2
- Adjusted R-squared:自由度修正済みR2
モデル選択（モデリング）
pythonのソースコード
pythonによるステップワイズ法によって、モデル選択を行っていきます。
私の調べ方が悪く、pythonモジュールなどでステップワイズ法が行えるようなものが見つからなかったので、今回はAIC(赤池の情報規準)の値が最も小さくなるようなモデルを選択できるように実装します。
import pandas as pd
import statsmodels.formula.api as smf
high_male = pd.read_csv('/content/high_male_data.csv')
descriptors = ['立ち幅跳び', 'ハンドボール投げ', '握力得点','上体起こし得点']
model = smf.ols('X50m走 ~ ' + ' + '.join(descriptors), data = high_male)
results = model.fit()
print(results.summary())
best_aic = results.aic
best_model = model
dict_models = {}
while descriptors:
desc_selected = ''
flag = 0
#dict_fitsに辞書keys:削除対象変数 values:[モデル値,AIC]
for rm_desk in descriptors:
used_desks = descriptors.copy()
used_desks.remove(desk)
formula = 'X50m走 ~ ' + ' + '.join(used_desks)
resultmodel = smf.ols(formula = formula, data = high_male)
dict_models[rm_desk] = [resultmodel, resultmodel.fit().aic]
#AICが最小になる
min_k, min_v = min(dict_models.items(), key=lambda x: x[1][1])
if min_v[1] < best_aic:
best_model = min_v[0]
best_aic = min_v[1]
descriptors.remove(min_k)
else:
#削減してもAICの改善が行われなかったら終了
break
stepwise_model_fit = best_model.fit()
print(stepwise_model_fit.summary())
pythonによる出力結果
OLS Regression Results
==============================================================================
Dep. Variable: X50m走 R-squared: 0.525
Model: OLS Adj. R-squared: 0.510
Method: Least Squares F-statistic: 36.18
Date: Sat, 25 Jul 2020 Prob (F-statistic): 2.39e-20
Time: 02:25:36 Log-Likelihood: -41.639
No. Observations: 136 AIC: 93.28
Df Residuals: 131 BIC: 107.8
Df Model: 4
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 10.8194 0.325 33.331 0.000 10.177 11.462
立ち幅跳び -0.0120 0.002 -7.648 0.000 -0.015 -0.009
ハンドボール投げ -0.0144 0.006 -2.367 0.019 -0.026 -0.002
握力得点 -0.0402 0.024 -1.677 0.096 -0.088 0.007
上体起こし得点 -0.0255 0.020 -1.264 0.208 -0.065 0.014
==============================================================================
Omnibus: 3.064 Durbin-Watson: 1.813
Prob(Omnibus): 0.216 Jarque-Bera (JB): 2.955
Skew: 0.359 Prob(JB): 0.228
Kurtosis: 2.927 Cond. No. 2.61e+03
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.61e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
OLS Regression Results
==============================================================================
Dep. Variable: X50m走 R-squared: 0.519
Model: OLS Adj. R-squared: 0.508
Method: Least Squares F-statistic: 47.50
Date: Sat, 25 Jul 2020 Prob (F-statistic): 6.92e-21
Time: 02:25:36 Log-Likelihood: -42.463
No. Observations: 136 AIC: 92.93
Df Residuals: 132 BIC: 104.6
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 10.7121 0.314 34.114 0.000 10.091 11.333
立ち幅跳び -0.0121 0.002 -7.703 0.000 -0.015 -0.009
ハンドボール投げ -0.0169 0.006 -2.929 0.004 -0.028 -0.005
握力得点 -0.0439 0.024 -1.841 0.068 -0.091 0.003
==============================================================================
Omnibus: 2.668 Durbin-Watson: 1.820
Prob(Omnibus): 0.263 Jarque-Bera (JB): 2.632
Skew: 0.334 Prob(JB): 0.268
Kurtosis: 2.867 Cond. No. 2.51e+03
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.51e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
コメント
立ち幅跳び , ハンドボール投げ , 握力得点 , 上体起こし得点はステップワイズ法によって、立ち幅跳び , ハンドボール投げ , 握力得点を使ったモデルが選択されたようです。おそらく、AIC(赤池の情報規準)もこの組合せのときが、一番小さくなると考えられます。
AICの値がRの値とpythonの値で異なっているのは、気になりますがpythonのものでも正しく機能しているように見えます。
pythonで実装したステップワイズ法でも、Rのものと同じものを選択できているようでした。
今後は、ほかのモデル選択法等も検討してみるといいかもしれません。
[参考]Rのソースコード(教材より)
ステップワイズ法（変数増減法）を使って、適切なモデルを選択していています。
# データの読み込み
high_male <- read.csv("high_male_data.csv")
# 全ての説明変数によるステップワイズ法（変数増減法）による重回帰分析
res <- lm(X50m走 ~ 立ち幅跳び + ハンドボール投げ + 握力得点 + 上体起こし得点, data = high_male)
step(res)
[参考]Rによる出力結果
Start: AIC=-292.67
X50m走 ~ 立ち幅跳び + ハンドボール投げ + 握力得点 +
上体起こし得点
Df Sum of Sq RSS AIC
- 上体起こし得点 1 0.1792 14.868 -293.02
<none> 14.689 -292.67
- 握力得点 1 0.3152 15.004 -291.79
- ハンドボール投げ 1 0.6281 15.317 -288.98
- 立ち幅跳び 1 6.5586 21.248 -244.47
Step: AIC=-293.02
X50m走 ~ 立ち幅跳び + ハンドボール投げ + 握力得点
Df Sum of Sq RSS AIC
<none> 14.868 -293.02
- 握力得点 1 0.3817 15.250 -291.58
- ハンドボール投げ 1 0.9663 15.835 -286.46
- 立ち幅跳び 1 6.6831 21.552 -244.54
Call:
lm(formula = X50m走 ~ 立ち幅跳び + ハンドボール投げ +
握力得点, data = high_male)
Coefficients:
(Intercept) 立ち幅跳び ハンドボール投げ 握力得点
10.71205 -0.01210 -0.01689 -0.04389
ソースコード
python版
https://gist.github.com/ereyester/5c6e5a9b8aa55ba826c7c96a4daf7814
R版
https://gist.github.com/ereyester/b42e8d36bf47f70e7092a0d0bd44e673