今回は R と Python の両方を使って重回帰分析をしてみる。
モチベーションとしては、できるだけ手に慣れた Python を使って分析をしていきたいという気持ちがある。
ただ、計算結果が意図通りのものになっているのかを R の結果と見比べて確かめておきたい。
また、分析にはボストンデータセットを用いる。
このデータセットはボストンの各地区ごとの不動産の平均価格と、それに付随する情報が入っている。
今回使った環境は次の通り。
R と Python は、あらかじめインストール済みであることを想定している。
$ sw_vers
ProductName: Mac OS X
ProductVersion: 10.11.6
BuildVersion: 15G1212
$ python --version
Python 3.5.2
$ R --version
R version 3.3.2 (2016-10-31) -- "Sincere Pumpkin Patch"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under the terms of the
GNU General Public License versions 2 or 3.
For more information about these matters see
http://www.gnu.org/licenses/.
もし、まだ R をインストールしていないのであれば、こちらの記事が参考になるかも。
blog.amedama.jp
R
まずは R のやり方から。
最初に R を起動する。
$ R
続いてボストンデータセットの入っているライブラリ MASS を読み込む。
> library(MASS)
データセットの中身は次のようになっている。
> head(Boston)
crim zn indus chas nox rm age dis rad tax ptratio black lstat
1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98
2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14
3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03
4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33
6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21
medv
1 24.0
2 21.6
3 34.7
4 33.4
5 36.2
6 28.7
各変数は次のような意味を持っている。
ただ、今回はボストンデータセットについて詳しく見ることが目的ではないので詳しくは立ち入らない。
1. CRIM: per capita crime rate by town
2. ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
3. INDUS: proportion of non-retail business acres per town
4. CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
5. NOX: nitric oxides concentration (parts per 10 million)
6. RM: average number of rooms per dwelling
7. AGE: proportion of owner-occupied units built prior to 1940
8. DIS: weighted distances to five Boston employment centres
9. RAD: index of accessibility to radial highways
10. TAX: full-value property-tax rate per $10,000
11. PTRATIO: pupil-teacher ratio by town
12. B: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
13. LSTAT: % lower status of the population
14. MEDV: Median value of owner-occupied homes in $1000's
(https://archive.ics.uci.edu/ml/datasets/Housing より)
今回はシンプルに線形回帰モデルを使うので lm() 関数を使う。
medv は目的変数である「その地区の不動産の平均価格」になっている。
次の処理では、目的関数をその他の全ての変数を使って重回帰している。
> lmfit = lm(medv ~ ., data=Boston)
重回帰した結果の詳細は summary() 関数で取得できる。
> summary(lmfit)
Call:
lm(formula = medv ~ ., data = Boston)
Residuals:
Min 1Q Median 3Q Max
-15.595 -2.730 -0.518 1.777 26.199
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
crim -1.080e-01 3.286e-02 -3.287 0.001087 **
zn 4.642e-02 1.373e-02 3.382 0.000778 ***
indus 2.056e-02 6.150e-02 0.334 0.738288
chas 2.687e+00 8.616e-01 3.118 0.001925 **
nox -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
rm 3.810e+00 4.179e-01 9.116 < 2e-16 ***
age 6.922e-04 1.321e-02 0.052 0.958229
dis -1.476e+00 1.995e-01 -7.398 6.01e-13 ***
rad 3.060e-01 6.635e-02 4.613 5.07e-06 ***
tax -1.233e-02 3.760e-03 -3.280 0.001112 **
ptratio -9.527e-01 1.308e-01 -7.283 1.31e-12 ***
black 9.312e-03 2.686e-03 3.467 0.000573 ***
lstat -5.248e-01 5.072e-02 -10.347 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16
自由度調整済み決定係数が 0.7338 ということで、あまり上手く説明できていなさそうな感じ。
Python
続いて Python でのやり方について。
今回は statsmodels というライブラリを使うことにする。
ちなみに、線形回帰モデルは別のライブラリにも色々な実装がある。
ただ、フィッティングしたときに得られるパラメータが statsmodels は細かく取れて R に近いので、今回はこちらを使ってみた。
まずは必要なライブラリをインストールする。
今回は回帰に使わない scikit-learn を入れているのはボストンデータセットを読み込むため。
$ pip install statsmodels scikit-learn
最初に Python の REPL を起動しよう。
$ python
起動したら scikit-learn からボストンデータセットを読み込む。
>>> from sklearn import datasets
>>> dataset = datasets.load_boston()
このデータセットには次のようなメンバがある。
>>> dir(dataset)
['DESCR', 'data', 'feature_names', 'target']
target メンバが目的変数となる各地区ごとの不動産物件の平均価格となっている。
>>> dataset.target[:10]
array([ 24. , 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9])
そして説明変数が data に入っている。
ボストンを 506 地区に分割して 13 の変数を与えたデータになっている。
>>> dataset.data.shape
(506, 13)
変数名は feature_names から取得できる。
>>> dataset.feature_names
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT'],
dtype='<U7')
データセットを読み込めたところで statsmodels をインポートしよう。
>>> from statsmodels import api as sm
線形回帰モデルを使うときは OLS クラスを使う。
第一引数に目的変数を、第二変数に説明変数を渡す。
>>> model = sm.OLS(dataset.target, dataset.data)
準備ができたら fit() メソッドを使って分析する。
>>> result = model.fit()
分析結果からは summary() メソッドを使って R の summary() 関数に近い内容が取得できる。
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.959
Model: OLS Adj. R-squared: 0.958
Method: Least Squares F-statistic: 891.1
Date: Thu, 22 Dec 2016 Prob (F-statistic): 0.00
Time: 22:59:04 Log-Likelihood: -1523.8
No. Observations: 506 AIC: 3074.
Df Residuals: 493 BIC: 3129.
Df Model: 13
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1 -0.0916 0.034 -2.675 0.008 -0.159 -0.024
x2 0.0487 0.014 3.379 0.001 0.020 0.077
x3 -0.0038 0.064 -0.059 0.953 -0.130 0.123
x4 2.8564 0.904 3.160 0.002 1.080 4.633
x5 -2.8808 3.359 -0.858 0.392 -9.481 3.720
x6 5.9252 0.309 19.168 0.000 5.318 6.533
x7 -0.0072 0.014 -0.523 0.601 -0.034 0.020
x8 -0.9680 0.196 -4.947 0.000 -1.352 -0.584
x9 0.1704 0.067 2.554 0.011 0.039 0.302
x10 -0.0094 0.004 -2.393 0.017 -0.017 -0.002
x11 -0.3924 0.110 -3.571 0.000 -0.608 -0.177
x12 0.0150 0.003 5.561 0.000 0.010 0.020
x13 -0.4170 0.051 -8.214 0.000 -0.517 -0.317
==============================================================================
Omnibus: 204.050 Durbin-Watson: 0.999
Prob(Omnibus): 0.000 Jarque-Bera (JB): 1372.527
Skew: 1.609 Prob(JB): 9.11e-299
Kurtosis: 10.399 Cond. No. 8.50e+03
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.5e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
"""
めでたしめでたし、と行きたいところだけどちょっと待ってほしい。
先ほど取得した R の結果と内容が異なっている。
実は R の分析では説明変数にバイアス項という変数が自動で挿入されていた。
これは、不動産物件のベースとなる価格と捉えることができる。
statsmodels を使っているときにバイアス項を説明変数に挿入するときは statsmodels.add_constant() 関数を使う。
バイアス項を挿入した説明変数を X という名前で保存しておこう。
ついでに目的変数も Y という変数で保存しておく。
>>> X = sm.add_constant(dataset.data)
>>> Y = dataset.target
バイアス項を追加した説明変数を用いて、先ほどと同じように回帰してみよう。
>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.741
Model: OLS Adj. R-squared: 0.734
Method: Least Squares F-statistic: 108.1
Date: Thu, 22 Dec 2016 Prob (F-statistic): 6.95e-135
Time: 23:00:34 Log-Likelihood: -1498.8
No. Observations: 506 AIC: 3026.
Df Residuals: 492 BIC: 3085.
Df Model: 13
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const 36.4911 5.104 7.149 0.000 26.462 46.520
x1 -0.1072 0.033 -3.276 0.001 -0.171 -0.043
x2 0.0464 0.014 3.380 0.001 0.019 0.073
x3 0.0209 0.061 0.339 0.735 -0.100 0.142
x4 2.6886 0.862 3.120 0.002 0.996 4.381
x5 -17.7958 3.821 -4.658 0.000 -25.302 -10.289
x6 3.8048 0.418 9.102 0.000 2.983 4.626
x7 0.0008 0.013 0.057 0.955 -0.025 0.027
x8 -1.4758 0.199 -7.398 0.000 -1.868 -1.084
x9 0.3057 0.066 4.608 0.000 0.175 0.436
x10 -0.0123 0.004 -3.278 0.001 -0.020 -0.005
x11 -0.9535 0.131 -7.287 0.000 -1.211 -0.696
x12 0.0094 0.003 3.500 0.001 0.004 0.015
x13 -0.5255 0.051 -10.366 0.000 -0.625 -0.426
==============================================================================
Omnibus: 178.029 Durbin-Watson: 1.078
Prob(Omnibus): 0.000 Jarque-Bera (JB): 782.015
Skew: 1.521 Prob(JB): 1.54e-170
Kurtosis: 8.276 Cond. No. 1.51e+04
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.51e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
"""
今度は R で得られた内容と一致している！
AIC
重回帰分析では、モデリングに用いる説明変数の選択基準として AIC という数値を小さくするやり方があるので、その計算方法について。
R
R では AIC() 関数が用意されているので、それに回帰結果を渡す。
> AIC(lmfit)
[1] 3027.609
Python
Python では statsmodels の回帰結果のインスタンスに aic というメンバがあって、そこに格納されている。
ちなみに AIC はサマリーの中にも表示されていた。
>>> result.aic
3025.6767200074596
標準化
標準化というのはデータセットを平均が 0 で標準偏差が 1 になるように加工すること。
説明変数が標準化されていると、それぞれの説明変数が目的変数に対してどれくらい影響を与えているかを見たりするのに分かりやすい。
R
R には標準化のための関数として scale() が用意されている。
> X <- scale(Boston[, -14])
> Y <- Boston[, 14]
説明変数を標準化して回帰してみよう。
> lmfit <- lm(Y ~ X)
> summary(lmfit)
Call:
lm(formula = Y ~ X)
Residuals:
Min 1Q Median 3Q Max
-15.595 -2.730 -0.518 1.777 26.199
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.53281 0.21095 106.814 < 2e-16 ***
Xcrim -0.92906 0.28269 -3.287 0.001087 **
Xzn 1.08264 0.32016 3.382 0.000778 ***
Xindus 0.14104 0.42188 0.334 0.738288
Xchas 0.68241 0.21884 3.118 0.001925 **
Xnox -2.05875 0.44262 -4.651 4.25e-06 ***
Xrm 2.67688 0.29364 9.116 < 2e-16 ***
Xage 0.01949 0.37184 0.052 0.958229
Xdis -3.10712 0.41999 -7.398 6.01e-13 ***
Xrad 2.66485 0.57770 4.613 5.07e-06 ***
Xtax -2.07884 0.63379 -3.280 0.001112 **
Xptratio -2.06265 0.28323 -7.283 1.31e-12 ***
Xblack 0.85011 0.24521 3.467 0.000573 ***
Xlstat -3.74733 0.36216 -10.347 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16
説明変数が標準化されたことで、どの変数がどれくらい影響を与えているのかが分かりやすくなった。
例えば lstat などは影響が大きそうだ。
Python
Python の場合は scikit-learn の preprocessing() 関数を使うと簡単に標準化できる。
>>> X = dataset.data
>>> from sklearn import preprocessing
>>> X_scaled = preprocessing.scale(X)
標準化した上でバイアス項を追加して回帰してみよう。
>>> X_scaled_with_constant = sm.add_constant(X_scaled)
>>> model = sm.OLS(Y, X_scaled_with_constant)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.741
Model: OLS Adj. R-squared: 0.734
Method: Least Squares F-statistic: 108.1
Date: Thu, 22 Dec 2016 Prob (F-statistic): 6.95e-135
Time: 23:50:36 Log-Likelihood: -1498.8
No. Observations: 506 AIC: 3026.
Df Residuals: 492 BIC: 3085.
Df Model: 13
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const 22.5328 0.211 106.807 0.000 22.118 22.947
x1 -0.9204 0.281 -3.276 0.001 -1.472 -0.368
x2 1.0810 0.320 3.380 0.001 0.453 1.709
x3 0.1430 0.421 0.339 0.735 -0.685 0.971
x4 0.6822 0.219 3.120 0.002 0.253 1.112
x5 -2.0601 0.442 -4.658 0.000 -2.929 -1.191
x6 2.6706 0.293 9.102 0.000 2.094 3.247
x7 0.0211 0.371 0.057 0.955 -0.709 0.751
x8 -3.1044 0.420 -7.398 0.000 -3.929 -2.280
x9 2.6588 0.577 4.608 0.000 1.525 3.792
x10 -2.0759 0.633 -3.278 0.001 -3.320 -0.832
x11 -2.0622 0.283 -7.287 0.000 -2.618 -1.506
x12 0.8566 0.245 3.500 0.001 0.376 1.338
x13 -3.7487 0.362 -10.366 0.000 -4.459 -3.038
==============================================================================
Omnibus: 178.029 Durbin-Watson: 1.078
Prob(Omnibus): 0.000 Jarque-Bera (JB): 782.015
Skew: 1.521 Prob(JB): 1.54e-170
Kurtosis: 8.276 Cond. No. 9.82
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
"""
結果が R の内容と一致していることが分かる。
ちなみに標準化の操作自体は単純なので、自前でやることもできる。
数式に直すと、こんな感じ。
各要素から平均 (
) を引いて標準偏差 (
) で割る。
コードにするとこんな感じ。
>>> X_scaled = (X - X.mean(axis=0)) / X.std(axis=0)
交互作用モデル
交互作用モデルというのは、複数の変数をまとめて説明変数として扱うやり方。
ある説明変数と別の説明変数の相乗作用があるときに使うと良い。
R
R では lm() 関数を実行するときに説明変数を ^2 するだけで交互作用モデルを使える。
> lmfit = lm(medv ~ .^2, data=Boston)
> summary(lmfit)
Call:
lm(formula = medv ~ .^2, data = Boston)
Residuals:
Min 1Q Median 3Q Max
-7.9374 -1.5344 -0.1068 1.2973 17.8500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.579e+02 6.800e+01 -2.323 0.020683 *
crim -1.707e+01 6.554e+00 -2.605 0.009526 **
zn -7.529e-02 4.580e-01 -0.164 0.869508
indus -2.819e+00 1.696e+00 -1.663 0.097111 .
chas 4.451e+01 1.952e+01 2.280 0.023123 *
nox 2.006e+01 7.516e+01 0.267 0.789717
rm 2.527e+01 5.699e+00 4.435 1.18e-05 ***
age 1.263e+00 2.728e-01 4.630 4.90e-06 ***
dis -1.698e+00 4.604e+00 -0.369 0.712395
rad 1.861e+00 2.464e+00 0.755 0.450532
tax 3.670e-02 1.440e-01 0.255 0.798978
ptratio 2.725e+00 2.850e+00 0.956 0.339567
black 9.942e-02 7.468e-02 1.331 0.183833
lstat 1.656e+00 8.533e-01 1.940 0.053032 .
crim:zn 4.144e-01 1.804e-01 2.297 0.022128 *
crim:indus -4.693e-02 4.480e-01 -0.105 0.916621
crim:chas 2.428e+00 5.710e-01 4.251 2.63e-05 ***
crim:nox -1.108e+00 9.285e-01 -1.193 0.233425
crim:rm 2.163e-01 4.907e-02 4.409 1.33e-05 ***
crim:age -3.083e-03 3.781e-03 -0.815 0.415315
crim:dis -1.903e-01 1.060e-01 -1.795 0.073307 .
crim:rad -6.584e-01 5.815e-01 -1.132 0.258198
crim:tax 3.479e-02 4.287e-02 0.812 0.417453
crim:ptratio 4.915e-01 3.328e-01 1.477 0.140476
crim:black -4.612e-04 1.793e-04 -2.572 0.010451 *
crim:lstat 2.964e-02 6.544e-03 4.530 7.72e-06 ***
...(省略)...
rad:tax 3.131e-05 1.446e-03 0.022 0.982729
rad:ptratio -4.379e-02 8.392e-02 -0.522 0.602121
rad:black -4.362e-04 2.518e-03 -0.173 0.862561
rad:lstat -2.529e-02 1.816e-02 -1.392 0.164530
tax:ptratio 7.854e-03 2.504e-03 3.137 0.001830 **
tax:black -4.785e-07 1.999e-04 -0.002 0.998091
tax:lstat -1.403e-03 1.208e-03 -1.162 0.245940
ptratio:black 1.203e-03 3.361e-03 0.358 0.720508
ptratio:lstat 3.901e-03 2.985e-02 0.131 0.896068
black:lstat -6.118e-04 4.157e-04 -1.472 0.141837
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.852 on 414 degrees of freedom
Multiple R-squared: 0.9212, Adjusted R-squared: 0.9039
F-statistic: 53.18 on 91 and 414 DF, p-value: < 2.2e-16
結果から説明変数の数が一気に増えてることが分かる。
これは 13 種類の説明変数から 2 種類を組み合わせで取り出して、それらを交互作用モデルの説明変数として追加しているため。
また、自由度調整済み決定係数が 0.9039 となっており、ただ 13 種類の変数を使うときよりも上手くモデリングできていることが分かる。
まあ、この値だけを見ていても過適合の恐れがあるけど。
Python
statsmodels には回帰するときのやり方を R のような書き方 (R-style らしい) で指定できる。
このやり方を使うとさくっと相互作用モデルも使えそうな感じ。
Fitting models using R-style formulas — statsmodels 0.7.0 documentation
ただまあ、これを使うよりも自前で用意する方が楽しそう。
なので、今回は自分で計算してみることにする。
Python で組み合わせを計算するときは itertools モジュールを使うと良い。
この中にある combinations() 関数を使えば簡単に組み合わせが取れる。
例えば説明変数の組み合わせを計算したいなら、こう。
>>> import itertools
>>> list(itertools.combinations(dataset.feature_names, 2))
[('CRIM', 'ZN'), ('CRIM', 'INDUS'), ('CRIM', 'CHAS'), ('CRIM', 'NOX'), ('CRIM', 'RM'), ('CRIM', 'AGE'), ('CRIM', 'DIS'), ('CRIM', 'RAD'), ('CRIM', 'TAX'), ('CRIM', 'PTRATIO'), ('CRIM', 'B'), ('CRIM', 'LSTAT'), ('ZN', 'INDUS'), ('ZN', 'CHAS'), ('ZN', 'NOX'), ('ZN', 'RM'), ('ZN', 'AGE'), ('ZN', 'DIS'), ('ZN', 'RAD'), ('ZN', 'TAX'), ('ZN', 'PTRATIO'), ('ZN', 'B'), ('ZN', 'LSTAT'), ('INDUS', 'CHAS'), ('INDUS', 'NOX'), ('INDUS', 'RM'), ('INDUS', 'AGE'), ('INDUS', 'DIS'), ('INDUS', 'RAD'), ('INDUS', 'TAX'), ('INDUS', 'PTRATIO'), ('INDUS', 'B'), ('INDUS', 'LSTAT'), ('CHAS', 'NOX'), ('CHAS', 'RM'), ('CHAS', 'AGE'), ('CHAS', 'DIS'), ('CHAS', 'RAD'), ('CHAS', 'TAX'), ('CHAS', 'PTRATIO'), ('CHAS', 'B'), ('CHAS', 'LSTAT'), ('NOX', 'RM'), ('NOX', 'AGE'), ('NOX', 'DIS'), ('NOX', 'RAD'), ('NOX', 'TAX'), ('NOX', 'PTRATIO'), ('NOX', 'B'), ('NOX', 'LSTAT'), ('RM', 'AGE'), ('RM', 'DIS'), ('RM', 'RAD'), ('RM', 'TAX'), ('RM', 'PTRATIO'), ('RM', 'B'), ('RM', 'LSTAT'), ('AGE', 'DIS'), ('AGE', 'RAD'), ('AGE', 'TAX'), ('AGE', 'PTRATIO'), ('AGE', 'B'), ('AGE', 'LSTAT'), ('DIS', 'RAD'), ('DIS', 'TAX'), ('DIS', 'PTRATIO'), ('DIS', 'B'), ('DIS', 'LSTAT'), ('RAD', 'TAX'), ('RAD', 'PTRATIO'), ('RAD', 'B'), ('RAD', 'LSTAT'), ('TAX', 'PTRATIO'), ('TAX', 'B'), ('TAX', 'LSTAT'), ('PTRATIO', 'B'), ('PTRATIO', 'LSTAT'), ('B', 'LSTAT')]
組み合わせは 
通り。
>>> len(list(itertools.combinations(dataset.feature_names, 2)))
78
それぞれの説明変数を取り出すための添字にするとこんな感じ。
>>> list(itertools.combinations(range(dataset.data.shape[1]), 2))
[(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10), (0, 11), (0, 12), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9), (1, 10), (1, 11), (1, 12), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10), (2, 11), (2, 12), (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10), (3, 11), (3, 12), (4, 5), (4, 6), (4, 7), (4, 8), (4, 9), (4, 10), (4, 11), (4, 12), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10), (5, 11), (5, 12), (6, 7), (6, 8), (6, 9), (6, 10), (6, 11), (6, 12), (7, 8), (7, 9), (7, 10), (7, 11), (7, 12), (8, 9), (8, 10), (8, 11), (8, 12), (9, 10), (9, 11), (9, 12), (10, 11), (10, 12), (11, 12)]
交互作用モデルで使うそれぞれの説明変数は、説明変数同士を乗算することで計算できる。
>>> b = np.array([X[:, 0] * X[:, 1]])
計算した説明変数を、既存の説明変数にカラムとして追加するには次のようにすれば良い。
>>> a = X
>>> np.concatenate((a, b.T), axis=1)
array([[ 6.32000000e-03, 1.80000000e+01, 2.31000000e+00, ...,
3.96900000e+02, 4.98000000e+00, 1.13760000e-01],
[ 2.73100000e-02, 0.00000000e+00, 7.07000000e+00, ...,
3.96900000e+02, 9.14000000e+00, 0.00000000e+00],
[ 2.72900000e-02, 0.00000000e+00, 7.07000000e+00, ...,
3.92830000e+02, 4.03000000e+00, 0.00000000e+00],
...,
[ 6.07600000e-02, 0.00000000e+00, 1.19300000e+01, ...,
3.96900000e+02, 5.64000000e+00, 0.00000000e+00],
[ 1.09590000e-01, 0.00000000e+00, 1.19300000e+01, ...,
3.93450000e+02, 6.48000000e+00, 0.00000000e+00],
[ 4.74100000e-02, 0.00000000e+00, 1.19300000e+01, ...,
3.96900000e+02, 7.88000000e+00, 0.00000000e+00]])
これで説明変数が一つ増える。
>>> np.concatenate((a, b.T), axis=1).shape
(506, 14)
上記にもとづいて相互作用モデルで使う全ての変数を揃える。
>>> X = dataset.data
>>> for i, j in itertools.combinations(range(X.shape[1]), 2):
... interaction_column = np.array([X[:, i] * X[:, j]])
... X = np.concatenate((X, interaction_column.T), axis=1)
...
>>> X.shape
(506, 91)
そして、バイアス項を追加する。
>>> X = sm.add_constant(X)
>>> X.shape
(506, 92)
上記の説明変数を使って回帰する。
>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.921
Model: OLS Adj. R-squared: 0.904
Method: Least Squares F-statistic: 53.21
Date: Thu, 22 Dec 2016 Prob (F-statistic): 5.85e-181
Time: 19:16:55 Log-Likelihood: -1197.3
No. Observations: 506 AIC: 2579.
Df Residuals: 414 BIC: 2967.
Df Model: 91
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const -159.4473 67.980 -2.346 0.019 -293.076 -25.819
x1 -16.9621 6.549 -2.590 0.010 -29.835 -4.089
x2 -0.0769 0.458 -0.168 0.867 -0.977 0.823
x3 -2.8163 1.695 -1.662 0.097 -6.148 0.516
x4 44.6397 19.519 2.287 0.023 6.271 83.008
x5 22.1561 75.150 0.295 0.768 -125.567 169.879
x6 25.4191 5.700 4.459 0.000 14.214 36.624
x7 1.2569 0.273 4.612 0.000 0.721 1.793
x8 -1.6398 4.604 -0.356 0.722 -10.690 7.410
x9 1.8056 2.462 0.733 0.464 -3.034 6.646
x10 0.0374 0.144 0.260 0.795 -0.246 0.320
x11 2.7520 2.849 0.966 0.335 -2.849 8.353
x12 0.1008 0.074 1.354 0.177 -0.046 0.247
x13 1.6587 0.853 1.944 0.053 -0.018 3.336
x14 0.4124 0.180 2.287 0.023 0.058 0.767
x15 -0.0467 0.448 -0.104 0.917 -0.927 0.834
x16 2.4377 0.571 4.271 0.000 1.316 3.560
x17 -1.2206 0.920 -1.327 0.185 -3.029 0.588
x18 0.2117 0.049 4.345 0.000 0.116 0.307
x19 -0.0031 0.004 -0.828 0.408 -0.011 0.004
x20 -0.1956 0.106 -1.848 0.065 -0.404 0.012
x21 -0.6599 0.581 -1.135 0.257 -1.803 0.483
x22 0.0349 0.043 0.815 0.416 -0.049 0.119
x23 0.4895 0.333 1.471 0.142 -0.164 1.143
x24 -0.0004 0.000 -2.513 0.012 -0.001 -9.71e-05
x25 0.0295 0.007 4.512 0.000 0.017 0.042
...(省略)...
x82 3.151e-05 0.001 0.022 0.983 -0.003 0.003
x83 -0.0439 0.084 -0.523 0.601 -0.209 0.121
x84 -0.0004 0.003 -0.165 0.869 -0.005 0.005
x85 -0.0251 0.018 -1.384 0.167 -0.061 0.011
x86 0.0078 0.003 3.135 0.002 0.003 0.013
x87 -2.303e-06 0.000 -0.012 0.991 -0.000 0.000
x88 -0.0014 0.001 -1.167 0.244 -0.004 0.001
x89 0.0012 0.003 0.349 0.727 -0.005 0.008
x90 0.0038 0.030 0.128 0.898 -0.055 0.062
x91 -0.0006 0.000 -1.510 0.132 -0.001 0.000
==============================================================================
Omnibus: 121.344 Durbin-Watson: 1.524
Prob(Omnibus): 0.000 Jarque-Bera (JB): 628.079
Skew: 0.942 Prob(JB): 4.11e-137
Kurtosis: 8.123 Cond. No. 1.18e+08
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.18e+08. This might indicate that there are
strong multicollinearity or other numerical problems.
"""
上記を見ると R で相互作用モデルを使って計算したときの内容と一致していることが分かる。
まとめ
今回は R と Python (statsmodels) で線形モデルを使った重回帰分析をしてみた。
![日本統計学会公式認定 統計検定 2級 公式問題集[2013〜2015年]](https://images-fe.ssl-images-amazon.com/images/I/51PGSetFISL._SL160_.jpg "日本統計学会公式認定 統計検定 2級 公式問題集[2013〜2015年]")
