Pythonと機械学習で東京23区の地価（住宅地）を予測する①：公的データとCatBoostで挑む不動産評価のDX

# 1. はじめに
Pythonと機械学習を習得しました。
本稿では、地価公示・地価調査といった公的データを用い、今回は東京23区住宅地の地価をテーマに分析を実施。
その結果を解説し、不動産鑑定評価を補完する簡易検証ツールの可能性を検討します。

# 2. 解決したい課題
目的は、土地価格評価の客観性と効率性向上を目指します。
現状の不動産鑑定は専門性に依存しますが、データ活用による改善余地が大きいと考えます。
これらの課題に対し、Pythonと機械学習で地価データから価格形成要因をモデル化し、データ分析を通じ、鑑定評価の初期段階における客観的根拠を提供することを目指します。

# 3. 実行環境
パソコン: Windows
開発環境: Google Colaboratory
言語: Python
ライブラリ: Pandas, scikit-learn (sklearn.model_selection, sklearn.metrics), Matplotlib, Seaborn, CatBoost, NumPy, Folium

# 4. 利用データ
本モデル構築には、以下の公開データソースを利用しました。

**・地価公示・地価調査データ:** 
　国土交通省の国土数値情報ダウンロードサイトより取得した、公的な土地の価格情報です。
**・犯罪認知件数:** 
　警視庁サイトより取得した、地域別の犯罪発生状況データです。
**・標高データ:** 
　国土地理院より取得した、土地の物理的な高さに関するデータです。 
**・東京駅からの直線距離:**
　Google Mapsを利用して計算し取得した、対象地からの地理的距離データです。
**・固定資産税路線価・相続税路線価・借地権割合:**
　全国地価マップから取得した、土地の公的な評価額と権利割合に関するデータです。
**・新聞地域:**
　東京新聞の配達地域を参考に分類した地域区分データです。

## データの冒頭部分（df.head()）
```
| chimei   |     ido |   keido |   chiseki | tochi_keijo   | douro_syurui   | douro_houi   |   douro_fukuin | sokudo   |   eki_kyori | youto_chiki   | bouka_chiki   |   kenpei_yoseki |   kakaku_r07 |   kakaku_r06 |   kakaku_r05 |   hyoukou |   hanzai_R6 |   TT_eki_kyori | shinbun_chiki   |   rosenka |   kotei |   syakutiken |
|:---------|--------:|--------:|----------:|:--------------|:---------------|:-------------|---------------:|:---------|------------:|:--------------|:--------------|----------------:|-------------:|-------------:|-------------:|----------:|------------:|---------------:|:----------------|----------:|--------:|-------------:|
| 千代田   | 35.6901 | 139.745 |       969 | 長方形        | 区道           | 南           |           11   | 一方路   |         500 | 2住居         | 防火          |           24000 |      3960000 |      3600000 |      3340000 |      30.1 |        2532 |           2248 | 都心版          |   2880000 | 2330000 |           70 |
| 千代田   | 35.6812 | 139.738 |       535 | 正方形        | 区道           | 北           |            5.4 | 一方路   |         300 | 2住居         | 防火          |           24000 |      2530000 |      2260000 |      2110000 |      29.7 |        2532 |           2680 | 都心版          |   1810000 | 1470000 |           70 |
| 千代田   | 35.6881 | 139.733 |      2266 | 正方形        | 区道           | 北西         |           15   | 一方路   |         330 | 2住居         | 防火          |           24000 |      4830000 |      4390000 |      4280000 |      11.9 |        2532 |           3194 | 都心版          |   3510000 | 2840000 |           70 |
| 千代田   | 35.6986 | 139.746 |       142 | 長方形        | 区道           | 南東         |            8   | 角地     |         360 | 1住居         | 防火          |           24000 |      1970000 |      1740000 |      1620000 |      24.7 |        2532 |           2688 | 都心版          |   1350000 | 1070000 |           70 |
| 千代田   | 35.6961 | 139.746 |      1511 | 台形          | 区道           | 北           |           12.5 | 角地     |         450 | 1住居         | 防火          |           24000 |      3680000 |      3340000 |      3110000 |      23.4 |        2532 |           2513 | 都心版          |   2470000 | 2000000 |           70 |
```

## データ型の確認と欠損値の有無（df.info()）
```
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1222 entries, 0 to 1221
Data columns (total 23 columns):
 #   Column         Non-Null Count  Dtype  
---  ------         --------------  -----  
 0   chimei         1222 non-null   object 
 1   ido            1222 non-null   float64
 2   keido          1222 non-null   float64
 3   chiseki        1222 non-null   int64  
 4   tochi_keijo    1222 non-null   object 
 5   douro_syurui   1222 non-null   object 
 6   douro_houi     1222 non-null   object 
 7   douro_fukuin   1222 non-null   float64
 8   sokudo         1222 non-null   object 
 9   eki_kyori      1222 non-null   int64  
 10  youto_chiki    1222 non-null   object 
 11  bouka_chiki    1222 non-null   object 
 12  kenpei_yoseki  1222 non-null   int64  
 13  kakaku_r07     1222 non-null   int64  
 14  kakaku_r06     1222 non-null   int64  
 15  kakaku_r05     1222 non-null   int64  
 16  hyoukou        1222 non-null   float64
 17  hanzai_R6      1222 non-null   int64  
 18  TT_eki_kyori   1222 non-null   int64  
 19  shinbun_chiki  1222 non-null   object 
 20  rosenka        1222 non-null   int64  
 21  kotei          1222 non-null   int64  
 22  syakutiken     1222 non-null   int64  
dtypes: float64(4), int64(11), object(8)
memory usage: 219.7+ KB
```

## 統計情報の概要（df.describe()）
```
|       |         ido |        keido |   chiseki |   douro_fukuin |   eki_kyori |   kenpei_yoseki |   kakaku_r07 |    kakaku_r06 |    kakaku_r05 |   hyoukou |   hanzai_R6 |   TT_eki_kyori |       rosenka |         kotei |   syakutiken |
|:------|------------:|-------------:|----------:|---------------:|------------:|----------------:|-------------:|--------------:|--------------:|----------:|------------:|---------------:|--------------:|--------------:|-------------:|
| count | 1222        | 1222         |  1222     |     1222       |    1222     |         1222    |   1222       |   1222        |   1222        | 1222      |     1222    |        1222    |   1222        |   1222        |   1222       |
| mean  |   35.6973   |  139.72      |   285.1   |        6.31162 |     720.548 |        12457.1  | 749891       | 688783        | 654135        |   21.7318 |     3463.31 |       10515.5  | 549542        | 443056        |     64.5581  |
| std   |    0.059523 |    0.0852161 |   618.067 |        4.06468 |     435.469 |         6578.51 | 583989       | 521756        | 490404        |   17.6415 |     1079.46 |        3727.74 | 400416        | 326000        |      4.99888 |
| min   |   35.5435   |  139.569     |    47     |        2       |       0     |         3200    | 174000       | 169000        | 165000        |   -3      |     1194    |        1752    | 105000        |  18600        |     60       |
| 25%   |   35.654    |  139.654     |   110     |        4.4     |     430     |         9000    | 438000       | 416000        | 398000        |    2.2    |     2479    |        7801    | 340000        | 272250        |     60       |
| 50%   |   35.7035   |  139.708     |   148     |        5.5     |     600     |        12000    | 586000       | 546000        | 523000        |   24.5    |     3662    |       10804.5  | 440000        | 356000        |     60       |
| 75%   |   35.745    |  139.782     |   199     |        6.075   |     900     |        18000    | 819000       | 752750        | 717750        |   37.4    |     4370    |       13185    | 610000        | 487000        |     70       |
| max   |   35.8141   |  139.915     | 11464     |       40       |    3700     |        48000    |      5.9e+06 |      5.35e+06 |      5.12e+06 |   55      |     6025    |       19822    |      4.28e+06 |      3.58e+06 |     80       |
```
# 5. 分析の流れ
機械学習を活用した分析の流れは、以下のとおりです。
・データ確認
・特徴量エンジニアリングと前処理
・データ分割と可視化
・モデル構築と評価（初回）
・さらなる前処理と再評価
・データフィルタリングとモデル再構築
・予測結果の可視化

