scikit-learn機械学習④多変量線形回帰：実践編

今回は、複数の変数から値を予測する多変量線形回帰をscikit-learn を使って実装します。

データセットとしてカルフォルニアの住宅の値段に関する california_housing を使います。これは、scikit-learn に付属するもので手軽に使えるので便利です。このデータには、家の値段とその家の属性に関するデータが入っています。

例えば、築年数や部屋数などの属性が含まれており、これら複数の変数から家の値段を予測することを考えます。仮に、築年数と部屋数だけで予測するとすると、多変量線形回帰は次のような式になります。

$$
家の値段 = a  \times 築年数 + b \times 部屋数 + c + 誤差
$$

もちろん、家の値段が単純に築年数と部屋数だけで決まるものではないですが、ここでは説明のために変数を二つに絞っています。後で実装する時にはもっと多くの変数を扱います。

この式の要点は、複数の変数によって予測する値が変わってくることです。また、それぞれの変数が値段に与える影響度が a や b によって決まります。つまり、a の値が b よりも比較的に大きい場合は、築年数の方が部屋数よりも値段に与える影響が大きいということです。よって、多変量線形回帰のモデルを構築することで属性それぞれの重要度も見えてきます。

多変量線形回帰のモデルがデータを的確に説明できるようにするためには、予測の誤差がなるべく小さくなるようにする必要があります。よって、上式のパラメータである a, b, c を調節して誤差が小さくなるようにするのが多変量線形回帰による機械学習の目的となります。

それでは、データの読み込みや分析から始めて、実際に多変量線形回帰をscikit-learn を使って実装していきましょう。

多変量による予測

多変量線形回帰による家の値段の予測式をもう一度見てみましょう。

$$
家の値段 = a  \times 築年数 + b \times 部屋数 + c + 誤差
$$

上式では、変数の数だけ a や b といった係数（傾き）があります。これは、それぞれの属性を表す独立変数が、家の値段である従属変数（目的変数）に与える影響の大きさを表しています。

あと、c は切片です。これは、理論上は全ての属性が０の時の家の値段となりますが、現実的には築年数が０で部屋数が０の家など売りに出ていないので土地の値段とでも考えても良いのですが、単純にデータに線形回帰モデルをフィットさせるときに必要な値として捉えて問題ありません。

なお、誤差があるのは、このモデルで捉えきれないランダムな値段の変動を表しています。線形回帰では、この誤差をできるだけ小さくなるように傾きと切片を調整します。

線形回帰の理論編で解説しましたが、単純線形回帰と多変量線形回帰は独立変数が一つであるか複数であるかの違いだけで、行列を使って表現すると全く同じ式になります。

まとめると、一つの属性だけで説明できないことを扱うのが多変量線形回帰モデルだと考えて大方間違いありません。もちろん、多変量を扱うモデルは線形回帰だけではありません。実際にはデータを読み込んで、まずは探索的な分析を行います。その上で、どのようなモデルを使うのかを考察するといった流れになります。

ただし、この記事では多変量線形回帰を実装することを目的としているので、そのようなモデル選択のための考察は行いません。

Python環境の設定

Pythonの仮想環境を作ってscikit-learnとJupyterなど必要なライブラリをインストールします。今回は、pandas も使います。

mkdir multlinear_regression
cd multilinear_regression

python3 -m venv venv
source venv/bin/activate

# pip をアップグレードしておく
pip install --upgrade pip

# 必要なライブラリをインストール
pip install scikit-learn jupyter matplotlib pandas

Jupyterノートブックを立ち上げてPython３のノートブックを作成してください。Jupyterノートブックに関しては、こちらを参照してください。

なお、ノートブックの閲覧編集にはVSCode（Visual Studio Code）も使えます。私は、どちらかというとVSCodeをよく使います。これについてもこちらで簡単に解説しています。

データの読み込み

カルフォルニアの住宅データは、scikit-learn から次のように読み込みます。

from sklearn.datasets import fetch_california_housing


housing = fetch_california_housing()

この housing の型を見てみると Bunch となっています。

この Bunch オブジェクトには、以下のような属性や関数があります。

DESCR（description）は、カルフォルニア住宅のデータセットの説明のテキストです。その内容はプリントすれば見ることができます。

説明は続きますが、ここでは主要な部分だけを表示しています。次のような情報があるのがわかります。

