なんのコンペ？

皆さんはNishikaというデータ分析コンペについてご存じでしょうか？
かくいう私も、Udemyの講座を受け始めて興味を持ち、その講座内でNishikaのトレーニングコンペのデータをLightGBMを用いて分析したことをきっかけに参加してみようと決意しました。
Nishikaについて興味ある方は下記のURLから覗いてみてください。

さて、コンペの内容はと言いますと、画像のような中古マンションの価格予測です。
現在、2023年夏の部の中古マンションの価格予測のコンペも開催されているのでちょうどいいですね。
こちらの予測モデルを練習がてら作成し、開催されているコンペに参加してみようと思います！

データ分析で大事なこと

わたしがデータ分析を行う上で注意していること。
それは CRISP-DMのフレームワークに従う！ ということ。

CRISP-DMとは？

https://www.collidu.com/presentation-crisp-dm

Business Understanding（ビジネスの理解）
Data Understanding（データの理解）
Data Preparation（データの準備）
Modeling（モデリング）
Evaluation（評価）
Deployment（実装）
上記6つを受け取ったデータに対して実行していくことです。

どのようなデータがあるのか、データはどのような場所に保管されているのか、何を目的にデータを扱うのかデータについて理解を深めます。
データの状態（型や欠損値の有無など）を確認し、データを眺める
予測するために必要なデータに変換していく
データ予測に適したモデルを選択し、モデルを作成する
モデルの精度を評価する
実装

という流れになります。
画像にもあるようにビジネス理解とデータ理解、データの準備とモデリングは互いに相互関係にあるため、随時確認作業が必要です。
また、モデルの評価を行いますが、精度が悪かった場合は一番最初のビジネスの理解から振り返る必要があります。
今回はコンペの内容が決まってますのでデータを眺めるところから始めます。

データを見ていく

まずはこのトレーニング用のコンペに参加し、trainデータとtestデータをダウンロードします。
ではコードを書いていきましょう。
必要なライブラリをインポートします。

import glob
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

「%matplotlib inline」を追加することで

グラフがアウトプット行に出力される
plt.show()を省略できる
plt.show()でアウトプット行に2つ以上のグラフを表示可能

というメリットがあるのでとりあえず書いておきます。

ダウンロードしたtrainデータはtrainフォルダ内にたくさんのcsvファイルが格納されています。
このcsvファイルを一気に読み込みそれを結合しましょう。

files = glob.glob('train/train/*.csv')
data_list = []
for file in files:
  data_list.append(pd.read_csv(file,index_col=0))
df = pd.concat(data_list)

globライブラリを使用することで、指定したフォルダ内のファイルを一気に読み込むことができます。
また*はワイルドカードといい、すべて対象にとる的な認識でいいと思います。

では、データを見ていきましょう。

df.shape

データの数を見てみましょう。
（637351，27）
約64万行、27列のデータセットであることがわかります。
次にデータの情報についてみてみましょう。

df.info()

ちらほらと欠損値を多く含むカラムがありますね（non-nullが少ないほど欠損値が多い）。
また、「最寄り駅：距離（分）」など数値データでほしいカラムに文字列が入力されている状態です。
「市区町村コード」と「市区町村名」は同義ですよね。同じ相関を持つ特徴量を入れたままモデルを作成すると 多重共線性 が発生し精度に影響しますので片方消します。

ちょっとデータが汚いので綺麗にしますか。

nonnull_list = []
for col in df.columns:
    nonnull = df[col].count()
    if nonnull == 0:
        nonnull_list.append(col)
nonnull_list
df= df.drop(nonnull_list,axis=1)
df= df.drop('市区町村名',axis=1)
df = df.drop('種類',axis=1)

欠損値を含むカラムは削除できましたね。

しかし、「最寄り駅：距離（分）」や「面積（㎡）」、「築年数」、「取引時点」は数値としてほしい。
ではデータを見ていきましょう。

df['最寄駅：距離（分）'].value_counts()

なにやら怪しい文字列がありますね…
こいつらを数値に直していきましょう。

dis = {
  '30分?60分':45,
  '1H?1H30':75,
  '2H?':120,
  '1H30?2H':105
}
df['最寄駅：距離（分）'] = df['最寄駅：距離（分）'].replace(dis).astype(float)

怪しい値を辞書を用いてkeyに指定し、変換したい値をvalueに指定します。
あとはreplaceを用いて、型を指定してあげればデータを変換できます。

変換できていますね。
同様に「面積（㎡）」についてもデータを確認し、変換しましょう。

df['面積（㎡）'] = df['面積（㎡）'].replace('2000㎡以上',2000).astype(float)

