[solafune] 夜間光データから土地価格を予測ベースモデル

# はじめに

衛星データに特化したデータ分析コンペサイト[solafune](https://solafune.com)にて、主題のコンペが2021/2/15~2021/4/5の間で開催されています。

衛星データのコンペと言っても画像データを扱うわけではなく、衛星によって取得された夜間光のデータ+地域ごとの地価データを1ファイル1MB以下のcsvファイルにしたデータを扱うもので、初心者にも非常に参加しやすいコンペとなっています。列数でいえばあのデータ分析の登竜門Titanic以下の4列（予測対象の1列除く）しかありませんし、国内コンペでは珍しく5位入賞まで賞金が出るおまけつきです。

この記事ではそんなコンペ開催を記念して、~~そして最もコンペに貢献したコンテンツを書いた人に贈られれるSolafune賞を狙って~~以下の様な方を対象にベースラインモデルを紹介します。

- Titanicはやったのでデータ分析コンペに参加してみたい！
- [Kaggleスタートブック本](https://www.kspub.co.jp/book/detail/5190067.html)以上[Kaggle本](https://gihyo.jp/book/2019/978-4-297-10843-4)未満程度の知識がある。
  - 記事中では例えば以下の様な言葉は説明なしに登場し、触れたとしてもごく少量の説明で各種ライブラリの使用等を行います。また、ipynbファイルでの作業を想定しており、ライブラリのインストール等を含めた環境構築も完了していることを前提としています。
    - 訓練（train）データ、評価（test）データ
    - EDA
    - 特徴量エンジニアリング
    - 交差検証
    - LightGBM
    - 評価指標
    - などなど...

[この記事よりもハイレベルな（そしてスコアが高い）ベースライン](https://qiita.com/tubo/items/f83a97f2488cc1f40088)も別の参加者の方が投稿されていますので、そちらも是非ご確認ください！

# ベースラインモデル

まずは何はともあれSolafuneへの登録、コンペへの参加登録、配布データのダウンロードを行っておいてください。ディレクトリ構成は以下を想定しています。

```
├── data
│   ├── TrainDataSet.csv       <- 配布データ1
│   ├── EvaluationData.csv     <- 配布データ2
│   └── UploadFileTemplate.csv <- 配布データ3
│
└── notebook                    <- .ipynbファイルの保存場所
```

## 基本ライブラリのインポート＆データ読み込み


```python
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme(style="ticks")
```


```python
# データを保存しているフォルダのディレクトリ指定
DATA_DIR = '../data/'

# データ読み込み
train = pd.read_csv(os.path.join(DATA_DIR, 'TrainDataSet.csv'))
test = pd.read_csv(os.path.join(DATA_DIR, 'EvaluationData.csv'))
submission = pd.read_csv(os.path.join(DATA_DIR, 'UploadFileTemplate.csv'))
```

1つ、補足します。`EvaluationData.csv`をtestとして読み込んでいますが、これはデータ分析コンペサイトの中でもSolafune独特のファイル名が付けられていることによります。evaluationと名付けられていますがこれが他で言うtestデータのことであり、交差検証の時に出てくるvalidationデータとは別物ですのでご注意ください。

簡単にデータのサイズや中身を見てみます。


```python
# shapeの表示
print(f'train     : {train.shape}')
print(f'test      : {test.shape}')
print(f'submission: {submission.shape}')
```

    train     : (21883, 5)
    test      : (13860, 4)
    submission: (13860, 3)



```python
# 先頭3行の表示
print('# train')
display(train.head(3))

print('# test')
display(test.head(3))

print('# submission')
display(submission.head(3))
```

    # train



<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }


    .dataframe tbody tr th {
        vertical-align: top;
    }
    
    .dataframe thead th {
        text-align: right;
    }

</style>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>PlaceID</th>
      <th>Year</th>
      <th>AverageLandPrice</th>
      <th>MeanLight</th>
      <th>SumLight</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1128</td>
      <td>1993</td>
      <td>740.909091</td>
      <td>57.571430</td>
      <td>403.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1128</td>
      <td>1994</td>
      <td>739.390909</td>
      <td>62.714287</td>
      <td>439.0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1128</td>
      <td>1995</td>
      <td>739.390909</td>
      <td>61.857143</td>
      <td>433.0</td>
    </tr>
  </tbody>
</table>

</div>


    # test



<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }


    .dataframe tbody tr th {
        vertical-align: top;
    }
    
    .dataframe thead th {
        text-align: right;
    }

</style>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>PlaceID</th>
      <th>Year</th>
      <th>MeanLight</th>
      <th>SumLight</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1387</td>
      <td>1992</td>
      <td>28.775862</td>
      <td>3338.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1387</td>
      <td>1993</td>
      <td>36.948277</td>
      <td>4286.0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1387</td>
      <td>1994</td>
      <td>30.706896</td>
      <td>3562.0</td>
    </tr>
  </tbody>
</table>

</div>


    # submission



<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }


    .dataframe tbody tr th {
        vertical-align: top;
    }
    
    .dataframe thead th {
        text-align: right;
    }