• 住宅データの数：20640

• ターゲットの数：1（これは値段のことで、$100,000単位）

• 数値の属性の数：8

  • MedInc：地域の収入の中央値（メジアン、$10,000単位）

  • HouseAge：家の築年数

  • AveRooms：地域の部屋数の平均

  • AveBedrms：地域の寝室数の平均

  • Population：地域の人口数

  • AveOccup：地域の世帯数の平均

  • Latitude：緯度

  • Longitude：経度

仮に、実際にデータの分析を行うことを想定すると、これら全てが家の値段に影響するのか、また変数と値段の関係は、多変量線形回帰である程度説明できるのかといった疑問が湧くところですね。さらに、どの変数がより重要なのかも気になるところです。

よく見てみると住宅に特有の情報としては、築年数、緯度、経度だけで後は地域の平均の情報（部屋数、寝室数など）であるのが分かります。よって、このデータで住宅の値段を予測するとなると、地域全体の住宅の値段の傾向が大きく影響することになりそうです。

さらに、調べるとターゲットである MedHouseVal とは、この住宅のある地域の住宅価格の中央値（メジアン）なので、この家の値段ではありません。よって、このデータを使っても各住宅の値段を予測するのには役に立ちません。むしろ、その地域の平均的な住宅価格を知ることができることになります。

そう考えると緯度経度の情報は案外と地域の住宅価格を予測するのに役立つかもしれません。結局のところ住宅の値段って周りの売買価格に左右されるものなのでしょう。などなど、データを見ていると色々と考えを巡らせることができます。

データの分析

もう少しデータを詳しく見ていきましょう。

なお、線形改軌における独立変数を機械学習では特徴量と呼んだりもします。この住宅データのオブジェクト（Bunch）では、特徴量の名前のリストを次のように確認することもできます。

また、目的変数の名前も次のように得ることができます。

このデータセットでは、目的変数は一つしかないですが、リスト形式で返されています。

次に、データの型を確認します。

これは、NumPyの配列型です。シェイプを見てみます。

先ほど説明のテキストで見たようにデータ数は 20640 で、属性の数は 8 であるのが確認できます。

最初のデータを見てみましょう。

これらの値は、先ほど見た特徴量の名前（housing.feature_names）に対応しています。

• MedInc：地域の収入の中央値（メジアン）、8.3252（83,252 USドル）

• HouseAge：家の築年数、41

• AveRooms：地域の部屋数の平均、6.98

• AveBedrms：地域の寝室数の平均、1.02

• Population：地域の人口数、322

• AveOccup：地域の世帯数の平均、2.56

• Latitude：緯度、37.88

• Longitude：経度、-122.23

次に目的変数のデータ型を見てみましょう。

おそらく予測がついたと思いますが、これもNumPyの配列型になっています。シェイプを見てみます。

20640個の値がある１次元配列（リスト、ベクトル）です。これも予想通りでしょう。一番最初の値を見ます。

つまり、4.526 x 100,000 = 452,600 USドルということでカルフォルニアの住宅は高いですね。 

Pandasによる分析

カリフォルニア住宅データセットを pandas のデータフレームとして読み込む事もできます。pandas はテーブル形式のデータをPythonで分析する場合によく使われるライブラリです。Excel や Google Sheet で扱えるようなデータを扱うのは pandas が得意とするところです。

次のように、データを読み込みます。

housing2 = fetch_california_housing(as_frame=True)

as_frame=True と指定することで、データが pandas のデータフレームとして読み込まれます。

housing2.frame にデータフレームがあります。

これで、属性の値と名前が簡単に見分けられますし、データのシェイプ（20640 x 9）確認できます。９となっているのは属性とターゲット（一番右の列、MedHouseVal）を含んでいるからです。

次のようにすると属性だけを取り出すことができます。

housing2.data

また、ターゲットだけを次のように取り出すことも可能です。

housing2.target

ちなみに、frame と data と target のデータ型は次のようになっています。

type(housing2.frame), type(housing2.data), type(housing2.target)

つまり、frame と data は、pandas.core.frame.DataFrame でデータフレーム型です。それに対して、target は pandas.core.series.Series とシリーズ型になっています。データフレームが表のデータを２次元の行列として扱うのに対し、シリーズは表のデータの１列を取り出したものと考えると分かりやすいでしょう。