# 6. データ確認
まずは必要なライブラリのインポートと、上記のデータのCSVファイルを読込みます。
```python 
# 0. 初めに:
# 必要なライブラリをインストール
pip install catboost
pip install japanize-matplotlib
``` 
```python 
# 1. ライブラリのインポート:
# データ処理、可視化、機械学習モデル構築に必要なライブラリを読み込む
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import japanize_matplotlib
import seaborn as sns

from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score # 同じモジュールからはまとめてインポート
from catboost import CatBoostRegressor

import folium
```
```python 
# 2. CSVファイルの読み込み:
# 分析に使用するCSVファイルを指定したエンコーディングで読み込む。
# エンコーディングエラー発生時は、原因を特定しやすくするためメッセージを表示し終了。
try:
    data = pd.read_csv('R07tikako_tokyo.csv', encoding='Shift-JIS')
    print("CSVファイルを 'Shift-JIS' で読み込みました。")
except UnicodeDecodeError:
    # Shift-JISでの読み込みに失敗した場合のみエラーメッセージを表示
    print("エラー: CSVファイルを 'Shift-JIS' で読み込めませんでした。エンコーディングを確認してください。")
    exit() # プログラムを終了
``` 


# 7. 特徴量エンジニアリングと前処理
地価予測モデルの精度を最大限に引き出すため、データを機械学習モデルが理解しやすい形に加工・変換する「特徴量エンジニアリング」と、欠損値などの問題に対処する「前処理」は極めて重要です。
このセクションでは、具体的な処理内容と、いくつかの試行錯誤の過程を解説します。

### 欠損値の処理
データには、情報が欠けている欠損値が含まれます。これらを適切に処理しないと、モデルの学習が不安定になったり、予測精度が低下したりする原因となります。

### 数値型特徴量の欠損値処理
hanzai_R6（令和6年の犯罪発生件数）、TT_eki_kyori（東京駅からの直線距離）、rosenka（相続税路線価）、kotei（固定資産税路線価）、syakutiken（借地権割合）といった数値型特徴量の欠損値は、各列の平均値で補完しました。
```python 
# 3. 数値特徴量の欠損値処理:
# 数値型の特徴量について、欠損値（NaN）を各列の平均値で補完する。
# これにより、欠損値によるモデル学習の阻害を防ぐ。
numerical_new_features = ['hanzai_R6', 'TT_eki_kyori', 'rosenka', 'kotei', 'syakutiken']
for col in numerical_new_features:
    if col in data.columns and data[col].isnull().any():
        data[col].fillna(data[col].mean(), inplace=True)

# 処理後の欠損値の確認（例として、処理対象の特徴量の欠損値が0になったことを示す）
print("\n数値特徴量処理後の欠損値数:")
print(data[numerical_new_features].isnull().sum())
print("-" * 50)
``` 
```text
数値特徴量処理後の欠損値数:
hanzai_R6       0
TT_eki_kyori    0
rosenka         0
kotei           0
syakutiken      0
dtype: int64
``` 
### カテゴリカル特徴量の欠損値処理
shinbun_chiki（東京新聞地域）、douro_houi（道路方位）、douro_syurui（道路種類）、tochi_keijo（土地形状）、bouka_chiki（防火地域）、youto_chiki（用途地域）といったカテゴリカル型の特徴量に対しては、その列で最も頻繁に出現する値である最頻値で欠損値を補完しました。
```python 
# 4. カテゴリカル特徴量の欠損値処理:
# カテゴリカル型の特徴量について、欠損値（NaN）を各列の最頻値で補完する。
# 数値型とは異なる補完方法を用いることで、データの特性に合わせた適切な前処理を行う。
print("4. 特徴量の欠損値処理と順序エンコーディング:")
categorical_features_for_mode = ['shinbun_chiki', 'douro_houi', 'douro_syurui', 'tochi_keijo', 'bouka_chiki', 'youto_chiki']

for col in data.columns: # 全カラムをループ
    if col in categorical_features_for_mode and data[col].isnull().any(): # 指定したカラムで欠損値があれば
        data[col].fillna(data[col].mode()[0], inplace=True) # 最頻値で補完

# 処理後の欠損値の確認
print("\nカテゴリカル特徴量処理後の欠損値数:")
print(data[categorical_features_for_mode].isnull().sum())
print("-" * 50)
``` 
```text
カテゴリカル特徴量処理後の欠損値数:
shinbun_chiki    0
douro_houi       0
douro_syurui     0
tochi_keijo      0
bouka_chiki      0
youto_chiki      0
dtype: int64
``` 

## 7.1. 各カテゴリカル変数の処理について（重要）
データに含まれるカテゴリカル変数（分類やグループを示す非数値データ）は、機械学習モデルにそのままでは利用できません。
そのため、これらを数値に変換する「エンコーディング」という処理が必要です。
このステップはモデルの性能に大きく影響するため、慎重な検討と試行錯誤を重ねました。

### 地名を地価ランキングで数値化
chimei（地名）は、東京23区という広い範囲を扱う上で非常に重要な情報です。
しかし、地名の種類が多岐にわたるため、そのままカテゴリカル変数として扱うとモデルが複雑になりすぎたり、過学習を招く可能性があります。

そこで、各地名ごとの平均地価を算出し、その平均地価が高い順にランキングを付与して数値化することにしました。
このchimei_rank（地名ランキング）という新しい特徴量は、地価と非常に高い相関を示すことが期待され、実際にモデルにとって非常に有効な情報となりました。
```python 
# 5. 地名のランキング化:
# 各地名の平均地価（kakaku_r07）に基づき、地名を地価が高い順にランキング化して数値特徴量を作成する。
# これは地価予測において非常に重要な特徴量となる。
chimei_avg_price = data.groupby('chimei')['kakaku_r07'].mean().sort_values(ascending=False)
chimei_ranking = {chimei: rank + 1 for rank, chimei in enumerate(chimei_avg_price.index)}
data['chimei_rank'] = data['chimei'].map(chimei_ranking)

print("地名ランキングの確認（上位5件）:")
print(data[['chimei', 'chimei_rank', 'kakaku_r07']].head()) # chimerank と kakaku_r07 を比較
print("-" * 50)
``` 
```text
地名ランキングの確認（上位5件）:
  chimei  chimei_rank  kakaku_r07
0    千代田            1     3960000
1    千代田            1     2530000
2    千代田            1     4830000
3    千代田            1     1970000
4    千代田            1     3680000
``` 
:::note info
このように、地名が地価に応じて数値化され、モデルに地価の地域特性を伝える重要な手がかりとなりました。
:::
### その他のカテゴリカル変数へのアプローチ（CatBoostの活用）
douro_houi（道路方位）、douro_syurui（道路種類）、tochi_keijo（土地形状）、bouka_chiki（防火地域）、youto_chiki（用途地域）といったカテゴリカル変数については、地価との関係性を考慮した数値変換を試みることが一般的です。
私も当初、これらの一部に順序エンコーディング（特定の順序を付与した数値への変換）を適用するアプローチを検討しました。

しかし、今回採用したCatBoostモデルは、カテゴリカル変数を直接扱う能力に非常に優れているという強力な特徴を持っています。
CatBoostは、内部的にこれらのカテゴリ間の複雑な関係性や、目的変数との関連性を自動で学習するため、無理に手動で順序を付与して数値化する必要がありません。

複数の方法を試した結果、これらのカテゴリカル変数をそのままの形でCatBoostに渡した方が、最終的なモデルの予測精度（特に元のスケールでのRMSE）が最も高くなることが判明しました。
これは、私が直感的に設定した順序よりも、CatBoostの内部アルゴリズムがデータからより最適な関連性を見つけ出したためだと考えられます。

この知見に基づき、douro_houi（道路方位）、douro_syurui（道路種類）、tochi_keijo（土地形状）、bouka_chiki（防火地域）、youto_chiki（用途地域）については、欠損値処理（最頻値補完）のみを行い、数値変換は行わずにCatBoostに任せる方針を採用しました。