</style>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>PlaceID</th>
      <th>Year</th>
      <th>LandPrice</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1387</td>
      <td>1992</td>
      <td>0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1387</td>
      <td>1993</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1387</td>
      <td>1994</td>
      <td>0</td>
    </tr>
  </tbody>
</table>

</div>


trainにのみある`AverageLandPrice`が予測対象の列です。submissionの列名は`LandPrice`になっていますが特に意味はなさそうです…とにかくtrainで学習し、testの`AverageLandPrice`を予測、submissionの形式で提出することを目指します。

各列の説明をコンペサイトから引用しておきます。

| 列               | 内容             | 備考                                 |
| :--------------- | :--------------- | :----------------------------------- |
| PlaceID          | 地域ごとのデータ | 地域固有のID                         |
| Year             | 年代             | 1992~2013年まで                      |
| AverageLandPrice | 土地の平均価格   | 1992~2013年まで                      |
| MeanLight        | 夜間光量の平均値 | 0~63までのレンジでその地域の平均光量 |
| SumLight         | 夜間光量の合計値 | その地域の合計光量                   |

コンペサイトには各情報の引用元が記されており、`AverageLandPrice`については以下の通り説明されています。

> 地価公示（国土交通省）を元に弊社が集計、加工したものを使用しています。

公式に「この記載はデータが日本国内の地域に絞られていることを意味するか」と質問したところ「回答不可」とのことでしたが、引用元に海外の各地域の地価情報を示すものは無かったため、確証はありませんが日本国内の地域のデータに絞られている可能性が高いと勝手に思っています。

`MeanLight`と`SumLight`については引用元を辿ってもそれらしい記述が見当たらなかったのですが、同じデータを引用していると思われる[こちらの論文](https://dept.sophia.ac.jp/econ/econ_cms/wp-content/uploads/2016/11/62-2.pdf)に以下の記載があります。

> まず夜間光は、先行研究と同じく DMSP-OLS データを用いる（NGDC 2016）。ラスタ形式で配布されて
> いる同データは、各地域の 1km 四方のセルに対して夜間光の強さを 0（光が無い状態）から 63（最も光
> が強い状態）までの整数値で分類するものである。図 1 は 2010 年におけるバングラデシュ周辺の状況を
> 示しており、中央に位置するバングラデシュの首都ダッカや、左下に位置するインドの都市コルカタなど
> の地区で夜間光が強いことが確認できる。また、バングラデシュの 527 郡について夜間光の平均を集計し、
> その対数をヒストグラムとして示したのが図 2 である。なお同図が双峰性の分布を示しているのは、夜間
> 光量の最大値（63）がトップ・コーディングされた値であることに由来する。

`MeanLight`は恐らく上記の方法で算出した地域別の値を平均したもの、`SumLight`は足し合わせたものと思われます。となると、`SumLight`を`MeanLight`で割った値はその地域のおおよその面積を表しそうです。※面積が分かるからと言って地域を特定し、その情報を学習に使うのはコンペの規約違反（未承認外部データの利用）にあたりますのでご注意ください。

## EDA

もう少し詳しくデータを見ていきましょう。Kaggleスタートブック本に則って[pandas-profiling](https://github.com/pandas-profiling/pandas-profiling)やよく比較対象にあがる[sweetviz](https://github.com/fbdesignpro/sweetviz)を使用しても良いのですが、せっかくなので自力で見ていきます。

まずは基本のデータ型と欠損の確認から。


```python
train.info()
```

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 21883 entries, 0 to 21882
    Data columns (total 5 columns):
     #   Column            Non-Null Count  Dtype  
    ---  ------            --------------  -----  
     0   PlaceID           21883 non-null  int64  
     1   Year              21883 non-null  int64  
     2   AverageLandPrice  21883 non-null  float64
     3   MeanLight         21883 non-null  float64
     4   SumLight          21883 non-null  float64
    dtypes: float64(3), int64(2)
    memory usage: 854.9 KB



```python
test.info()
```

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 13860 entries, 0 to 13859
    Data columns (total 4 columns):
     #   Column     Non-Null Count  Dtype  
    ---  ------     --------------  -----  
     0   PlaceID    13860 non-null  int64  
     1   Year       13860 non-null  int64  
     2   MeanLight  13860 non-null  float64
     3   SumLight   13860 non-null  float64
    dtypes: float64(2), int64(2)
    memory usage: 433.2 KB


train、testともにすべての列がfloat64もしくはint64の数値型で、欠損がないことがわかりました。扱いやすいですね！

続いて、列毎に更に細かく見ていきます。

### AverageLandPrice

いきなりですが今回の予測対象であるこちらです。

```python
# 普通のヒストグラム
sns.histplot(data=train, x='AverageLandPrice', bins=100)
plt.show()
```
​![output_16_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/834c5a2a-7d1c-47dd-1f57-d8dd0eea745c.png)
ヒストグラムを書いてみましたが0～10,000以下？がほとんどみたいですがx軸が100,000以上まであって良く分からないグラフになりました。このままでは良く分からないので対数変換して可視化してみます。

```python
# 常用対数をとったヒストグラム
sns.histplot(data=train, x='AverageLandPrice', log_scale=10, bins=30)
plt.show()
```
![output_18_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/5e03a167-7222-d2b7-55f5-a4b4d6773b90.png)

良い感じに正規分布っぽくなりました。ほとんどが$10^3$程度以下の様ですが、数は少ないながらも大きいものだと$10^5$くらいまであるようです。※こちらの単位は公表されておらず、土地の何の平均価格かは不明です。恐らく単位千円の坪単価だとは思うのですが…確証はありません。


### ベン図

残りの列は一つ一つ見ていく前に、trainとtestで共通している列なのでどのくらいデータが被っているかをベン図を書いてみます。


```python
from matplotlib_venn import venn2
```


```python
# trainとtestの共通4列の要素がどれだけ被っているかをベン図で描画
fig, axes = plt.subplots(figsize=(16, 4), ncols=4)

for c, ax in zip(test.columns, axes.ravel()):
    venn2(
        subsets=(set(train[c].unique()), set(test[c].unique())),
        set_labels=('Train', 'Test'),
        ax=ax
    )
    ax.set_title(c)

fig.tight_layout()
plt.show()
```
![output_22_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/c48afd1e-1167-11e2-8c16-6ac2a36addba.png)


各データの数（ユニーク数）とともに、以下の様なことが分かります。

- PlaceIDはtrainとtestで全く被っていない。つまり同じ地域の情報を学習して予測するということができません。例えば東京都中央区のデータで学習して、別の地域である東京都千代田区の地価を予想するようなことが必要。（あくまで例であり、各`PlaceID`が具体的にどの地域を指しているかは不明です）
- Yearはすべて共通。列名の説明で記載した通り1992~2013年の22年間分。
- MeanLightはほとんど被っていない（float64の連続値なので当たり前）
- SumLightは半分以上は被っている（float64の連続値なのに…？偶然？）

以降ベン図以外も確認していきますが`AverageLandPrice`以外の列はtrainとtest両方に含まれますので、一つのDataFrameに結合しておき適宜いっぺんに確認していきたいと思います。また、結合時にはどの行がtrainのデータでどの行がtestのデータかが分かる様にフラグを付けておきます。


```python
# フラグ付け（trainには1、testには0）
train['is_train'] = 1
test['is_train'] = 0

# 結合（trainにしかない列は削除しつつ）
df = pd.concat([train.drop('AverageLandPrice', axis=1), test], ignore_index=True)
```

### PlaceID


```python
# is_trainフラグ毎にPlaceIDの最小値・四分位点・最大値を確認する
df.groupby('is_train')['PlaceID'].describe().iloc[:,3:]
```

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }


    .dataframe tbody tr th {
        vertical-align: top;
    }
    
    .dataframe thead th {
        text-align: right;
    }

</style>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>min</th>
      <th>25%</th>
      <th>50%</th>
      <th>75%</th>
      <th>max</th>
    </tr>
    <tr>
      <th>is_train</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>3.0</td>
      <td>446.0</td>
      <td>867.0</td>
      <td>1279.0</td>
      <td>1662.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1.0</td>
      <td>395.5</td>
      <td>810.0</td>
      <td>1223.0</td>
      <td>1660.0</td>
    </tr>
  </tbody>
</table>

</div>



train、testともに最小値・各四分位点・最大値の値がおおよそ一緒で満遍なく分割されている様です。

※`PlaceID`は数値データといってもただのIDなので、最大最小や四分位数に意味がある数値ではありません。しかしながら、極端な例ですとtrainには1~1000のID、testには1001以降のIDが割り振られている様な場合は何か意味のある分割（例えばもしかするとtrainは本州でtestは本州以外で分割されているかも、など）がされている可能性があるためdescribeを使って確認しています。またもちろん、今回の様にIDがばらついていても実は意味のある分割がされていることもあります。

### Year


```python
# trainとtestで色分けしてヒストグラムの描画
sns.histplot(data=df, x='Year', hue='is_train', bins=22)
plt.show()
```
![output_30_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/1f341974-3ebd-0351-895a-d398a2aaf022.png)
※積み上げ棒グラフの様に見えますがオレンジ(train)のヒストグラムの上に青(test)のヒストグラムが重なっています。
testはすべての年度が約600件ずつ均一にありそうですが、trainは少しだけ～1995年頃までの情報が少ない様です。

### MeanLight


```python
# trainとtestで色分けしてヒストグラムの描画
sns.histplot(data=df,x='MeanLight',hue='is_train',bins=64)
plt.show()
```
![output_33_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/9f5a7f31-fd0c-c249-552c-6faa7c71833f.png)
※こちらも積み上げ棒グラフの様に見えますが二つのヒストグラムが重なっています。
5近辺と最大値である63近辺に膨らみがある、双峰性の分布です。63が極端に多いのはデータの説明のところで触れた通りトップコーティングの影響と思われます。

### SumLight


```python
# trainとtestで色分けしてヒストグラムの描画
sns.histplot(data=df,x='SumLight',hue='is_train')
plt.show()
```
![output_36_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/29afb106-59c2-18ae-c34c-b82c54223f57.png)
ターゲットである`AverageLandPrice`と似たような形状をしており、対数変換すると正規分布に近い形になります。

### AverageLandPriceと各種列

※`AverageLandPrice`はtestには情報がないためtrainのデータのみを使っています。

#### with PlaceID


```python
sns.scatterplot(data=train, x='AverageLandPrice', y='PlaceID')
plt.show()
```
![output_40_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/70359007-15b2-70d0-905b-9881a8576993.png)
バラバラしています。やはり`PlaceID`はランダムに割り振られているようです。これは例えば「同じ都道府県なら近いID」などになっていた場合、「地価が高い東京地域のIDと高い地価がもっと密集した分布になるはずだ」という想定によります。図でいうと3つ`AverageLandPrice`が特に大きいものがありますが、これらは例えば東京都中央区や港区などの地価が高いエリアであると想定されます。

#### with Year

```python
sns.scatterplot(data=train, x='AverageLandPrice', y='Year')
plt.show()
```
![output_42_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/47cbb6c1-be7a-476a-de13-1bf8bef1e4f4.png)
影響の大小はありそうですが概ね年次が上がるほど下がる傾向にあるようです。1992年あたりはバブルだったことを考えると、当然でしょうか。

#### with MeanLight

```python
sns.scatterplot(data=train, x='MeanLight', y='AverageLandPrice')
plt.show()
```
​![output_45_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/015be68c-4272-fb40-c798-74e288e4f725.png)

完全ではありませんが傾向として`MeanLight`が大きいほど地価も高そうです。特に最大値63近辺で顕著です。夜間光量の平均値が高い（＝電気がギラギラ）なほど地価が高いというのは一般的な感覚とも合致します。

#### with SumLight

```python
sns.scatterplot(data=train, x='SumLight', y='AverageLandPrice')
plt.show()
```
![output_48_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/edb22f83-8aa8-11af-3353-3366bf36ba21.png)
今度はどちらかというと`MeanLight`とは逆の結果になりました。合計光量が少ないほど地価が高く、多いほど地価が安い傾向が見られます。`SumLight`が大きい->`MeanLight`が大きいとも考えられますが、別の場合を考えると平均的な光度は小さくとも単に面積が広い場合が考えられます。都市ほど地価は高くても面積が狭い傾向にありますので、その傾向の現れと思われます。

#### with area

データの概要で触れた通り、`SumLight`を`MeanLight`で割ったものはその地域の面積を表すと思われます。その値と`AverageLandPrice`の関係も見てみます。


```python
train['area'] = train['SumLight'] / train['MeanLight']
```


```python
sns.scatterplot(data=train, x='AverageLandPrice', y='area')
plt.show()
```
![output_53_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/21f7d98e-7767-ed6d-bd90-1018ea1e6d6e.png)
やはり基本的には面積が小さいほど地価が高い傾向にありそうです。ただし島村部など面積も地価も小さい箇所も混ざっているはずなので、必ずしもそうでもないため扱いに悩むところです。今回は実装しませんが、例えば`MeanLight`と合わせて考えてみることで「狭いけど栄えている街」を表す特徴量が作れるかも知れません。

## 特徴量エンジニアリング
### 方針
いきなり特徴量を作成する前に、EDAを通して見えてきたことをまとめつつ、どんな特徴量を作成するかまとめておきます。

#### PlaceID
- trainとtestで被っていないため、そのまま使用するのは危険そうです。
  - 特定できる方法で特徴量化してしまうと、交差検証の際に不当に高いスコアがでてしまう可能性があります。
    - 詳しくは交差検証の際に後述します。
  - `PlaceID`を特定せずに特徴量化する方法として、ここではCountEncodingで変換することにします。
    - ほとんどがデータがある年次分の値（=22）になってしまいますが、一部全年次分のデータがない地域があるため多少の差が付きます。
    - 反対に、仮にCountした数がバラバラですと結局CountEncodingをしてもその数の違いによって`PlaceID`を特定できてしまいます。CountEncodingが常に元の値を特定できないようにする変換方法ではないことに注意してください。


#### Year
- 影響の大小はありそうですが、傾向としてはバブルの影響か年次が上がるほど地価が下がる様なので、そのまま特徴量化します。
  - 今回使用するアルゴリズムは決定木系統のため、大小関係がついていれば変換の必要なく使用できます。
    - 決定木系統のアルゴリズムの場合は大小関係の変わらない変換は不要という点（と、そうとも限らない点）については、[こちらのブログ記事](https://amalog.hateblo.jp/entry/decision-tree-scaling)などをご参考ください。
  - 今回は実装しませんが、`Year`の影響を受けやすい地域とそうでない地域（バブルで地価が崩壊した地域とそこまで影響がなかった地域）を何かしらの方法で分類し、その分類別に年次をスコア化した特徴量なども効果的かもしれません。

#### MeanLight
- そのまま使います。
- いくつか、`PlaceID`単位で各種統計値を特徴量化します。
- また、63の値となる地域は特に地価が高い傾向にあったため、この情報も抜き出したいと思います。
  - 色々やり方あると思いますが、今回は「値が63になった回数」で表現したいと思います。
    - 22年間63を取り続ける地域と、何回か63もあるけど62...もある地域との差を表現したいと考え、一つの方法として選びました。

#### SumLight
- そのまま使います。`MeanLight`と`area`と情報が被っている気もしますが、今回使用するGBDT系のモデルは多少無意味な特徴量が含まれていても精度が落ちにくい傾向にあるため（詳細はkaggle本6.2などを参照ください）、とりあえず入れてみます。
- `MeanLight同様`、いくつか`PlaceID`単位で各種統計値を特徴量化します。

#### area
- `SumLight`を`MeanLight`で割って算出します。
- 面積のはずなのに、同じ地域でも小数点以下程度でズレが生じるので地域ごとの中央値に置き換えます（おそらく配布データの桁数の問題？）。

### 特徴量作成
ではいよいよ特徴量を作っていきましょう！
まずは、特徴量作成に便利なライブラリをimportしつつ、EDAの時に色々といじくって余計な情報を加えてしまったのでtrainとtestを読み込みなおしておき、比較的簡単な特徴量を作成していきます。このとき、今回はtrainとtestを別々に処理してもまとめて処理しても結果が変わらない特徴量しか作成していないため、コードを短く済ますために最初にtrainとtestを結合してまとめて特徴量を作成しています。しかし場合によってはtrainとtestは別々のまま作るべき特徴量もあれば、まとめて作るべき特徴量もあります。

※category_encodersの使い方は[こちらのQiita記事](https://qiita.com/Hyperion13fleet/items/afa49a84bd5db65ffc31)などを参照ください。


```python
import category_encoders as ce # CountEncodingなどカテゴリ変換に便利なライブラリ

train = pd.read_csv(os.path.join(DATA_DIR, 'TrainDataSet.csv'))
test = pd.read_csv(os.path.join(DATA_DIR, 'EvaluationData.csv'))
df = pd.concat([train.drop('AverageLandPrice', axis=1), test], ignore_index=True)
```


```python
# PlaceIDのCountEncoding
enc = ce.CountEncoder(cols='PlaceID')
df['count_PlaceID'] = enc.fit_transform(df)['PlaceID']

# 面積の計算
df['area'] = df['SumLight'] / df['MeanLight']

# 面積なのに同じ地域でも年次によって小数点以下でズレが生じるため、地域ごとの中央値に置き換える
_mapping = df.groupby('PlaceID')['area'].median()
df['area'] = df['PlaceID'].map(_mapping)

# 面積について、SumLightとMeanLightがともに0な地域があり、NaNになってしまうので0で補完する
# 完全に光が無い闇の地域？引用元でのデータ取得時のエラーが原因？
df['area'] = df['area'].fillna(0)

# 各地域でMeanLightが63をとった回数を特徴量にする
_mapping = df[df['MeanLight']==63].groupby('PlaceID').size()
df['count63'] = df['PlaceID'].map(_mapping)

# MeanLightが一回も63になっていないとNaNになってしまうので0で補完する
df['count63'] = df['count63'].fillna(0)
```
---
#### 2/19 19:26追記

上記のコードで`area`と`count63`の欠損値を0で補完していますが、LightGBMは欠損値があっても処理をして訓練・予測を行ってくれます。「LightGBMも欠損値処理しないと訓練できないのか」と誤解されませぬと幸いです。「その処理とは何ぞや」については大変申し上げにくいのですが私も完璧には理解できておらず、英語ですがこちらの[LightGBMのissue](https://github.com/microsoft/LightGBM/issues/2921)が参考になるかと思います。

とてもざっくり言うと欠損値があってもエラーを出さないように欠損値なりの分岐をするようになっており、今回のコードで補完している箇所の値を変えると以下の様にCVスコア（後述の予測の項で算出しています）が変化します。黄色の両方0で補完しているのが本記事での結果です。
![無題.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/b1f4639c-df01-dfda-3ebc-35652a0c1f84.png)

上記の通り実は両方0で補完した場合とNaNで補完した結果は変わらないのですが他の組み合わせでよりCVスコアが良いものもあります（あくまで今回の状況においてはこうなったという話で常にこうなるわけではありません）。じゃあなんでお前はわざわざ0で補完したのかというと…

- `count63`：`MeanLight`が63になった回数なので「0回」と明示的に表すことに意味があると思った
- `area`：0除算なので計算できないが、`MeanLight`と`SumLight`が両方0なのはなにかエラーが起きていそう->元となった数値のエラーを0としている（と思われる）なら、その結果も0にしておくと良いのかも知れない

と考えたことによります。コンペにおいて上位に入ることを目指すのであれば、（Privateと相関していると信じられる）交差検証の結果を見ながら処理を決めるのが良いかと考えています。

---
これでシンプルな特徴量はできましたが、少し作成に手間がかかる特徴量も作ってみます。冒頭のメモで触れていた「`PlaceID`単位の`MeanLight`と`SumLight`の各種統計値」というもので、pandasのgroupbyメソッドを使って作成します。基本的な統計量はデフォルトの関数を利用できますが、少し凝った統計量も関数を自作して適用します。


```python
# 最小値と最大値の差
def range_diff(x):
    return x.max() - x.min()

# 第一/三四分位点とその差
def third_quartile(x):
    return x.quantile(0.75)
def first_quartile(x):
    return x.quantile(0.25)
def quartile_range(x):
    return x.quantile(0.75) - x.quantile(0.25)
```


```python
# 集計
cols = ['MeanLight', 'SumLight']
aggs = ['min', 'mean', 'median', 'max', 'std', 'skew',
        range_diff, third_quartile, first_quartile, quartile_range]
_df = df.groupby('PlaceID')[cols].agg(aggs).reset_index()

# 列名変更のため、aggsの各名称を_aggsに格納
_aggs = []
for agg in aggs:
    if isinstance(agg, str):
        _aggs.append(agg)
    else:
        _aggs.append(agg.__name__)
        
# 列名変更
_df.columns = ['PlaceID'] + ["_".join([c, agg]) for c in cols for agg in _aggs]
```

ここまでの処理で、以下の様なPlaceID毎の`MeanLight`、`SumLight`の各種統計量をまとめたDataFrameができました。あとはこれを元のDataFrameに`PlaceID`をキーにmergeします。


```python
_df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }


    .dataframe tbody tr th {
        vertical-align: top;
    }
    
    .dataframe thead th {
        text-align: right;
    }

</style>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>PlaceID</th>
      <th>MeanLight_min</th>
      <th>MeanLight_mean</th>
      <th>MeanLight_median</th>
      <th>MeanLight_max</th>
      <th>MeanLight_std</th>
      <th>MeanLight_skew</th>
      <th>MeanLight_range_diff</th>
      <th>MeanLight_third_quartile</th>
      <th>MeanLight_first_quartile</th>
      <th>...</th>
      <th>SumLight_min</th>
      <th>SumLight_mean</th>
      <th>SumLight_median</th>
      <th>SumLight_max</th>
      <th>SumLight_std</th>
      <th>SumLight_skew</th>
      <th>SumLight_range_diff</th>
      <th>SumLight_third_quartile</th>
      <th>SumLight_first_quartile</th>
      <th>SumLight_quartile_range</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>6.654412</td>
      <td>9.578523</td>
      <td>9.625000</td>
      <td>12.708955</td>
      <td>1.219428</td>
      <td>0.101592</td>
      <td>6.054543</td>
      <td>10.244485</td>
      <td>9.016956</td>
      <td>...</td>
      <td>905.0</td>
      <td>1299.818182</td>
      <td>1305.0</td>
      <td>1703.0</td>
      <td>163.160055</td>
      <td>0.004118</td>
      <td>798.0</td>
      <td>1393.25</td>
      <td>1221.25</td>
      <td>172.00</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2</td>
      <td>41.154762</td>
      <td>48.058121</td>
      <td>48.285715</td>
      <td>53.911243</td>
      <td>3.037256</td>
      <td>-0.471909</td>
      <td>12.756481</td>
      <td>49.413692</td>
      <td>47.144344</td>
      <td>...</td>
      <td>6914.0</td>
      <td>8080.590909</td>
      <td>8112.0</td>
      <td>9111.0</td>
      <td>515.614442</td>
      <td>-0.436319</td>
      <td>2197.0</td>
      <td>8301.50</td>
      <td>7920.25</td>
      <td>381.25</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3</td>
      <td>6.647059</td>
      <td>9.794118</td>
      <td>9.808824</td>
      <td>12.735294</td>
      <td>1.582701</td>
      <td>-0.064804</td>
      <td>6.088235</td>
      <td>10.963236</td>
      <td>8.669117</td>
      <td>...</td>
      <td>226.0</td>
      <td>333.000000</td>
      <td>333.5</td>
      <td>433.0</td>
      <td>53.811842</td>
      <td>-0.064804</td>
      <td>207.0</td>
      <td>372.75</td>
      <td>294.75</td>
      <td>78.00</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4</td>
      <td>28.584158</td>
      <td>36.124792</td>
      <td>35.789603</td>
      <td>44.330048</td>
      <td>3.660401</td>
      <td>-0.062634</td>
      <td>15.745890</td>
      <td>38.267328</td>
      <td>34.143565</td>
      <td>...</td>
      <td>5774.0</td>
      <td>7302.681818</td>
      <td>7229.5</td>
      <td>8999.0</td>
      <td>746.146761</td>
      <td>-0.032737</td>
      <td>3225.0</td>
      <td>7730.00</td>
      <td>6897.00</td>
      <td>833.00</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5</td>
      <td>41.631580</td>
      <td>49.547581</td>
      <td>49.368420</td>
      <td>55.777779</td>
      <td>3.309324</td>
      <td>-0.578802</td>
      <td>14.146199</td>
      <td>51.868420</td>
      <td>48.381578</td>
      <td>...</td>
      <td>791.0</td>
      <td>934.181818</td>
      <td>935.5</td>
      <td>1027.0</td>
      <td>57.405522</td>
      <td>-0.777810</td>
      <td>236.0</td>
      <td>972.50</td>
      <td>917.50</td>
      <td>55.00</td>
    </tr>
  </tbody>
</table>
<p>5 rows × 21 columns</p>

</div>




```python
# 元のDataFrameにmerge
df = df.merge(_df, on='PlaceID', how='left')
```

## 訓練・予測
いよいよ訓練＆予測です！ここまでに比べると長いコードですがkaggleスタートアップ本を参考に書いております。自分もコンペ初参加時はとても参考にした本ですので、お持ちでない方は是非是非読んでみてください！
ここではlightgbmの説明は省かせて頂きますが、2点補足いたします。

### y_trainを対数変換する理由
lightgbmのobjectiveには今回の評価指標であるRMSLEがなく、最適化するためには自分で関数を作成するなどの対応をとる必要があります。それはちょっと面倒なので「対数変換したものをy_trainとしてRMSEで最適化して、あとで元のスケールに戻す」という方法が良く取られます。対数変換せずにいきなりRMSEで最適化するのがよろしくない理由は[こちらの記事](https://aizine.ai/rmse-rmsle1114/)を、具体的にどうやって変換するかは以下のコードの通りですが[こちらのQiita記事](https://qiita.com/shnchr/items/a88b1eb2bba8002ecd5d)などもご参考ください。

### バリデーション方法をGroupKFoldにしている理由
EDAで見た通り、trainとtestでは`PlaceID`が一切被っていません。順当に考えれば`PlaceID`が同じ場合の地価は似たような値になるはずで、例えば訓練時に「各地域の1992年~～2010年の情報を使って学習して2011年～2013年を予測し評価する」としてしまうとtestの時にはできないカンニングをしていることになり、不当に高いスコアとなり正当な評価ができなくなってしまいます。そのため、できるだけ条件をtestを予測する際と合わせるために訓練データと検証データで`PlaceID`が混ざらないようにGroupKFoldを選んでいます。


```python
x_train = df[:len(train)].drop('PlaceID', axis=1).reset_index(drop=True)
x_test = df[len(train):].drop('PlaceID', axis=1).reset_index(drop=True)
y_train = np.log1p(train['AverageLandPrice'])
```


```python
import lightgbm as lgb
from sklearn.model_selection import GroupKFold
```


```python
y_preds = []
models = []
oof_train = np.zeros((len(x_train),))
cv = GroupKFold(n_splits=5)

params = {
    'objective': 'rmse',
    'learning_rate': 0.01,
}

for fold_id, (train_index, valid_index) in enumerate(cv.split(train, y_train, train['PlaceID'])):
    X_tr = x_train.loc[train_index, :]
    X_val = x_train.loc[valid_index, :]
    y_tr = y_train.loc[train_index]
    y_val = y_train.loc[valid_index]
    
    lgb_train = lgb.Dataset(X_tr, y_tr)
    lgb_eval = lgb.Dataset(X_val, y_val, reference=lgb_train)
    
    model = lgb.train(params, lgb_train,
                      valid_sets=[lgb_train, lgb_eval],
                      verbose_eval=100,
                      num_boost_round=10000,
                      early_stopping_rounds=100)
    
    oof_train[valid_index] = \
        model.predict(X_val, num_iterations=model.best_iteration)
    y_pred = model.predict(x_test, num_iterations=model.best_iteration)
    
    y_preds.append(y_pred)
    models.append(model)
```

    Training until validation scores don't improve for 100 rounds
    [100]	training's rmse: 0.615957	valid_1's rmse: 0.646416
    [200]	training's rmse: 0.445885	valid_1's rmse: 0.578616
    [300]	training's rmse: 0.371811	valid_1's rmse: 0.575487
    Early stopping, best iteration is:
    [270]	training's rmse: 0.389244	valid_1's rmse: 0.574996
    Training until validation scores don't improve for 100 rounds
    [100]	training's rmse: 0.619655	valid_1's rmse: 0.675378
    [200]	training's rmse: 0.451198	valid_1's rmse: 0.592721
    [300]	training's rmse: 0.377045	valid_1's rmse: 0.574313
    [400]	training's rmse: 0.333679	valid_1's rmse: 0.571923
    [500]	training's rmse: 0.300587	valid_1's rmse: 0.572853
    Early stopping, best iteration is:
    [460]	training's rmse: 0.313101	valid_1's rmse: 0.57166
    Training until validation scores don't improve for 100 rounds
    [100]	training's rmse: 0.60571	valid_1's rmse: 0.773037
    [200]	training's rmse: 0.439709	valid_1's rmse: 0.67063
    [300]	training's rmse: 0.368091	valid_1's rmse: 0.652486
    [400]	training's rmse: 0.324932	valid_1's rmse: 0.649166
    Early stopping, best iteration is:
    [393]	training's rmse: 0.327592	valid_1's rmse: 0.648902
    Training until validation scores don't improve for 100 rounds
    [100]	training's rmse: 0.613064	valid_1's rmse: 0.687336
    [200]	training's rmse: 0.446505	valid_1's rmse: 0.585084
    [300]	training's rmse: 0.372544	valid_1's rmse: 0.5676
    [400]	training's rmse: 0.329101	valid_1's rmse: 0.568454
    Early stopping, best iteration is:
    [329]	training's rmse: 0.358066	valid_1's rmse: 0.566874
    Training until validation scores don't improve for 100 rounds
    [100]	training's rmse: 0.607628	valid_1's rmse: 0.725103
    [200]	training's rmse: 0.438502	valid_1's rmse: 0.643018
    [300]	training's rmse: 0.36519	valid_1's rmse: 0.633235
    [400]	training's rmse: 0.322654	valid_1's rmse: 0.630308
    [500]	training's rmse: 0.289557	valid_1's rmse: 0.629517
    [600]	training's rmse: 0.2642	valid_1's rmse: 0.630763
    Early stopping, best iteration is:
    [510]	training's rmse: 0.286657	valid_1's rmse: 0.629432



```python
scores = [
    m.best_score['valid_1']['rmse'] for m in models
]
score = sum(scores) / len(scores)
print('===CV scores===')
print(scores)
print(score)
```

    ===CV scores===
    [0.5749961295452265, 0.5716597018095386, 0.648901926064152, 0.5668740744440104, 0.6294317310106663]
    0.5983727125747187


ここまではkaggleスタートアップ本の通りです。スタートアップ本の分類の例と異なり今回は回帰問題ですので、各modelの予測結果の平均を最終予測とします。※y_trainを対数変換したのでy_predも対数スケールでの予測結果です。スケールを戻すのを忘れずに！


```python
y_pred = 0
for pred in y_preds:
    y_pred += np.expm1(pred) / len(y_preds)
```

投稿の前にあまりに的外れな予測をしていないか、正解であるy_trainとout-of-foldの予測結果が似たような分布になっているか確認してみます。


```python
fig, ax = plt.subplots(figsize=(8, 8))
sns.histplot(y_train, label='y_train', kde=True ,stat="density", common_norm=False, color='orange', alpha=0.3)
sns.histplot(oof_train, label='oof_train', kde=True, stat="density", common_norm=False, alpha=0.3)
ax.legend()
ax.grid()
```
![output_79_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/786297/dcec2d2d-2269-9739-8ba6-7398bba69f4e.png)
予測結果の方が真ん中に寄った分布になってしまってはいますが、極端なズレはないためひとまずそれらしい予測はできたと見て予測ファイルを作成し投稿します。ちなみにこちらの結果を投稿するとPublicLBで0.551596のスコアになります。


```python
# 予測ファイルの作成
submission['LandPrice'] = y_pred
submission.to_csv('submission.csv', index=False)
```

# さいごに

以上、予測を作成し投稿する一連の流れの一例を紹介しました。データの捉え方、特徴量の作り方、アルゴリズムの選び方などなど…これだけが正解ではもちろんないですが、この記事がデータ分析コンペ参加にハードルを感じている方の後押しとなれば幸いです。