今回は 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 をインストールしていないのであれば、こちらの記事が参考になるかも。
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) で線形モデルを使った重回帰分析をしてみた。
スマートPythonプログラミング: Pythonのより良い書き方を学ぶ
- 作者: もみじあめ
- 発売日: 2016/03/12
- メディア: Kindle版
- この商品を含むブログを見る