一方で、shinbun_chiki（新聞地域）については、例外的に数値変換（順序エンコーディング）を適用しました。これは、「都心版」「山手版」「したまち版」といったカテゴリ間に、地価に対する明確な影響度の順序が存在すると考えたためです。
この順序付けにより、モデルが地域の特性をより効率的に学習できると判断しました。
```python 
# 6. 新聞地域の順序エンコーディング:
# 新聞地域カテゴリを地価への影響度を考慮して順序付けし、数値に変換する (都心版:0, 山手版:1, したまち版:2)
# 変換後、元の列は削除する。
shinbun_mapping = {'都心版': 0, '山手版': 1, 'したまち版': 2}
if 'shinbun_chiki' in data.columns:
    data['shinbun_chiki_rank'] = data['shinbun_chiki'].map(shinbun_mapping)
    data.drop('shinbun_chiki', axis=1, inplace=True) # 元の列を削除
``` 
以下に、数値変換を行わないと判断したカテゴリカル変数に対する処理を再掲します。
これらは、欠損値の最頻値補完のみを行い、モデルに直接渡されます。
```python 
# 7. 道路方位の欠損値処理:
# カテゴリカル変数である「道路方位」の欠損値を最頻値で補完する。
# CatBoostはカテゴリカル変数を直接効率的に扱えるため、ここでは数値への変換は行わない。
if 'douro_houi' in data.columns and data['douro_houi'].isnull().any():
    most_frequent_houi = data['douro_houi'].mode()[0]
    data['douro_houi'].fillna(most_frequent_houi, inplace=True)

# 8. 道路種類の欠損値処理:
# カテゴリカル変数である「道路種類」の欠損値を最頻値で補完する。
# CatBoostはカテゴリカル変数を直接効率的に扱えるため、ここでは順序エンコーディングは行わない。
if 'douro_syurui' in data.columns and data['douro_syurui'].isnull().any():
    most_frequent_syurui = data['douro_syurui'].mode()[0]
    data['douro_syurui'].fillna(most_frequent_syurui, inplace=True)

# 9. 土地形状の欠損値処理:
# カテゴリカル変数である「土地形状」の欠損値を最頻値で補完する。
# CatBoostはカテゴリカル変数を直接効率的に扱えるため、ここでは順序エンコーディングは行わない。
if 'tochi_keijo' in data.columns and data['tochi_keijo'].isnull().any():
    most_frequent_keijo = data['tochi_keijo'].mode()[0]
    data['tochi_keijo'].fillna(most_frequent_keijo, inplace=True)

# 10. 防火地域の欠損値処理:
# カテゴリカル変数である「防火地域」の欠損値を最頻値で補完する。
# CatBoostはカテゴリカル変数を直接効率的に扱えるため、ここでは順序エンコーディングは行わない。
if 'bouka_chiki' in data.columns and data['bouka_chiki'].isnull().any():
    most_frequent_bouka = data['bouka_chiki'].mode()[0]
    data['bouka_chiki'].fillna(most_frequent_bouka, inplace=True)

# 11. 用途地域の欠損値処理:
# カテゴリカル変数である「用途地域」の欠損値を最頻値で補完する。
# CatBoostはカテゴリカル変数を直接効率的に扱えるため、ここでは数値への変換は行わない。
if 'youto_chiki' in data.columns and data['youto_chiki'].isnull().any():
    most_frequent_youto = data['youto_chiki'].mode()[0]
    data['youto_chiki'].fillna(most_frequent_youto, inplace=True)
``` 
:::note info
このアプローチは、モデルの簡潔さを保ちつつ、高い予測精度を実現する上で非常に有効でした。
:::
# 8. データ分割と可視化
前処理が完了したデータセットを、機械学習モデルの学習と評価利用できる形に整えます。
ここでは、まずデータ全体の最終的な状態を確認し、その後、データを学習用とテスト用に分割し、モデル構築前にデータの特性を視覚的に把握するための可視化を行います。

### データの最終確認
前処理によってデータが適切に加工されているか、最初の数行を確認します。
```python 
# 12. データフレームの先頭行を確認:
# データが正しく読み込まれ、前処理が行われているかをsanity checkとして確認する。
print("12. データフレームの先頭行を確認:")
print(data.head())
print("-" * 50)
``` 
```text
12. データフレームの先頭行を確認:
  chimei        ido       keido  chiseki tochi_keijo douro_syurui douro_houi  \
0    千代田  35.690140  139.744815      969         長方形           区道          南   
1    千代田  35.681200  139.737520      535         正方形           区道          北   
2    千代田  35.688141  139.732873     2266         正方形           区道         北西   
3    千代田  35.698633  139.746461      142         長方形           区道         南東   
4    千代田  35.696080  139.746152     1511          台形           区道          北   

   douro_fukuin sokudo  eki_kyori  ... kakaku_r06 kakaku_r05  hyoukou  \
0          11.0    一方路        500  ...    3600000    3340000     30.1   
1           5.4    一方路        300  ...    2260000    2110000     29.7   
2          15.0    一方路        330  ...    4390000    4280000     11.9   
3           8.0     角地        360  ...    1740000    1620000     24.7   
4          12.5     角地        450  ...    3340000    3110000     23.4   

   hanzai_R6  TT_eki_kyori  rosenka    kotei  syakutiken  shinbun_chiki_rank  \
0       2532          2248  2880000  2330000          70                   0   
1       2532          2680  1810000  1470000          70                   0   
2       2532          3194  3510000  2840000          70                   0   
3       2532          2688  1350000  1070000          70                   0   
4       2532          2513  2470000  2000000          70                   0   

   chimei_rank  
0            1  
1            1  
2            1  
3            1  
4            1  

[5 rows x 24 columns]
``` 

### データの分割
モデルが未知のデータに対しても適切に予測できるかを評価するため、全データを学習用 (Train) とテスト用 (Test) に分割します。
今回は、全データの80%を学習用、20%をテスト用とし、stratify=data['chimei']と指定することで、地名が学習用・テスト用データセットに均等に分散されるようにしました。
これにより、特定の地域に偏ったデータ分割を防ぎ、モデルの汎化性能をより適切に評価できます。

また、モデルの目的変数であるkakaku_r07（地価）は、その分布が偏っている（高価格帯のデータが少ないなど）ため、そのまま学習に用いるとモデルの学習が不安定になることがあります。
これを防ぐため、地価には対数変換（np.log1p）を適用しました。
これにより、データの分布が正規分布に近づき、モデルの学習が安定し、予測精度を向上することが期待されます。

:::note warn
np.log1pは、x に 1を足してから自然対数を取るNumPy関数です。
xが非常に小さい値の場合、通常の対数関数より高精度に計算できるのが特徴。
データ分析で、値が0を含む、または0に近いデータの分布を整える（正規分布に近づける）際によく使われます。
:::
```python 
# 13. データを学習用とテスト用に分割:
# モデル学習のためのデータセットを、学習用(train)とテスト用(test)に8:2の比率で分割する。
# 'chimei'（地名）で層化抽出することで、各地域のデータが両セットに均等に配分されるようにする。
train_df, test_df = train_test_split(data, test_size=0.2, stratify=data['chimei'], random_state=42)
print("13. trainとtestへの分割:")
print(f"Trainデータの形状: {train_df.shape}")
print(f"Testデータの形状: {test_df.shape}")
print("-" * 50)

# 14. 特徴量と目的変数の分離、および目的変数の対数変換:
# 学習用とテスト用のデータセットから、目的変数('kakaku_r07')を分離し、
# 残りを特徴量(X)とする。目的変数には、分布を正規化しモデルの性能向上を目的として対数変換(log1p)を適用する。
X_train = train_df.drop('kakaku_r07', axis=1)
y_train = np.log1p(train_df['kakaku_r07']) # 目的変数を対数変換

X_test = test_df.drop('kakaku_r07', axis=1)
y_test = np.log1p(test_df['kakaku_r07']) # 目的変数を対数変換
print("14. 学習データの確認:")
print("14. 学習データの確認:")
print("X_trainの形状:", X_train.shape)
print("y_trainの形状:", y_train.shape)
print("X_trainの記述統計:")
print(X_train.describe())
print("-" * 50)
``` 
```text
13. trainとtestへの分割:
Trainデータの形状: (977, 24)
Testデータの形状: (245, 24)
--------------------------------------------------
14. 学習データの確認:
X_trainの形状: (977, 23)
y_trainの形状: (977,)
X_trainの記述統計:
              ido       keido       chiseki  douro_fukuin    eki_kyori  \
count  977.000000  977.000000    977.000000    977.000000   977.000000   
mean    35.696948  139.719945    284.828045      6.357625   725.598772   
std      0.059754    0.084906    642.518169      4.155892   433.058455   
min     35.543465  139.569372     47.000000      2.000000     0.000000   
25%     35.653426  139.655094    110.000000      4.400000   430.000000   
50%     35.703807  139.707821    149.000000      5.500000   620.000000   
75%     35.744704  139.781418    198.000000      6.000000   900.000000   
max     35.814096  139.914638  11464.000000     40.000000  3700.000000   

       kenpei_yoseki    kakaku_r06    kakaku_r05     hyoukou    hanzai_R6  \
count     977.000000  9.770000e+02  9.770000e+02  977.000000   977.000000   
mean    12485.977482  6.880890e+05  6.533664e+05   21.798158  3459.495394   
std      6569.818744  5.144199e+05  4.817951e+05   17.570655  1080.116346   
min      3200.000000  1.690000e+05  1.650000e+05   -3.000000  1194.000000   
25%      9000.000000  4.160000e+05  3.990000e+05    2.200000  2479.000000   
50%     12000.000000  5.470000e+05  5.230000e+05   24.700000  3662.000000   
75%     18000.000000  7.530000e+05  7.180000e+05   37.500000  4370.000000   
max     48000.000000  5.350000e+06  5.120000e+06   55.000000  6025.000000   

       TT_eki_kyori       rosenka         kotei  syakutiken  \
count    977.000000  9.770000e+02  9.770000e+02  977.000000   
mean   10504.253838  5.484186e+05  4.419584e+05   64.626407   
std     3717.803053  3.918397e+05  3.184403e+05    4.988577   
min     1752.000000  1.350000e+05  4.840000e+04   60.000000   
25%     7825.000000  3.400000e+05  2.720000e+05   60.000000   
50%    10796.000000  4.400000e+05  3.570000e+05   60.000000   
75%    13172.000000  6.100000e+05  4.900000e+05   70.000000   
max    19822.000000  4.280000e+06  3.580000e+06   70.000000   

       shinbun_chiki_rank  chimei_rank  
count          977.000000   977.000000  
mean             1.123849    14.420676  
std              0.689510     6.177871  
min              0.000000     1.000000  
25%              1.000000    11.000000  
50%              1.000000    14.000000  
75%              2.000000    20.000000  
max              2.000000    23.000000  
``` 