「建築年」について、データも同様に見てみると、 平成19年 のような形で入力されています。
この数字の部分だけ抜き出してみましょう。

y_list = {}
for i in df['建築年'].value_counts().keys():
    if '平成' in i :
        num = float(i.split('平成')[1].split('年')[0])
        year = 33 - num
    if '令和' in i :
        num = float(i.split('令和')[1].split('年')[0])
        year = 3 - num
    if '昭和' in i :
        num = float(i.split('昭和')[1].split('年')[0])
        year = 96 - num
    y_list[i] = year
    y_list['戦前'] = 74
df['建築年'] = df['建築年'].replace(y_list)

最後に「取引時点」についても変換していきましょう。
データは 2018年第３四半期 のように入力されています。
この 年第３四半期 の部分を数値に変換しましょう。

year = {
  '年第１四半期':'.25',
  '年第２四半期':'.50',
  '年第３四半期':'.75',
  '年第４四半期':'.99'
}
year_list = {}
for i in df['取引時点'].value_counts().keys():
  for k,j in year.items():
    if k in i:
      year_rep = i.replace(k,j)
  year_list[i] = year_rep
year_list
df['取引時点'] = df['取引時点'].replace(year_list).astype(float)

変換できていますね。
あとは文字列のデータが入力されたカラムに対してカテゴリー変数に変換します。

for col in ["都道府県名", "地区名", "最寄駅：名称", "間取り", "建物の構造", "用途", "今後の利用目的", "都市計画", "改装", "取引の事情等"]:
    df[col] = df[col].astype("category")

では、データの統計値を見ていきましょう。

df.describe()

可視化していく

「最寄駅：距離（分）」、「面積（㎡）」、「建築年」、「取引価格（総額）_log」のヒストグラムを作成してみましょう。

fig,axes = plt.subplots(2,2,figsize=(20,10))
axes[0][0].hist(df['最寄駅：距離（分）'],bins=20)
axes[0][1].hist(df['面積（㎡）'],bins=200)
axes[0][1].set_xlim(0,250)
axes[1][0].hist(df['建築年'],bins=20)
axes[1][1].hist(df['取引価格（総額）_log'],bins=20)
plt.show()

それぞれのグラフから様々なことが読み取れます。
最寄り駅まで近ければ近いほど契約数が多い。
1人暮らしをしている人、ファミリー層が多い。
などが想像できます。
建築年は綺麗な正規分布となってますね。取引価格（総額）_logについては対数を取っていますので、実際にはここまで綺麗な正規分布ではないと思います。

今回のコンペでは「取引価格（総額）_log」がターゲットになっていますので、こいつに対する最寄駅：距離（分）」、「面積（㎡）」、「建築年」の散布図を作成します。

fig,axes = plt.subplots(3,1,figsize=(10,10))
axes[0].scatter(df['最寄駅：距離（分）'],df['取引価格（総額）_log'],alpha=0.1)
axes[1].scatter(df['面積（㎡）'],df['取引価格（総額）_log'],alpha=0.1)
axes[2].scatter(df['建築年'],df['取引価格（総額）_log'],alpha=0.1)
plt.show()

面積は若干ですが正の相関がありそうですね。
では、それぞれの相関係数を算出し、ヒートマップを作成してみましょう。

df[['取引価格（総額）_log','最寄駅：距離（分）','面積（㎡）','建築年']].corr()
import matplotlib
matplotlib.rcParams['font.family'] = 'AppleGothic'
sns.heatmap(df[['取引価格（総額）_log','最寄駅：距離（分）','面積（㎡）','建築年']].corr())

文字化けしてますごめんなさい。
matplotlib.rcParams['font.family'] = 'AppleGothic'を追加したんですけどうまくいきませんでした。
Japanize_matplotlibが必要か？
まあ、相関係数の表を一緒に載せるので許してください。
やはり、面積は若干ですが正の相関がありますね。
このとき、あまりにも相関の高い特徴量があった場合、それの影響が大きいためモデルの精度に影響が考えられるので特徴量から外すなどの処置が必要となる可能性もあることを頭の片隅に置いておいてください。
最後にここまでの処理を関数化しておきます。

