はじめに

今回は文部科学省のページで公開されている情報Ⅰの教員研修用教材の「第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