### データの可視化
モデル構築前に、データの特性や特徴量間の関係性を理解することは、予測モデルの解釈や改善に役立ちます。

### 目的変数（地価）の分布
対数変換を施した目的変数y_trainの分布をヒストグラムで確認します。
これにより、地価の分布がより正規分布に近づいていることを視覚的に確認できます。
```python 
# 15. 目的変数の分布を可視化:
# 対数変換された目的変数（地価）の分布をヒストグラムとKDE（カーネル密度推定）で確認する。
# これにより、データの分布が正規分布に近いか、外れ値がないかなどを把握できる。
# 目的変数 (kakaku_r07) の分布
plt.figure(figsize=(8, 6))
sns.histplot(y_train, kde=True)
plt.title('目的変数 (kakaku_r07) の分布')
plt.xlabel('価格 (R07)')
plt.ylabel('頻度')
plt.show()
``` 
![#15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4091959/6cad2778-7bc5-4704-9480-a0b412d2de76.png)


### 数値変数の散布図行列

主要な数値特徴量と目的変数との関係性を一度に確認するため、散布図行列を作成します。
ここでは、地価形成に特に大きな影響を与えると考えられる主要な数値特徴量を選び、地価（kakaku_r07）との関係性を散布図行列で視覚的に確認します。
これにより、特徴量間の相関の有無や、地価との線形関係の傾向などを視覚的に把握できます。
```python 
#16. 主要な数値変数の関係性可視化:
# 地価予測において最も影響力の高い主要な数値特徴量を選び、地価（kakaku_r07）との関係性を散布図行列で可視化します。
# これにより、地価を決定づける核となる要因間の相関と分布を直感的に把握できます。

# 厳選された主要な数値特徴量を選択
selected_numerical_features = [
    'kakaku_r06',      # 昨年度価格
    'rosenka',         # 相続税路線価
    'kotei',           # 固定資産税路線価
    'eki_kyori',       # 最寄り駅までの距離
    'TT_eki_kyori'     # 東京駅からの直線距離
]

# 散布図行列を作成
# train_dfから選択した特徴量と目的変数'kakaku_r07'を抽出
sns.pairplot(train_df[selected_numerical_features + ['kakaku_r07']])
plt.suptitle('主要な数値変数の散布図行列', y=1.02)
plt.show()
print("-" * 50)
``` 
![#16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4091959/c856ae5f-ac91-4bdc-818f-3c3ca0c442f7.png)

## 散布図行列から読み取れる相関関係

#### 価格・評価指標の圧倒的な影響:
kakaku_r06（昨年度価格）、rosenka（相続税路線価）、kotei（固定資産税路線価）は、kakaku_r07（当年度価格）と極めて強い正の相関を示し、地価の主要な決定要因であることが視覚的に確認できます。
ヒストグラムからも、データが低い価格帯に集中している傾向が見られます。

#### 立地条件の重要性:
eki_kyori（最寄り駅までの距離）やTT_eki_kyori（東京駅からの直線距離）とkakaku_r07の間には負の相関が見られ、距離が遠くなるほど地価が下がる傾向を確認できます。
これらの距離は、地価の多様性に大きく寄与しています。

### 欠損値の最終確認と相関係数
これまでの前処理で欠損値は対処されているはずですが、モデル学習の直前には、データに欠損値が残っていないかを再度確認することが重要です。
また、散布図行列で視覚的に捉えた相関関係をより厳密に分析するため、数値特徴量間の相関関係をヒートマップで可視化します。
これにより、多重共線性（特徴量同士が高い相関を持つこと）の可能性や、地価との関係性を数値として明確に把握できます。
```python 
# 17. 欠損値の最終確認:
# 前処理後の学習用データセットに欠損値が残っていないかを最終確認する。
# これにより、モデル学習前のデータ品質を保証する。
print("17. 欠損値確認:")
print(train_df.isnull().sum())
print("-" * 50)
``` 
```text
17. 欠損値確認:
chimei                0
ido                   0
keido                 0
chiseki               0
tochi_keijo           0
douro_syurui          0
douro_houi            0
douro_fukuin          0
sokudo                0
eki_kyori             0
youto_chiki           0
bouka_chiki           0
kenpei_yoseki         0
kakaku_r07            0
kakaku_r06            0
kakaku_r05            0
hyoukou               0
hanzai_R6             0
TT_eki_kyori          0
rosenka               0
kotei                 0
syakutiken            0
shinbun_chiki_rank    0
chimei_rank           0
dtype: int64
``` 
```python 
# 18. 数値特徴量の相関係数を確認:
# 各数値特徴量間の線形な関係性（相関）をヒートマップで可視化する。
# これにより、多重共線性の有無や、目的変数との関連性の強さを把握できる。
print("18. 相関係数:")
correlation_matrix = train_df.corr(numeric_only=True)
plt.figure(figsize=(16, 14))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt=".2f")
plt.title('相関係数行列')
plt.show()
print("-" * 50)
``` 
![#18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4091959/08603983-75be-406e-babe-ebfc8f809640.png)

## ヒートマップから読み取れる相関関係
数値特徴量の相関係数ヒートマップからは、各変数間の関係性が視覚的に把握できます。
ここでは、特に気になる相関関係について解説します。
#### chiseki (地積):
kakaku_r07との相関は0.28と弱いながら正の相関が見られます。
都心に近くなるほど地価は高い傾向にあります。そして、都心に近いほどマンション等の中高層集合住宅の割合が高まり、地積の大きい土地への需要が増すためと考えられます。
#### douro_fukuin (道路幅員): 
kakaku_r07との相関は0.30と弱いながら正の相関が見られます。
地積の場合と同様に、都心に近いほどマンション等の中高層集合住宅の割合が高く、道路幅員の広い需要が増すためと考えられます。
#### eki_kyori (最寄り駅までの距離): 
kakaku_r07との相関は-0.33と負の相関が見られます。
どの地域も駅から遠くなるほど地価が安くなるという予想通りの結果になったと考えられます。
#### kenpei_yoseki (建ぺい率×容積率): 
kakaku_r07との相関は0.41と中程度の正の相関が見られます。
都心に近いほど中高層住宅を建築するメリットが大きく、価値が高いということを示していると考えられます。
#### kakaku_r06 (前年度価格)、kakaku_r05 (2年前の価格):
それぞれ極めて高い正の相関を示しており、地価が過去の価格推移と強く連動していることが明確にわかります。これは直感とも合致する結果です。
#### rosenka (相続税路線価)、kotei (固定資産税路線価):
それぞれ高い正の相関を示しており、公的な評価額が地価に大きく影響していることを裏付けていることがわかります。
#### TT_eki_kyori (東京駅からの直線距離):
kakaku_r07との相関は-0.56と中程度の負の相関が見られます。
都心からの距離が遠くなるほど地価が安くなる傾向を示しています。
#### shinbun_chiki_rank (新聞地域ランキング):
kakaku_r07との相関は-0.42と中程度の負の相関が見られます。
東京23区は、「したまち」「山手」「都心」という呼称に象徴されるように、地価に及ぼす影響が他都市と比較して大きいのではないかと推測されます。
#### chimei_rank (地名ランキング):
kakaku_r07との相関は-0.73と強い負の相関が見られます。
これは、ランキング値が低い（例: 1位、2位）ほど地価が高い地域であることを示唆していると考えられます。

```python 
# 19. X_train, X_testの欠損値処理:
# train_test_splitで分割された特徴量データ（X_train, X_test）について、
# 数値型の欠損値をそれぞれの列の平均値で最終的に補完する。
# これにより、モデル学習前にデータに欠損がない状態を確実にする。
print("19. 前処理:")
# 欠損値処理 (平均値で補完)
for col in X_train.select_dtypes(include=np.number).columns:
    X_train.fillna({col: X_train[col].mean()}, inplace=True)
for col in X_test.select_dtypes(include=np.number).columns:
    X_test.fillna({col: X_test[col].mean()}, inplace=True)

# 学習後の欠損値確認
print("学習データの前処理後の欠損値:")
print(X_train.isnull().sum())
print("学習データのカラム:")
print(X_train.columns)
print("-" * 50)
``` 
```text
19. 前処理:
学習データの前処理後の欠損値:
chimei                0
ido                   0
keido                 0
chiseki               0
tochi_keijo           0
douro_syurui          0
douro_houi            0
douro_fukuin          0
sokudo                0
eki_kyori             0
youto_chiki           0
bouka_chiki           0
kenpei_yoseki         0
kakaku_r06            0
kakaku_r05            0
hyoukou               0
hanzai_R6             0
TT_eki_kyori          0
rosenka               0
kotei                 0
syakutiken            0
shinbun_chiki_rank    0
chimei_rank           0
dtype: int64
学習データのカラム:
Index(['chimei', 'ido', 'keido', 'chiseki', 'tochi_keijo', 'douro_syurui',
       'douro_houi', 'douro_fukuin', 'sokudo', 'eki_kyori', 'youto_chiki',
       'bouka_chiki', 'kenpei_yoseki', 'kakaku_r06', 'kakaku_r05', 'hyoukou',
       'hanzai_R6', 'TT_eki_kyori', 'rosenka', 'kotei', 'syakutiken',
       'shinbun_chiki_rank', 'chimei_rank'],
      dtype='object')
``` 
:::note info
これらの可視化と最終確認を通じて、モデル学習に進むためのデータが適切に準備されたことを確認しました。
:::
# 9. モデル構築と評価（初回）
前処理とデータ分割が完了した学習用データを用いて、いよいよ機械学習モデルの構築と評価を行います。
ここでは、CatBoostを回帰モデルとして採用し、その初期性能と、地価予測においてどの特徴量が重要であるかを確認します。

### CatBoostモデルの学習と初期評価
CatBoostは、カテゴリカル変数を直接扱えるという強力な特徴を持つ勾配ブースティング決定木モデルです。
前処理のセクションで述べたように、私はこのCatBoostの特性を最大限に活かすため、多くのカテゴリカル変数を数値変換せずにそのままモデルに渡しました。

モデルの学習後、テストデータに対する予測を行い、RMSE（二乗平均平方根誤差）を計算してモデルの予測精度を評価します。

:::note warn
RMSEは予測値と実際の値との差の平均的な大きさを表し、値が小さいほど高精度なモデルであることを意味します。
:::
```python 
# 20. CatBoostモデルの学習と評価:
# CatBoostRegressorモデルを初期化し、学習データ（X_train, y_train）で学習させます。
# その際、X_train内のカテゴリカル変数（object型）をCatBoostに正しく認識させます。
# 学習後、テストデータで予測を行い、RMSE（二乗平均平方根誤差）を計算・表示し、モデルの性能を評価します。
print("20. CatBoostによるスコア確認:")
categorical_features_indices = [X_train.columns.get_loc(col) for col in X_train.select_dtypes(include='object').columns]

model = CatBoostRegressor(iterations=100, random_state=42, verbose=0)
model.fit(X_train, y_train, cat_features=categorical_features_indices)

# y_trainとy_testを元のスケールに戻す（計算は一度でOK）
y_train_original = np.expm1(y_train)
y_test_original = np.expm1(y_test)

# 学習データでの予測値を元のスケールに戻す
y_train_pred = model.predict(X_train)
y_train_pred_original = np.expm1(y_train_pred)

# テストデータでの予測値を元のスケールに戻す
y_pred = model.predict(X_test)
y_pred_original = np.expm1(y_pred)

# RMSEの計算と表示
rmse_train_log = np.sqrt(mean_squared_error(y_train, y_train_pred))
rmse_test_log = np.sqrt(mean_squared_error(y_test, y_pred))

rmse_train_original_scale = np.sqrt(mean_squared_error(y_train_original, y_train_pred_original))
rmse_test_original_scale = np.sqrt(mean_squared_error(y_test_original, y_pred_original))


print(f"CatBoost RMSE (学習データ - 対数スケール): {rmse_train_log:.5f}")
print(f"CatBoost RMSE (テストデータ - 対数スケール): {rmse_test_log:.5f}")
print(f"CatBoost RMSE (学習データ - 元のスケール): {int(rmse_train_original_scale):,} 円") # 小数点以下を切り捨て
print(f"CatBoost RMSE (テストデータ - 元のスケール): {int(rmse_test_original_scale):,} 円") # 小数点以下を切り捨て
print("-" * 50)

# RMSEは二乗平均平方根誤差(Root Mean Squared Error)を示し、予測値と実際の値との差の平均的な大きさを表す。
# 値が小さいほどモデルの予測精度が高いことを意味する。
# RMSEの数値は小数点以下を分かりやすくするために、適切な桁数（ここでは5桁）に丸めて表示した
``` 
```text
20. CatBoostによるスコア確認:
CatBoost RMSE (学習データ - 対数スケール): 0.02290
CatBoost RMSE (テストデータ - 対数スケール): 0.03975
CatBoost RMSE (学習データ - 元のスケール): 26,256 円
CatBoost RMSE (テストデータ - 元のスケール): 44,250 円
``` 
:::note info
対数変換後のRMSEは、**学習データで0.02290**、**テストデータで0.03975**と低い値を示しています。
これを元のスケールに戻して実測すると、**学習データでは26,256円**、**テストデータでは44,250円**の予測誤差となることが示されました。
このRMSEは、地価予測において非常に高い精度が期待できることを意味します。

また、学習データとテストデータのRMSE（特に元のスケールでの誤差）に大きな乖離が見られないことから、**モデルが特定のデータに過度に適合する「過学習」が抑制されている**ことも確認できます。
これは、**未知のデータに対しても安定した予測性能が期待できることを意味します。**
:::

### 特徴量の重要度分析
次に、モデルが地価予測を行う上で、どの特徴量を特に重視しているのかを特徴量重要度から分析します。
これは、モデルの予測根拠を理解し、地価形成要因を考察する上で**非常に重要なステップ**です。
```python 
# 20.x. 特徴量の重要度を表示:
# モデルが学習した結果、どの特徴量が地価予測に最も貢献しているかを数値とグラフで可視化する。
# これにより、モデルの予測根拠を理解し、重要な特徴量を特定できる。
print("20.x. 特徴量の重要度:")

# モデルから特徴量の重要度を取得
feature_importances = model.get_feature_importance()

# 特徴量名を取得
feature_names = X_train.columns

# 特徴量の重要度をDataFrameにまとめる
importance_df = pd.DataFrame({
    '特徴量': feature_names,
    '重要度': feature_importances
})

# 重要度で降順にソート
importance_df = importance_df.sort_values(by='重要度', ascending=False)

# 上位の特徴量を表示（例: 上位20個）
print(importance_df.head(20))

# 可視化（オプション）
plt.figure(figsize=(10, 8))
sns.barplot(x='重要度', y='特徴量', data=importance_df.head(20))
plt.title('CatBoost 特徴量重要度 (上位20)') 
plt.xlabel('重要度') 
plt.ylabel('特徴量') 
plt.tight_layout() # レイアウトの調整
plt.show()

print("-" * 50)
``` 
```text
20.x. 特徴量の重要度:
            特徴量        重要度
13          kakaku_r06   27.503953
18             rosenka   26.508388
14          kakaku_r05   22.285387
19               kotei    9.303130
9            eki_kyori    4.063079
1                  ido    2.377384
16           hanzai_R6    1.715996
3              chiseki    1.682931
20          syakutiken    1.456102
22         chimei_rank    0.839577
17        TT_eki_kyori    0.666860
12       kenpei_yoseki    0.440774
11         bouka_chiki    0.295508
2                keido    0.279953
7         douro_fukuin    0.132287
15             hyoukou    0.131238
21  shinbun_chiki_rank    0.095295
8               sokudo    0.084399
10         youto_chiki    0.047766
4          tochi_keijo    0.040159
``` 
![#20x.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4091959/81435f31-21a8-4c03-9af1-2669a3e1044b.png)

## 特徴量重要度からの考察:
この結果から、地価予測において以下の点が明らかになりました。

### 過去の価格データと公的評価指標が圧倒的な重要性:
:::note info
kakaku_r06（昨年度価格）、rosenka（相続税路線価）、kakaku_r05（2年前の価格）、そしてkotei（固定資産税路線価）の**4つの特徴量だけで重要度全体の約85%を上位が占めました**。
これは、**地価が過去の価格推移や公的な評価基準に強く連動しているという直感的な理解を、モデルが明確に裏付けていること**を示しています。
:::

### One-Hot Encodingを用いたモデルとの比較（試行錯誤の経緯）
CatBoostがカテゴリカル変数を直接扱える利点を最大限に活かすため、多くのカテゴリカル変数をそのまま使用する方針を採りましたが、念のため一般的なOne-Hot Encodingを施したデータでもモデルを構築し、比較を行いました。
```python 
# 21. One-Hot Encodingの実施とカラム調整:
# CatBoost以外のモデルを試す場合（またはCatBoostに数値特徴量として扱わせる比較のため）、
# 残りのカテゴリカル変数（object型）をOne-Hot Encodingにより数値データに変換する。
# 学習用とテストデータで生成されるカラムの不一致を修正し、モデル学習に備える。

# カテゴリカル変数のOne-Hot Encoding (今回は 'shinbun_chiki' は順序エンコーディング済み)
categorical_cols_to_encode = X_train.select_dtypes(include='object').columns.tolist()

X_train_encoded = pd.get_dummies(X_train, columns=categorical_cols_to_encode, drop_first=True)
X_test_encoded = pd.get_dummies(X_test, columns=categorical_cols_to_encode, drop_first=True)

# One-Hot Encoding後のカラムのずれを修正
train_cols = set(X_train_encoded.columns)
test_cols = set(X_test_encoded.columns)

missing_in_test = list(train_cols - test_cols)
for col in missing_in_test:
    X_test_encoded[col] = 0
missing_in_train = list(test_cols - train_cols)
for col in missing_in_train:
    X_train_encoded[col] = 0

X_test_encoded = X_test_encoded[X_train_encoded.columns]

print("One-Hot Encoding後の学習データ形状:", X_train_encoded.shape)
print("One-Hot Encoding後のテストデータ形状:", X_test_encoded.shape)
print("-" * 50)
``` 
```text
One-Hot Encoding後の学習データ形状: (977, 68)
One-Hot Encoding後のテストデータ形状: (245, 68)
``` 

### 元のスケールでのRMSE比較（重要な判断基準）
モデルの評価指標として、対数変換後のRMSEだけでなく、元のスケール（円単位）に戻したRMSEも計算しました。
これは、実際に地価予測ツールとして利用する際に、どれくらいの誤差が生じるかを直感的に把握するためです。
```python 
# 22. CatBoostモデルの再学習と評価 (One-Hot Encodingされたデータを使用):
# One-Hot Encodingによって数値変換されたカテゴリカル変数を含むデータセットを用いて、
# 再度CatBoostモデルを学習・評価する。これは、元のカテゴリカル変数を直接扱ったモデルとの
# 性能比較のために行われる。
# 予測結果を元のスケールに戻し、実際の地価予測におけるRMSE（円単位の誤差）も計算・表示する。
model_encoded = CatBoostRegressor(iterations=100, random_state=42, verbose=0)
model_encoded.fit(X_train_encoded, y_train)
y_pred_encoded = model_encoded.predict(X_test_encoded)
rmse_encoded = np.sqrt(mean_squared_error(y_test, y_pred_encoded))
print(f"One-Hot Encoding後のCatBoost RMSE (対数スケール): {rmse_encoded:.5f}") # 小数点以下5桁に調整
print("-" * 50)

# y_test は対数変換されたテストデータの 'kakaku_r07'
# y_pred はモデルによる対数変換された 'kakaku_r07' の予測値

# y_test を元のスケールに戻す（※これは元モデルのy_test_originalと同じなので再計算不要だが、念のため記載）
# y_test_original = np.expm1(y_test) # 既に計算済みと仮定

# y_pred を元のスケールに戻す（これは元のモデルのy_predなので注意。One-Hot Encoding版のy_pred_encodedを使う）
y_pred_encoded_original = np.expm1(y_pred_encoded) # One-Hot Encoding版の予測値を元のスケールに戻す

# 元のスケールでの RMSE を計算
rmse_original_scale_encoded = np.sqrt(mean_squared_error(y_test_original, y_pred_encoded_original)) # y_test_originalは最初のモデル評価時に計算されたものを使用
print(f"One-Hot Encoding後のCatBoost RMSE (元のスケール): {int(rmse_original_scale_encoded):,} 円") # 小数点以下を切り捨て、カンマ区切りで表示
print("-" * 50)
``` 
```text
One-Hot Encoding後のCatBoost RMSE (対数スケール): 0.04377
--------------------------------------------------
One-Hot Encoding後のCatBoost RMSE (元のスケール): 54,417 円
``` 
CatBoostにカテゴリカル変数を直接扱わせたモデルとの性能比較のため、One-Hot Encodingを施したデータで再学習したモデルの評価を行いました。
One-Hot Encoding後のCatBoostモデルのRMSEは、対数スケールで0.04377、元のスケールでは54,417円となりました。

元のスケールおよび対数スケールでのRMSE（円単位の予測誤差）においては、CatBoostにカテゴリカル変数を直接扱わせた方が**より優れた性能を発揮しました**。
最終的な実用性を考慮し、本分析ではカテゴリカル変数を数値変換せず、CatBoostに直接扱わせるアプローチを主モデルとして採用しました。

### その他のモデル評価指標
RMSEに加え、MAE（平均絶対誤差）とR²（決定係数）も計算し、モデルの性能を多角的に評価します。

:::note warn
MAE(平均絶対誤差)は予測値と実際の値の差の平均です。モデルの予測精度を測る指標で、値が小さいほど良いです。
R²スコア(決定係数)はモデルがデータの変動をどれだけ説明できているかを示す指標です。0から1の間で、1に近いほどモデルの当てはまりが良いことを意味します。
:::
```python 
# 23. モデル性能の評価 (MAEとR²スコア):
# モデルの予測精度を多角的に評価するため、平均絶対誤差 (MAE) と決定係数 (R²スコア) を計算し表示する。
# MAEは予測誤差の平均的な絶対値を示し、R²はモデルが目的変数の変動をどれだけ説明できているかを示す。
# MAEは値が小さいほど、R²は1に近いほど（最大1）モデルの性能が良いことを意味する。
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)

print(f"CatBoost MAE (対数スケール): {mae:.5f}") # 小数点以下5桁に調整
print(f"CatBoost R^2: {r2:.5f}") # 小数点以下5桁に調整
print("-" * 50)
``` 
```text
CatBoost MAE (対数スケール): 0.02877
CatBoost R²: 0.99449
``` 
:::note info
**MAEは0.02877（対数スケール）と非常に小さく**、**R²スコアは0.99449と1に極めて近い値**を示しています。R²スコアが0.99を超えることは、**モデルが地価の変動の約99.4%を説明できていることを意味し、これは非常に高い予測精度であることを示しています**。
:::

### 予測結果の可視化
最後に、モデルが予測した地価と実際の地価がどれくらい一致しているかを視覚的に確認します。特に、対数変換を元に戻した元のスケール（円単位）で比較することで、モデルの予測精度をより直感的に評価できます。

散布図では、横軸に実際の地価、縦軸にモデルの予測地価を取ります。理想的な予測であれば、全ての点がy=xの対角線上に並ぶはずです。
```python 
# 24. 予測結果の可視化 (実測値 vs 予測値):
# モデルによる予測値と実際の値の関係性を散布図で可視化する。
# 対数変換を戻した元のスケール（円単位）で比較することで、モデルの予測性能を直感的に評価できる。
# グラフ上にy=xの対角線を引き、点群がこの線に近ければ近いほど、予測精度が高いことを示す。
plt.figure(figsize=(8, 8))
sns.scatterplot(x=y_test_original , y=y_pred_original )
plt.plot([y_test_original .min(), y_test_original .max()], [y_test_original .min(), y_test_original .max()], '--k')
plt.xlabel("実測値") 
plt.ylabel("予測値") 
plt.title("実測値 vs. 予測値") 

# x軸のフォーマットを変更 (日本円のような3桁区切りに)
plt.gca().xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, p: format(int(x), ',')))

# y軸のフォーマットを変更 (日本円のような3桁区切りに)
plt.gca().yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, p: format(int(x), ',')))

plt.show() # グラフを表示

print("-" * 50)
print("300万円以上のデータ数が少ないため、300万円以下でモデルを再度作成する。")
``` 
![#24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4091959/5fc3d40e-a706-4a3b-add9-f402cb4433ea.png)


この散布図から、地価予測モデルが非常に高い精度で地価を捉えられていることが視覚的にも確認できます。一方で、データ数が少ない高価格帯の物件では、予測がややばらつく傾向も見て取れます。
この課題に対応するため、次に300万円以下のデータに限定してモデルを再構築するアプローチを試みました。

:::note warn
300万円以上のデータ数が少ないため、300万円以下でモデルを再度作成する。
:::


# 10. データフィルタリングとモデル再構築
地価データには、高価格帯の物件は数が少なく、データ分布が偏っているという特徴があります。このような少数の高価格帯データが、モデル全体の学習に悪影響を与えたり、予測精度を低下させたりする可能性があります。

そこで、データセットをさらに改善する試みとして、価格が300万円以下の物件に限定してモデルを再構築しました。これにより、比較的データが豊富な価格帯に焦点を当て、その範囲での予測精度を最大化することを目指します。

### データのフィルタリング
kakaku_r07が300万円以下のデータのみを抽出します。
```python 
# 25. 特定の価格帯でデータをフィルタリング:
# データセットの中で、地価（kakaku_r07）が300万円以下の物件に絞り込んで新たなデータセットを作成する。
# これは、データ数が少ない高価格帯の物件を除外することで、モデルの予測精度を向上させることを目的としている。
filtered_data = data[data['kakaku_r07'] <= 3e+06].copy()
print(f"300万円以下のデータでフィルタリング後のデータ数: {len(filtered_data)}")
print("-" * 50)
``` 
```text
300万円以下のデータでフィルタリング後のデータ数: 1206
``` 

### フィルタリング後のデータでのモデル再構築と評価
フィルタリングされたデータセットを再度学習用とテスト用に分割し、同様にCatBoostモデルを学習させ、その性能を評価します。
```python 
# 26. フィルタリング後のデータで学習用とテスト用に分割:
# 300万円以下にフィルタリングされたデータセットを、再度学習用とテスト用に8:2の比率で分割する。
# 同様に「地名」で層化抽出し、目的変数には対数変換を適用する。
print("26. 特徴量と目的変数を分離、目的変数にlog1p変換を適用:")
train_df_filtered, test_df_filtered = train_test_split(
    filtered_data, test_size=0.2, stratify=filtered_data['chimei'], random_state=42
)

X_train_filtered = train_df_filtered.drop('kakaku_r07', axis=1)
y_train_filtered = np.log1p(train_df_filtered['kakaku_r07'])
X_test_filtered = test_df_filtered.drop('kakaku_r07', axis=1)
y_test_filtered = np.log1p(test_df_filtered['kakaku_r07'])
``` 
```text
26. 特徴量と目的変数を分離、目的変数にlog1p変換を適用
``` 
```python 
# 27. CatBoostモデルの学習と評価 (フィルタリング後のデータを使用):
# 地価300万円以下のフィルタリングされたデータセットを用いて、再度CatBoostモデルを学習・評価する。
# このモデルのRMSEを計算し、全体データで学習したモデルと比較することで、
# 特定の価格帯に特化したモデルの予測性能を検証する。
categorical_features_indices_filtered = [
    X_train_filtered.columns.get_loc(col) for col in X_train_filtered.select_dtypes(include='object').columns
]

model_filtered = CatBoostRegressor(iterations=100, random_state=42, verbose=0)
model_filtered.fit(X_train_filtered, y_train_filtered, cat_features=categorical_features_indices_filtered)

# 学習データでの予測とRMSEの計算（フィルタリング後、対数スケール）
y_train_pred_filtered = model_filtered.predict(X_train_filtered)
rmse_train_filtered = np.sqrt(mean_squared_error(y_train_filtered, y_train_pred_filtered))
print(f"300万円以下のデータのみで学習した CatBoost RMSE (学習データ): {rmse_train_filtered:.5f}") # 小数点以下5桁に調整

# テストデータでの予測とRMSEの計算（フィルタリング後、対数スケール）
y_pred_filtered = model_filtered.predict(X_test_filtered)
rmse_filtered = np.sqrt(mean_squared_error(y_test_filtered, y_pred_filtered))
print(f"300万円以下のデータのみで学習した CatBoost RMSE (テストデータ): {rmse_filtered:.5f}") # 小数点以下5桁に調整
print("-" * 50)
``` 
```text
300万円以下のデータのみで学習した CatBoost RMSE (学習データ): 0.02213
300万円以下のデータのみで学習した CatBoost RMSE (テストデータ): 0.04355
``` 
:::note info
300万円以下のデータに限定した場合のRMSE（対数スケール）は、**学習データで0.02213、テストデータで0.04355**となりました。
これは、**全体データで学習したモデルと非常に近い性能**です。学習データとテストデータのRMSEが近い値であり、過学習が起きていないことを示唆しています。
これにより、比較的データが豊富な価格帯に特化したモデルでも、高い汎化性能が期待できます。
:::

### フィルタリング後のモデルの最終評価
フィルタリング後のモデルについても、元のスケールでのRMSE、MAE、R²スコアを確認します。
```python 
# 28. フィルタリング後のデータでの元のスケールでのRMSE:
# 地価300万円以下のデータで学習したモデルの予測結果を元のスケールに戻し、
# 実際の円単位でのRMSEを計算・表示する。
# これにより、特定の価格帯に特化したモデルの予測誤差を実測値と比較して評価できる。
print("28. RMSE:")
y_test_original_filtered = np.expm1(y_test_filtered)
y_pred_original_filtered = np.expm1(y_pred_filtered)
rmse_original_scale_filtered = np.sqrt(
    mean_squared_error(y_test_original_filtered, y_pred_original_filtered)
)
# ここでround()で四捨五入する
pprint(f"300万円以下のデータのみで学習した CatBoost RMSE (元のスケール): {int(round(rmse_original_scale_filtered)):,} 円") # 四捨五入、カンマ区切りで表示
print("-" * 50)

# 29. フィルタリング後のデータでのMAEとR²スコア:
# 地価300万円以下のデータで学習したモデルの平均絶対誤差 (MAE) と決定係数 (R²スコア) を計算し表示する。
# これにより、特定の価格帯に特化したモデルの予測精度をより詳細に評価できる。
print("29. MAEとR2 スコア:")
mae_filtered = mean_absolute_error(y_test_filtered, y_pred_filtered)
r2_filtered = r2_score(y_test_filtered, y_pred_filtered)
print(f"300万円以下のデータのみで学習した CatBoost MAE: {mae_filtered:.5f}") # 小数点以下5桁
print(f"300万円以下のデータのみで学習した CatBoost R^2: {r2_filtered:.5f}") # 小数点以下5桁
print("-" * 50)
``` 
```text
28. RMSE:
300万円以下のデータのみで学習した CatBoost RMSE (元のスケール): 45,305 円
--------------------------------------------------
29. MAEとR² スコア:
300万円以下のデータのみで学習した CatBoost MAE: 0.03004
300万円以下のデータのみで学習した CatBoost R²: 0.99209
``` 

### フィルタリング後のモデルの評価:
:::note info
300万円以下のデータに限定した場合のRMSEは、**元のスケールで45,305円、R²は0.99209**となりました。これは、全体データで学習したモデル（RMSE 44,250円）と非常に近い性能です。この結果は、高価格帯の少数のデータが全体モデルの精度に大きな悪影響を与えていなかったことを示唆しています。

どちらのモデルも高い精度を達成しているため、分析の目的や対象とする価格帯に応じて、どちらのモデルを採用するかを柔軟に選択できることが分かりました。

**MAEは0.03004（対数スケール）と非常に小さく、R²スコアは0.99209と1に極めて近い値を示しています。R²スコアが0.99を超えることは、モデルが地価の変動の約99.2%を説明できていることを意味し、これは非常に高い予測精度であることを示しています。**
:::

### 予測結果の可視化 (フィルタリング後のデータ)
300万円以下のデータに限定して再構築したモデルについても、元のスケールでの実測値と予測値の散布図を作成し、視覚的にその性能を確認します。
```python 
# 30. 予測結果の可視化 (フィルタリング後のデータを使用):
# 地価300万円以下のデータで学習したモデルの予測結果を、元のスケールで散布図として可視化する。
# これにより、特定の価格帯におけるモデルの予測性能を視覚的に評価し、実際の値と予測値の一致度を確認できる。
plt.figure(figsize=(8, 8))
sns.scatterplot(x=y_test_original_filtered, y=y_pred_original_filtered)
plt.plot([y_test_original_filtered.min(), y_test_original_filtered.max()], [y_test_original_filtered.min(), y_test_original_filtered.max()], '--k')
plt.xlabel("実測値 (元のスケール)") 
plt.ylabel("予測値 (元のスケール)") 
plt.title("実測値 vs. 予測値 (300万円以下データ - 元のスケール)") 
plt.gca().xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, p: format(int(x), ',')))
plt.gca().yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, p: format(int(x), ',')))
plt.show()

print("-" * 50)
print("データ分析は以上です")
``` 
![#30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4091959/ca5bc783-0148-4f1f-9c3e-ed83c9617410.png)


:::note warn
データ分析は以上です
:::


この散布図から、**地価300万円以下のデータに限定したモデルも、非常に高い精度で地価を予測できていることが視覚的にも確認できます**。フィルタリング前のモデルと比べても、この価格帯における予測のばらつきが抑えられている様子が伺えます。


# 傾向
今回の地価予測モデルを構築する中で、地価形成に関するいくつかの明確な傾向が見えてきました。

:::note info
**最も顕著なのは、過去の価格データ (kakaku_r06、kakaku_r05) と公的な評価指標 (rosenka、kotei) が、地価予測において圧倒的に高い重要度を示した点です**。
これらの特徴量だけで、モデルの予測に約85%もの貢献をしていました。
これは、**地価が過去の市場動向や国・地方自治体による公式な評価に強く影響されるという、不動産評価の基本的な考え方をデータが裏付けていることを示しています**。
:::
:::note info
また、chimei_rank（地名ランキング）やeki_kyori（駅距離）といった地理的・地域的な特徴も、地価に一定の影響を与える要因であることが確認できました。
**特に、地名を平均地価に基づいてランキング化するという特徴量エンジニアリングの手法は、カテゴリカルな地域情報を数値化する有効なアプローチでした**。
:::
:::note info
そして、**CatBoostモデルの活用は、この分析において大きな成功要因となりました**。
CatBoostがカテゴリカル変数を直接扱う能力により、煩雑な手動での数値変換（例：One-Hot Encodingや順序エンコーディング）を行うことなく、複雑なカテゴリ間の関係性をモデルが自動的に学習し、非常に高い予測精度を達成できました。
これは、**データの前処理の手間を削減しつつ、モデルの性能を最大化する上で重要な知見となりました**。
:::

# 考察
:::note warn
結果としては、
**今回の地価予測モデルは、不動産鑑定評価を補完する簡易検証ツールとして非常に大きな可能性を秘めていると考えます**。
:::

### 客観的根拠の提供:
経験や専門性に依存しがちな不動産評価に対し、本モデルは大量の公的データに基づいた客観的な価格の目安を迅速に提供できます。
これにより、**鑑定評価の初期段階における市場価格のスクリーニングや、評価の妥当性を検証する上での客観的な根拠となり得ます**。

### 効率性の向上:
不動産鑑定士や不動産投資家が、多数の物件の概算価値を短時間で把握する必要がある場合、**本モデルは作業効率を大幅に向上させられると考えられます**。
特に、事前調査や簡易査定のツールとして活用することで、より専門的な鑑定作業にリソースを集中できるようになります。

# 課題
今回の分析で高い予測精度を達成できた一方で、いくつかの課題も見えてきました。

:::note info
まず、**データ数の不足による高価格帯での予測のばらつきが挙げられます**。
地価データは、比較的安価な物件のデータが豊富である一方、高額物件のデータは限られています。
これにより、**モデルが高価格帯の地価形成要因を十分に学習しきれず、予測の精度が低下する傾向が見られました**。

これを補うため、300万円以下のデータに限定したモデルも試みましたが、**根本的な解決にはデータ量の拡充が必要です**。
:::
:::note info
次に、**本モデルは「地価を予測する」ことには優れていますが、「地価がなぜその価格になるのか」という詳細な因果関係の解明には限界があります**。
例えば、rosenka（相続税路線価）やkotei（固定資産税路線価）が高い重要度を示したとしても、なぜその数値がその地域の地価に影響を与えるのか、といった深い理由まではモデル自体は教えてくれません。これは、**あくまで統計的な相関関係に基づいた予測であるためです**。
:::
:::note info
また、**利用した地価公示・地価調査データは公的なため高い信頼性がありますが、公開までに一定のタイムラグがあります**。
そのため、**モデルがリアルタイムな市場の変動や、直近の経済情勢を完全に反映した地価を予測することは難しいという課題も存在します**。
:::


# まとめ
Pythonと機械学習（CatBoost）を活用し、国土交通省の公的データを用いて東京23区住宅地の地価予測モデルを構築しました。

### 高い予測精度:
:::note info
CatBoostモデルにより、**元のスケールでRMSE約4.4万円、R²スコア約0.994という非常に高い予測精度を達成しました。
これは、モデルが地価の変動の約99.4%を説明できることを意味します**。
:::
### 重要な地価形成要因の特定: 
:::note info
過去の地価データ、相続税路線価、固定資産税路線価が地価予測に最も寄与する主要な要因であることを、特徴量重要度分析を通じて明確にしました。
:::

### CatBoostの有効性:
:::note info
カテゴリカル変数を直接扱うCatBoostの強みを活かすことで、複雑な手動エンコーディングなしに高精度を達成できることを実証しました。
:::
### 不動産鑑定評価の補完:
:::note info
本モデルは、**不動産鑑定評価の初期段階における客観的な参考価格を迅速に提供する簡易検証ツールとして、実務における客観性と効率性向上に貢献する可能性を秘めていることがわかりました**。

この簡易検証ツールの最も直接的な活用イメージは、Webアプリケーションとしての実装です。
ユーザーが基本的な土地情報（緯度・経度、地積、駅距離など）を入力するだけで、概算地価を即座に提示するツールとして機能します。
これは、不動産鑑定の初期段階でのスクリーニングや、市場価格の迅速な目安提供、さらには関連業務の効率化に貢献する**強力なツールとなり得ます**。
:::

今後は、より多様な地域への適用、時系列データの活用による価格変動予測、さらには地理空間情報や景観データといった非構造化データの組み込みなど、さらなるモデルの改善と応用可能性を探求していく予定です。

:::note warn
**この結果が、データ駆動型不動産評価の発展の一助となれば幸いです。**
:::

以上


# 出典
・国土交通省国土数値情報ダウンロードサイト https://nlftp.mlit.go.jp/
・国土地理院 https://www.gsi.go.jp/
・一般財団法人 資産評価システム研究センター（全国地価マップ） https://www.recpas.or.jp/
・警視庁（統計） https://www.npa.go.jp/index.html


:::note alert
**本モデルはあくまでデータに基づいた統計的予測であり、最終的な不動産鑑定評価に取って代わるものではありません**。
不動産鑑定評価には、個別具体的な物件の物理的特性、法的規制、周辺環境の将来性、さらには専門的な知見や経験に基づく総合的な判断が不可欠です。
:::


次の記事では、グリッドサーチを使って最適なハイパーパラメータを見つける方法を詳しく解説します。

https://qiita.com/samuraiyongt/items/5782e2b216d3c8673548