重回帰

サンプル:bostonデータ

統計ソフトの多くに付属しているボストンデータを例に重回帰について試してみる。

ボストンデータとは、ボストン市内の住宅物件の価格と様々な社会的データをまとめたものであり、それぞれの社会的データから価格を予測(評価)する手法を考察しようとするものである。

データの内容

以下のような項目(Attribute)と物件の価格のデータが500件あまり用意されている。

Attribute description
CRIM 人口 1 人当たりの犯罪発生数
ZN 25,000平方フィート以上の住居区画の占める割合
INDUS 小売業以外の商業が占める面積の割合
CHAS チャールズ川によるダミー変数 (1: 川の周辺, 0: それ以外)
NOX NOx の濃度
RM 住居の平均部屋数
AGE 1940 年より前に建てられた物件の割合
DIS 5 つのボストン市の雇用施設からの距離 (重み付け済)
RAD 環状高速道路へのアクセスしやすさ
TAX $10,000 ドルあたりの不動産税率の総計
PTRATIO 町毎の児童と教師の比率
B 町毎の黒人 (Bk) の比率を次の式で表したもの。 1000(Bk – 0.63)2
LSTAT 給与の低い職業に従事する人口の割合 (%)

これらの項目を説明変数(feature 特徴量)、住宅物件の価格を目的変数(target)として、両者が線形な関係にあると仮定して、その係数を求めることを線形重回帰と呼ぶ。

説明変数を$(x_1, x_2, \cdots, x_M)$、目的変数を$y$で表せば、

$$ y = w_1 x_1 + w_2 x_2 + \cdots + w_M x_M $$

参考:"scikit-learn boston 重回帰"くらいで検索すると山のようにページがでてくるが、上位のものとしては

そのほか、ランダムフォレストを使った分析もある。(後の回にふれるかも)

重回帰は、複数の 以下では、最初にあげたページに沿って実行(訓練)し、係数を決めた後、テストデータを使って予測の良さを調べてみる。

データの準備(読み込み)

In [2]:
from sklearn.datasets import load_boston
boston = load_boston() # データセットの読み込み
In [3]:
# Bostonデータの項目確認
print(boston.DESCR)
.. _boston_dataset:

Boston house prices dataset
---------------------------

**Data Set Characteristics:**  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

    :Attribute Information (in order):
        - CRIM     per capita crime rate by town
        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
        - INDUS    proportion of non-retail business acres per town
        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
        - NOX      nitric oxides concentration (parts per 10 million)
        - RM       average number of rooms per dwelling
        - AGE      proportion of owner-occupied units built prior to 1940
        - DIS      weighted distances to five Boston employment centres
        - RAD      index of accessibility to radial highways
        - TAX      full-value property-tax rate per $10,000
        - PTRATIO  pupil-teacher ratio by town
        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
        - LSTAT    % lower status of the population
        - MEDV     Median value of owner-occupied homes in $1000's

    :Missing Attribute Values: None

    :Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/


This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic
prices and the demand for clean air', J. Environ. Economics & Management,
vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics
...', Wiley, 1980.   N.B. Various transformations are used in the table on
pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression
problems.   
     
.. topic:: References

   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.
   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.

データを訓練データとテスト用データに分割

以下で、randome_state=0として分割の仕方を固定している。実際には、何種類もの分割を試して、モデルの良さを確かめる、交差検定を行う。

In [4]:
from sklearn.model_selection import train_test_split
# 入出力の切り分け
x = boston['data']  # 物件の情報
t = boston['target']  # 家賃
x_train, x_test, t_train, t_test = train_test_split(x, t, test_size=0.3, random_state=0) # データの70%を訓練用に、残りをテスト用に