pandas の形式でデータを取り出すと、簡単な統計やグラフを素早く作ることができます。

例えば、データフレームの describe() を呼び出してみます。

housing2.frame.describe()

このように列ごとに統計が表示されます。

• 平均（mean）

• 標準偏差（std）

• 最小値（min）

• 25パーセンタイル（25%）

• 50パーセンタイル（50%、メジアン）

• 75パーセンタイル（75%）

• 最大値（max）

これらの値は、データの分布などをみるためのものです。

また、HouseAge：家の築年数の統計だけを見ることもできます。

housing2.frame["HouseAge"].describe()

平均とメジアンがそれぞれ28.6年と29年となっています。また、一番古くても52年です。

また、次のようにデータの分布を図にすることもできます。

housing2.frame.hist(figsize=(12, 10), bins=30, edgecolor="k")

これは、データフレームの各列のデータをヒストグラムにしたものです。

こうしてみると HouseAge が50以上のものが結構あるのが分かります。また、MedHouseVal（地域の住宅価格の中央値）でも最大値のところで頻度が大きくなっています。これは地域によって価格に大きな差があることの表れでしょう。

緯度軽度で住宅価格の分布を見てみましょう。

housing2.frame.plot.scatter(
    x='Longitude', 
    y='Latitude', 
    s='MedHouseVal', 
    label='Median House Value', 
    alpha=0.6, 
    figsize=(12, 10), 
    c='MedHouseVal', 
    cmap='jet', 
    colorbar=True)

うっすらとカルフォルニア州の形が見えます。海岸沿いの住宅の値段が高いようです。特にサンフランシスコとロサンジェルスあたりですかね。

なお、pandas は中で matplotlib を使ってグラフを表示しているので matplotlib が得意な方は直接そちらを使っても同じようなグラフが描けるはずです。

モデルの学習

目的変数である MedHouseVal は、個々の住宅の値段ではなく、その地域の住宅の値段の中央値であることがわかりました。よって、住宅の値段を予測するのではなくなりましたが、多変量線形回帰モデルでこのデータセットにフィットさせることができるかやってみましょう。

ここでもう一つのデータの読み込み方を使います。

X, y = fetch_california_housing(return_X_y=True)

これで独立変数のデータは X に、目的変数のデータは y に格納されます。

X と y のそれぞれのシェイプを見てみます。

両方とも 20640個のデータがあり、X には 8個の変数があるのがわかります。これは期待通りです。

X のデータは NumPy配列になっています。

y のデータもNumPy配列です。

では、X と y を訓練用とテスト用に分けます。前回の記事と同様ですね。

from sklearn.model_selection import train_test_split

# データをトレーニングセットとテストセットに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, 
                                                    random_state=0)

次に線形回帰モデル（LinearRegression）のオブジェクトを生成して、データにフィットさせます。

from sklearn.linear_model import LinearRegression

# 線形回帰モデルのインスタンスを作成
model = LinearRegression()

# モデルをトレーニングデータにフィットさせる
model.fit(X_train, y_train)

このモデルを使ってテストデータを使った予測を行います。

# テストデータで予測
y_pred = model.predict(X_test)

この予測値を使って平均２乗誤差を計算します。

from sklearn.metrics import mean_squared_error

# モデルの性能を評価
mse = mean_squared_error(y_test, y_pred)
print(f'Mean Squared Error: {mse:.2f}')

0.53と出ました。前回同様、この値が良いか悪いかは、そもそものデータの分散などにもよります。MSEがゼロになれば、全ての点が綺麗に線形モデルで説明できることになりますが、実際のデータの関係が線形であるとは限らず、また誤差などもあり、ほとんどの場合はゼロにはなりません。

とはいっても、さらに深くモデルの評価を行うことはできます。

モデルの評価

係数の評価

まず、モデルが選んだ係数の値を見てみましょう。model.coef_ から係数を読むことができます。なお、coef は、coefficient （係数）の略です。

# 線形回帰モデルの係数を表示
coefficients = model.coef_

print(coefficients)

ただし、このままだと読みにくいので、次のように特徴量の名前と一緒に表示させます。

# 特徴量の名前と係数を対応付けて表示
feature_names = fetch_california_housing().feature_names

for feature_name, coefficient in zip(feature_names, coefficients):
    print(f'{feature_name:12}: {coefficient: .2f}')

これらの係数は、線形回帰モデルから得られたものであり、各特徴量が住宅価格の中央値（`MedHouseVal`）に与える影響を示しています。以下は各係数の簡単な分析（感想）です。

• MedInc (0.43): 地域の所得の中央値が高いほど住宅価格が高くなる傾向がある。

• HouseAge (0.01): 住宅の年齢が住宅価格に与える影響は比較的小さい。

• AveRooms (-0.10): 平均部屋数が多いほど住宅価格が低くなる傾向がある。部屋数が多い住宅は部屋が極端に小さいなどの理由があるのでしょうか。これはデータからはわからないでしょう。

• AveBedrms (0.59): 平均ベッドルーム数の値が大きいほど住宅価格が高くなる傾向がある。まあ、これは妥当な感じがします。

• Population (-0.00), AveOccup (-0.00): これらの特徴量の係数は非常に小さい負の値です。これは人口や一世帯当たりの平均居住者数が住宅価格にほとんど影響を与えないことを示しています。

• Latitude (-0.42), Longitude (-0.43): これをみると北方または西方に位置する地区が低い住宅価格と関連しているようです。逆にいうと海沿いや南の方の値段が高いので、やはり、地理的位置が住宅価格に大きな影響を与えるのが見えます。

これらの結果から、地域の中央所得（`MedInc`）と平均寝室数（AveBedrms）と地理的位置（緯度と経度）が住宅価格の中央値に最も強い影響を与えていることがわかります。

R2スコア（決定係数）

モデルの性能評価として、R2スコア（決定係数、R-squared）を使います。R2スコアは、モデルがデータをどれだけ説明できているかを示す指標で、1に近いほどモデルの予測が実際のデータに適合していることを意味します。逆に、0に近い場合は、モデルがデータの変動をほとんど説明できていないことを示します。

次のようにR2スコアを計算します。

from sklearn.metrics import r2_score

# R2スコアを計算
r2 = r2_score(y_test, y_pred)

print(f'R2 Score: {r2:.2f}')

R2スコアが0.59であることは、モデルがデータの変動性の約59%を説明できていることを意味します。この数値を見ると、モデルがデータをある程度は捉えているとは言えますが、そんなに良いわけでもありません。逆にいうと、残りの41%はモデルが説明できないということです。

この一因として、線形モデル自体がこのデータセットには向いていない可能性があります。また、利用可能な特徴量だけではデータセット内のすべてのパターンや関係を完全に捉えきれていないのかもしれません。データセットがこれしかないのであれば、モデルの改善の余地があることを意味します。

なお、R2スコアの定義は以下になります。

$$
R^2 = 1 - \dfrac{\sum\limits_{i=1}^{n} \bigl(y_i - \hat{y}_i\bigr)^2}{\sum\limits_{i=1}^{n} \bigl(y_i - \bar{y}\bigr)^2} 
$$

ここで、

• $${ y_i }$$ はデータセットの各サンプルの実際の目的変数の値です。

• $${ \hat{y}_i }$$ はモデルによる予測値です。

• $${\bar{y} }$$ は実際の目的変数の値の平均です。

• $${n }$$ はサンプルの総数です。

分子 $${ \sum\limits_{i=1}^{n} (y_i - \hat{y}_i)^2 }$$ は、残差平方和（RSS: Residual Sum of Squares）と呼ばれ、モデルの予測値と実際の値の差（誤差）の二乗和です。これは、モデルがどれだけデータポイントから外れているかを示します。

分母 $${ \sum\limits_{i=1}^{n} (y_i - \bar{y})^2 }$$ は、全変動（TSS: Total Sum of Squares）と呼ばれ、実際の値がその平均からどれだけ離れているかの二乗和です。これは、データの変動性の全体的な尺度です。

簡単に言うと、R2スコアは、モデルがデータの変動をどれだけ説明できているかを示します。

よって、R2スコアは、平均二乗誤差と共に回帰モデルの性能を評価するために使用される重要な指標です。平均二乗誤差は予測誤差の大きさの評価であり、R2スコアはモデルの説明力の指標となります。両方とも他のモデルとの比較などにも使うことができます。

次回予告

次回は、分類を行うためのロジスティック回帰について解説します。

お楽しみに！