def data_pre(df):
    nonnull_list = []
    for col in df.columns:
        nonnull = df[col].count()
        if nonnull == 0:
            nonnull_list.append(col)
    df = df.drop(nonnull_list, axis=1)

    df = df.drop("市区町村名", axis=1)

    df = df.drop("種類", axis=1)

    dis = {
        "30分?60分":45,
        "1H?1H30":75,
        "2H?":120,
        "1H30?2H":105
    }
    df["最寄駅：距離（分）"] = df["最寄駅：距離（分）"].replace(dis).astype(float)

    df["面積（㎡）"] = df["面積（㎡）"].replace("2000㎡以上", 2000).astype(float)


    y_list = {}
    for i in df["建築年"].value_counts().keys():
        if "平成" in i:
            num = float(i.split("平成")[1].split("年")[0])
            year = 33 - num
        if "令和" in i:
            num = float(i.split("令和")[1].split("年")[0])
            year = 3 - num
        if "昭和" in i:
            num = float(i.split("昭和")[1].split("年")[0])
            year = 96 - num
        y_list[i] = year
    y_list["戦前"] = 76
    df["建築年"] = df["建築年"].replace(y_list)

    year = {
        "年第１四半期": ".25",
        "年第２四半期": ".50",
        "年第３四半期": ".75",
        "年第４四半期": ".99"
    }
    year_list = {}
    for i in df["取引時点"].value_counts().keys():
        for k, j in year.items():
            if k in i:
                year_rep = i.replace(k, j)
        year_list[i] = year_rep
    df["取引時点"] = df["取引時点"].replace(year_list).astype(float)
    
    for col in ["都道府県名", "地区名", "最寄駅：名称", "間取り", "建物の構造", "用途", "今後の利用目的", "都市計画", "改装", "取引の事情等"]:
                df[col] = df[col].astype("category")
    
    return df

df = data_pre(df)

LightGBMを使ってモデルを作成

いよいよモデルを作成します。
今回は決定木モデルの勾配ブースティング手法であるLightGBMを用いていきます。
では早速構築していきましょう。

import lightgbm as lgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error as mae

df_train,df_val = train_test_split(df,test_size=0.2)

col = '取引価格（総額）_log'
train_y = df_train[col]
train_x = df_train.drop(col,axis=1)

val_y = df_val[col]
val_x = df_val.drop(col,axis=1)

trains = lgb.Dataset(train_x,train_y)
valids = lgb.Dataset(val_x,val_y)

params = {
  'objective':'regression',
  'metrics':'mae',
}

model = lgb.train(params, trains, valid_sets=valids, num_boost_round=1000, callbacks=[lgb.early_stopping(stopping_rounds=100)])

今回のコンペでは評価指標はMAE（平均絶対値誤差）と指定があるので、scikit-learnからインポートしてきましょう。
まずはデータ整形したデータフレームをトレーニングデータとバリデーションデータに分割します。
その後、トレーニングデータ、バリデーションデータにおいてターゲットと特徴量を指定します。
次にLightGBMのパラメータを設定します。Optunaなどを用いてハイパーパラメータのチューニングを行ったりしますが、今回はobjectiveとmetricsの指定だけします。
最後にモデルを実行します。ここではnum_boost_roundで決定木の更新回数、callbacksで精度が向上しなかった場合途中でモデルを停止する設定を行っています。

Did not meet early stopping. Best iteration is:
[993] valid_0's l1: 0.0765655

993回目の予測で精度が最も高く出たようですね。
一応、精度が正しいか確認しておきましょう。

vals = model.predict(val_x)
mae(vals, val_y)

0.0765655273764904

一致してますね。
LightGBMではターゲットに対しての特徴量の重要度を算出することができます。
確認してみましょう。

pd.DataFrame(model.feature_importance(), index=val_x.columns, columns=["importance"]).sort_values("importance", ascending=False)

地区名が最も重要度が高い結果となりました。
しかし、この結果が必ずしも予測に寄与しているとは限らないということは頭に入れておきましょう。
あとは、先ほど作成したモデルにダウンロードしておいたテストデータを入れて予測してみましょう。

df_test = pd.read_csv("test.csv", index_col=0)
df_test = data_pre(df_test)
predict = model.predict(df_test)
df_test["取引価格（総額）_log"] = predict
df_test[["取引価格（総額）_log"]].to_csv("submit_test.csv")

予測させたあとはcsvファイルとして出力し、コンペに提出します。
あとは結果を待つべし！！！

終わりに

お疲れさまでした！！！
長い時間お付き合いいただきありがとうございました。
今回はデータ分析におけるデータの理解、データの準備、モデリング、評価を行いました。
中身はどのような目的か、どのようなデータかなどによって異なってきます。
そのために、データ整形スキルや機械学習モデルについての勉強は必要ですね…
また、今回は紹介しませんでしたが、特徴量エンジニアリングや教師なし学習によって特徴量を生成したりなど、やれることはたくさんあります！
非常に難しい分野であると同時に奥深い分野でもあるので、おもしろいですね（めちゃくちゃ疲れるけど…）

データサイエンティストを目指して、日々精進してまいります！
また、予測モデルをつくったときは紹介します。
それでは！！！