print(x_train[:3,:]) # 試しに最初の3行を出力
[[1.62864e+00 0.00000e+00 2.18900e+01 0.00000e+00 6.24000e-01 5.01900e+00
  1.00000e+02 1.43940e+00 4.00000e+00 4.37000e+02 2.12000e+01 3.96900e+02
  3.44100e+01]
 [1.14600e-01 2.00000e+01 6.96000e+00 0.00000e+00 4.64000e-01 6.53800e+00
  5.87000e+01 3.91750e+00 3.00000e+00 2.23000e+02 1.86000e+01 3.94960e+02
  7.73000e+00]
 [5.57780e-01 0.00000e+00 2.18900e+01 0.00000e+00 6.24000e-01 6.33500e+00
  9.82000e+01 2.11070e+00 4.00000e+00 4.37000e+02 2.12000e+01 3.94670e+02
  1.69600e+01]]

訓練(scikit-learnの重回帰モジュールの適用)

順序は

  1. moduleの読み込み
  2. objectの用意 (下のプログラムではmodel)
  3. メソッドfitに説明変数、目的変数を与えて学習 (係数はmodelの中に保持)
    • 説明変数は、列として各種項目(feature)を、1データ1行の行列として与える。

その後、

  1. テスト用データとしてとっておいた、説明変数値を与えて予測値を算出し、テスト用の目的変数(住宅価格)と比べてみる。
In [5]:
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(x_train, t_train)

model.score(x_train, t_train) # 決定係数R^2の出力
Out[5]:
0.7645451026942549

決定係数

テストデータの実測値を$(t_1, t_2, \cdots, t_N)$、予測(推定)値を$(y_1, y_2, \cdots, y_N)$としたとき $$ R^2 = 1 - \frac{\sum_1^N (y_i - t_i)^2}{ \sum_1^N (t_i - \bar{t})^2} $$

In [6]:
# テストデータの予測値と実測値比較

x_predict = model.predict(x_test) # テストデータ(物件)の予測値

import matplotlib.pyplot as plt
plt.plot(x_predict, t_test, "o")
plt.xlabel("予測値")
plt.ylabel("実測値")
Out[6]:
Text(0, 0.5, '実測値')

今後の授業内容との関係

授業では、もっぱら説明変数が1つで目的変数との関係が非線形な場合を扱う。

ここで線形重回帰に触れたのは、重回帰の場合

$$ y = w_1 x_1 + w_2 x_2 + \cdots + w_M x_M $$

で使ったLinearRegressionは非線形な回帰で用いることができるという点である。

一つの説明変数$x$と目的変数$y$が非線形な関係にあると推量される場合、、例えばその関係が「べき関数」だとすれば、各項目(Attribute)を$x$のべき乗として $$ y = w_1 x + w_2 x^2 + \cdots w_M x^M $$

のようにかけば、係数$\{w_i\}$を線形重回帰で推定することと同じであうことがわかるだろう。

$N$個のデータがあるとき、説明変数としては、各行に$x^1, x^2,\cdots, x^{M-1}$の値をセットした$N$行の行列を与えることになる。

[補足] 一般に、1つの説明変数$x$、1つの目的変数$y$の間に非線形な関係がある場合、$x$に対するいくつかの関数を用意し、それの線形結合によって表す。

$$ y = w_1\phi_1 (x) + w_2 \phi_2 (x) + \cdots + w_M\phi_M (x) $$

のように書いて、線形重回帰と同様のスキームで扱う。このような場合、 $\{\phi_i\}$のセットを「基底」と呼び、データの性質によって、べき以外に、ガウシアンやシグモイド関数($\tanh$関数)、フーリエ基底(三角関数)を用いることが多い。フーリエ基底の場合はwaveletと呼ばれる。

補足 (pandasの利用)

データの可視化や様々な加工にはpandasを用いる方が便利である。

(例1) 結果を見やすく表示

In [7]:
# 回帰係数の出力
print(model.coef_)
[-1.21310401e-01  4.44664254e-02  1.13416945e-02  2.51124642e+00
 -1.62312529e+01  3.85906801e+00 -9.98516565e-03 -1.50026956e+00
  2.42143466e-01 -1.10716124e-02 -1.01775264e+00  6.81446545e-03
 -4.86738066e-01]
In [8]:
# 回帰係数とfeatureの名前を結合してみる
import numpy as np
coef_names = np.stack([boston.feature_names, model.coef_])
print(coef_names)
[['CRIM' 'ZN' 'INDUS' 'CHAS' 'NOX' 'RM' 'AGE' 'DIS' 'RAD' 'TAX' 'PTRATIO'
  'B' 'LSTAT']
 ['-0.12131040096834629' '0.04446642542889598' '0.011341694511898402'
  '2.5112464244913792' '-16.23125290285598' '3.8590680098251915'
  '-0.009985165654562847' '-1.500269563265005' '0.2421434660581058'
  '-0.011071612402085984' '-1.0177526384185978' '0.006814465447642474'
  '-0.48673806564491884']]
In [34]:
# 上記が見にくいので、pandasのDataFrameに直して出力してみる
import pandas as pd
result= pd.DataFrame(coef_names)
# NOXが負の大きな値に、RMが正の比較的大きな値になっていることがわかる。
result
Out[34]:
0 1 2 3 4 5 6 7 8 9 10 11 12
0 CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT
1 -0.12131040096834629 0.04446642542889598 0.011341694511898402 2.5112464244913792 -16.23125290285598 3.8590680098251915 -0.009985165654562847 -1.500269563265005 0.2421434660581058 -0.011071612402085984 -1.0177526384185978 0.006814465447642474 -0.48673806564491884
In [36]:
# 行と列を入れ替え
result = pd.DataFrame(coef_names).T
# 絶対値で降順にソート
result['abs_coef'] = np.abs(result[1].astype(float))
result.sort_values('abs_coef', ascending=False)
Out[36]:
0 1 abs_coef
4 NOX -16.23125290285598 16.231253
5 RM 3.8590680098251915 3.859068
3 CHAS 2.5112464244913792 2.511246
7 DIS -1.500269563265005 1.500270
10 PTRATIO -1.0177526384185978 1.017753
12 LSTAT -0.48673806564491884 0.486738
8 RAD 0.2421434660581058 0.242143
0 CRIM -0.12131040096834629 0.121310
1 ZN 0.04446642542889598 0.044466
2 INDUS 0.011341694511898402 0.011342
9 TAX -0.011071612402085984 0.011072
6 AGE -0.009985165654562847 0.009985
11 B 0.006814465447642474 0.006814

係数の大きさは、単純には寄与度を表すが、注意が必要 ⇒ 一番下の記述を参照!

散布図matrixの描画

その他はWeb上のいろいろなページを参照願いたい。

In [13]:
# 元のデータをpandasのDataFrameに入れて処理
import pandas as pd
boston_df = pd.DataFrame(boston.data, columns=boston.feature_names)# 説明変数(data) 
boston_df['PRICE'] = boston.target # 目的変数(target)をDataFrameの列として追加
boston_df.head()
Out[13]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE
0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0
1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6
2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7
3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4
4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2
In [16]:
import matplotlib.pyplot as plt

# 全データ項目の散布図
pd.plotting.scatter_matrix(boston_df,figsize=(20,20))
plt.plot()
Out[16]:
[]

一番下の散布図に注目してみよう。

縦軸は、目的変数の「価格」である。

RMと相関、NOXとLSTATが逆相関になっていることが目立つだろう。

RMとNOXの重回帰係数の絶対値が大きくなっていることと符合する。が、LSTATの係数が小さいのはなぜだろうか?  ⇒ データの「正規化」あるいは標準化が必要なことを示唆している。

(時間があれば実際に行ってみよう)

ボストンデータのような様々な属性のデータを使って、予測する場合には、木構造にデータを分割するランダムフォレストによる回帰分析が良く用いられるる。 たとえば、アルゴリズムの説明を抜きにシンプルなコードのみを説明したページとして

https://hinomaruc.hatenablog.com/entry/2019/11/14/200857

を上げて置く。

In [ ]: