[solafune] 夜間光データから土地価格を予測 BaseLine

20 min read

はじめに

こんにちは。mst8823という名前でいろんなコンペに参加している者です。
この記事は、コンペプラットフォームsolafuneで行われている夜間光コンペのベースラインとなります。(CV: 0.5353, LB: 0.4918)

directory structure

ディレクトリ構成は以下を想定しています。

  • input: TrainDataSet.csv などの配布データ置き場
  • notebooks: ベースラインノートブック置き場
  • submission: 提出用ファイル置き場
├─input
├─notebooks
└─submission

import

最初に必要なライブラリをインポートします

import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from lightgbm import LGBMModel
from sklearn import model_selection
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from tqdm import tqdm

config

INPUT = "../input/"
warnings.filterwarnings('ignore')
TARGET = "AverageLandPrice"

load data

使うデータたちを読み込みます

train = pd.read_csv(INPUT+"TrainDataSet.csv")
test = pd.read_csv(INPUT+"EvaluationData.csv")
sample_sub = pd.read_csv(INPUT+"UploadFileTemplate.csv")

feature engineering

前処理 & 特徴量作成をします
はじめに特徴量作成で使用する汎用的な関数を定義しておきます

# aggregationのagg_methodsで使用
def max_min(x):
    return x.max()-x.min()

def q75_q25(x):
    return x.quantile(0.75) - x.quantile(0.25)

# xfeat.aggregationとほぼ同じ
def aggregation(input_df, group_key, group_values, agg_methods):
    new_df = []
    for agg_method in agg_methods:
        for col in group_values:
            if callable(agg_method):
                agg_method_name = agg_method.__name__
            else:
                agg_method_name = agg_method
            new_col = f"agg_{agg_method_name}_{col}_grpby_{group_key}"
            df_agg = (input_df[[col] + [group_key]].groupby(group_key)[[col]].agg(agg_method))
            df_agg.columns = [new_col]
            new_df.append(df_agg)
            
    _df = pd.concat(new_df, axis=1).reset_index()
    output_df = pd.merge(input_df[[group_key]], _df, on=group_key, how="left")
    return output_df.drop(group_key, axis=1)

# group 内で diffをとる関数
def diff_aggregation(input_df, group_key, group_values, num_diffs):
    dfs = []
    for nd in num_diffs:
        _df = input_df.groupby(group_key)[group_values].diff(nd)
        _df.columns = [f'diff={nd}_{col}_grpby_{group_key}' for col in group_values]
        dfs.append(_df)
    output_df = pd.concat(dfs, axis=1)
    return output_df

# group 内で shiftをとる関数
def shift_aggregation(input_df, group_key, group_values, num_shifts):
    dfs = []
    for ns in num_shifts:
        _df = input_df.groupby(group_key)[group_values].shift(ns)
        _df.columns = [f'shift={ns}_{col}_grpby_{group_key}' for col in group_values]
        dfs.append(_df)
    output_df = pd.concat(dfs, axis=1)
    return output_df

つぎに実際に特徴量を作成する関数たちを定義します

# そのままの値の特徴量
def get_raw_features(input_df):
    cols = [
        "MeanLight",
        "SumLight",
        "Year"
    ]
    return input_df[cols].copy()

# 面積(?)
def get_area_feature(input_df):
    output_df = pd.DataFrame()
    output_df["Area"] = input_df["SumLight"] / (input_df["MeanLight"]+1e-3)
    return output_df

# PlaceIDをキーにした集約特徴量
def get_agg_place_id_features(input_df):
    _input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)
    group_key = "PlaceID"
    group_values = ["MeanLight", "SumLight", "Area"]
    agg_methods = ["min", "max", "median", "mean", "std", max_min, q75_q25]
    output_df = aggregation(_input_df, 
                            group_key=group_key, 
                            group_values=group_values,
                            agg_methods=agg_methods)
    return output_df

# Year をキーにした集約特徴量
def get_agg_year_features(input_df):
    _input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)   
    group_key = "Year"
    group_values = ["MeanLight", "SumLight", "Area"]
    agg_methods = ["min", "max", "median", "mean", "std", max_min, q75_q25]
    output_df = aggregation(_input_df, 
                            group_key=group_key, 
                            group_values=group_values,
                            agg_methods=agg_methods)
    return output_df

# PlaceID をキーにしたグループ内差分
def get_diff_agg_place_id_features(input_df):
    group_key = "PlaceID"
    group_values = ["MeanLight", "SumLight"]
    num_diffs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
    output_df = diff_aggregation(input_df, 
                                 group_key=group_key, 
                                 group_values=group_values, 
                                 num_diffs=num_diffs)
    return output_df

# PlaceID をキーにしたグループ内シフト
def get_shift_agg_place_id_features(input_df):
    group_key = "PlaceID"
    group_values = ["MeanLight", "SumLight"]
    num_shifts = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
    output_df = shift_aggregation(input_df, 
                                  group_key=group_key, 
                                  group_values=group_values, 
                                  num_shifts=num_shifts)
    return output_df

# pivot tabel を用いた特徴量
def get_place_id_vecs_features(input_df):
    _input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)
    # pivot table
    area_df = pd.pivot_table(_input_df, index="PlaceID", columns="Year", values="Area").add_prefix("Area=")
    mean_light_df = pd.pivot_table(_input_df, index="PlaceID", columns="Year", values="MeanLight").add_prefix("MeanLight=")
    sum_light_df = pd.pivot_table(_input_df, index="PlaceID", columns="Year", values="SumLight").add_prefix("SumLight=")
    all_df = pd.concat([area_df, mean_light_df, sum_light_df], axis=1)
    
    # PCA all 
    sc_all_df = StandardScaler().fit_transform(all_df.fillna(0))
    pca = PCA(n_components=64, random_state=2021)
    pca_all_df = pd.DataFrame(pca.fit_transform(sc_all_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_all_PCA={x:03}")
    # PCA Area
    sc_area_df = StandardScaler().fit_transform(area_df.fillna(0))
    pca = PCA(n_components=16, random_state=2021)
    pca_area_df = pd.DataFrame(pca.fit_transform(sc_area_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_Area_PCA={x:03}")
    # PCA MeanLight
    sc_mean_light_df = StandardScaler().fit_transform(mean_light_df.fillna(0))
    pca = PCA(n_components=16, random_state=2021)
    pca_mean_light_df = pd.DataFrame(pca.fit_transform(sc_mean_light_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_MeanLight_PCA={x:03}")
    # PCA SumLight
    sc_sum_light_df = StandardScaler().fit_transform(sum_light_df.fillna(0))
    pca = PCA(n_components=16, random_state=2021)
    pca_sum_light_df = pd.DataFrame(pca.fit_transform(sc_sum_light_df), index=all_df.index).rename(columns=lambda x: f"PlaceID_SumLight_PCA={x:03}")
    
    df = pd.concat([all_df, pca_all_df, pca_area_df, pca_mean_light_df, pca_sum_light_df], axis=1)
    output_df = pd.merge(_input_df[["PlaceID"]], df, left_on="PlaceID", right_index=True, how="left")
    return output_df.drop("PlaceID", axis=1)

# PlaceIDをキーにしたグループ内相関係数
def get_corr_features(input_df):
    _input_df = pd.concat([input_df, get_area_feature(input_df)], axis=1)
    group_key = "PlaceID"
    group_vlaues = [
        ["Year", "MeanLight"],
        ["Year", "SumLight"],
        ["Year", "Area"],
    ]
    dfs = []
    for gv in group_vlaues:
        _df = _input_df.groupby(group_key)[gv].corr().unstack().iloc[:, 1].rename(f"Corr={gv[0]}-{gv[1]}")
        dfs.append(pd.DataFrame(_df))
    dfs = pd.concat(dfs, axis=1)
    output_df = pd.merge(_input_df[[group_key]], dfs, left_on=group_key, right_index=True, how="left").drop(group_key, axis=1)
    return output_df
    
# count 63
def get_count63_feature(input_df):
    # 各地域でMeanLightが63をとった回数を特徴量にする
    _mapping = input_df[input_df['MeanLight']==63].groupby('PlaceID').size()
    
    output_df = pd.DataFrame()
    output_df['count63'] = input_df['PlaceID'].map(_mapping).fillna(0)
    return output_df
  • 今回提供されているデータは、学習用データ・テストデータともにほぼ完全なパネルデータの構造をしています(学習用データには欠損がある場合がある)。
    また、学習用データとテストデータには同じPlaceIDが存在しません。
  • このことから、あるPlaceID、つまり特定の場所の特性を特徴量として入れることが重要になりそうです。イメージとしてはPlaceIDごとの分散表現のようなものを与える感じでしょうか。
    そのイメージで拵えた特徴量がplace_id_vecs_featuresです。
  • この特徴量はPlaceIDYearで作成したピボットテーブルが主です。また、それをPCAで次元削減したものです。こうすることで、場所ごとで同一の特徴量を分散表現チックに獲得するとこができるのではないかという思惑です。
  • PlaceIDをキーとしたグループ内の相関係数を特徴量としているのも同様の考えです。

pivot tableで作った特徴量はこんな感じです。これをPlaceIDでマージします

Year	MeanLight=1992	MeanLight=1993	MeanLight=1994	MeanLight=1995	MeanLight=1996
PlaceID					
1	10.867647	8.985294	10.595589	10.794118	10.404411
2	48.815475	42.107143	50.172619	48.172619	48.398811
4	33.886139	29.480198	35.683167	35.009899	35.272278
8	32.555557	31.166666	35.027779	34.083332	33.861111
10	3.800000	5.400000	7.400000	5.400000	5.800000

最後に上記で定義した関数を実行します。
ここでは学習用データとテストデータを結合してから前処理していますが、気になる方は別々の処理に書き直してください<(_ _)>(PCA等の処理で少しleakyになります)

# 特徴量作成用の関数を実行する関数
def preprocess(train, test):
    input_df = pd.concat([train, test]).reset_index(drop=True)  # use concat data

    # load process functions
    process_blocks = [
        get_raw_features,
        get_area_feature,
        get_agg_place_id_features,
        get_agg_year_features,
        get_diff_agg_place_id_features,
        get_shift_agg_place_id_features,
        get_place_id_vecs_features,
        get_corr_features,
        get_count63_feature
    ]

    # preprocess
    output_df = []
    for func in tqdm(process_blocks):
        _df = func(input_df)
        output_df.append(_df)
    output_df = pd.concat(output_df, axis=1)

    # separate train and test
    train_x = output_df.iloc[:len(train)]
    test_x = output_df.iloc[len(train):].reset_index(drop=True)
    train_y = train[TARGET]
    
    return train_x, train_y, test_x
    
train_x, train_y, test_x = preprocess(train, test)
train_y = np.log1p(train_y)  # 目的変数のlogをとる

train & predict

最初に使うモデルの定義とパラメータの設定を行います
今回はお馴染みlightGBMです

model = LGBMModel  # 使うモデル
group = train["PlaceID"]  # 今回は group kfold を用いる!

# model の parameter
model_params = {  
    "n_estimators": 10000,
    "objective": 'rmse',
    "learning_rate": 0.01,
    "num_leaves": 31,
    "random_state": 2021,
    "n_jobs": -1,
    "importance_type": "gain",
    'colsample_bytree': .5,
    "reg_lambda": 5
    }

# fit 時の parameter
fit_params = {
    "early_stopping_rounds": 100,
    "verbose": False
}

次に cross validation 用の関数と評価関数を定義します
今回は、PlaceIDをキーとしたGroup K foldという方法を用います

class GroupKFold:
    """
    GroupKFold with random shuffle with a sklearn-like structure
    """

    def __init__(self, n_splits=4, shuffle=True, random_state=42):
        self.n_splits = n_splits
        self.shuffle = shuffle
        self.random_state = random_state

    def get_n_splits(self, X=None, y=None, group=None):
        return self.n_splits

    def split(self, X=None, y=None, group=None):
        kf = model_selection.KFold(n_splits=self.n_splits, shuffle=self.shuffle, random_state=self.random_state)
        unique_ids = group.unique()
        for tr_group_idx, va_group_idx in kf.split(unique_ids):
            # split group
            tr_group, va_group = unique_ids[tr_group_idx], unique_ids[va_group_idx]
            train_idx = np.where(group.isin(tr_group))[0]
            val_idx = np.where(group.isin(va_group))[0]
            yield train_idx, val_idx


# PlaceID をキーにした Group K fold
def make_gkf(X, y, n_splits=5, random_state=2020):
    gkf = GroupKFold(n_splits=n_splits, random_state=random_state)
    return list(gkf.split(X, y, group))


# 評価関数(log(目的変数)を使っているためrmseを用いる)
def root_mean_squared_error(y_true, y_pred):
    return mean_squared_error(y_true, y_pred) ** .5

学習を始める前に feature importance を見てみます

# GBDTのfeature importanceを出力&可視化
def tree_importance(X, y, model, model_params, fit_params, cv, folds):
    est = model(**model_params)
    feature_importance_df = pd.DataFrame()
    for i, (tr_idx, va_idx) in enumerate(cv(X, y, n_splits=folds)):
        tr_x, va_x = X.values[tr_idx], X.values[va_idx]
        tr_y, va_y = y.values[tr_idx], y.values[va_idx]
        est.fit(tr_x, tr_y,
                eval_set=[[va_x, va_y]],
                **fit_params) 
        _df = pd.DataFrame()
        _df['feature_importance'] = est.feature_importances_
        _df['column'] = X.columns
        _df['fold'] = i + 1
        feature_importance_df = pd.concat([feature_importance_df, _df], axis=0, ignore_index=True)
    order = feature_importance_df.groupby('column') \
                .sum()[['feature_importance']] \
                .sort_values('feature_importance', ascending=False).index[:50]
    fig, ax = plt.subplots(figsize=(12, max(4, len(order) * .2)))
    sns.boxenplot(data=feature_importance_df, y='column', x='feature_importance', order=order, ax=ax,
                  palette='viridis')
    fig.tight_layout()
    ax.grid()
    ax.set_title('feature importance')
    fig.tight_layout()
    plt.show()
    return fig, feature_importance_df
    
fig, _ = tree_importance(train_x, train_y, model, model_params, fit_params, cv=make_gkf, folds=5)

特徴量重要度

特徴量重要度をみるとplace_id_vecs_featuresで作った特徴量やPlaceIDをキーにした集約特徴量がけっこう効いているようです

学習の段階までたどり着けました。学習用の関数を定義してから、そこに上記で設定したモデルとパラメータをぶち込んで学習を回してみましょう
今回は 5fold, 3 seed average です

# cross validation で学習
def predict_cv(X, y, model, model_params, fit_params, cv, folds, seeds, metrics, name=""):
    oof_seeds = []; scores_seeds = []; models = {}
    for seed in seeds:
        oof = []; va_idxes = []; scores = []
        train_x, train_y = X.values, y.values
        fold_idx = cv(train_x, train_y, n_splits=folds, random_state=seed) 
        if "random_state" in model_params:
            model_params["random_state"] = seed
            
        # train and predict by cv folds
        for cv_num, (tr_idx, va_idx) in enumerate(fold_idx):
            tr_x, va_x = train_x[tr_idx], train_x[va_idx]
            tr_y, va_y = train_y[tr_idx], train_y[va_idx]
            va_idxes.append(va_idx)
            est = model(**model_params)

            # fitting - train
            est.fit(tr_x, tr_y,
                    eval_set=[[va_x, va_y]],
                    **fit_params) 
            model_name = f"{name}_SEED{seed}_FOLD{cv_num}_model"
            models[model_name] = est  # save model

            # predict - validation
            pred = est.predict(va_x)
            oof.append(pred)

            # validation score
            score = metrics(va_y, pred)
            scores.append(score)
            print(f"SEED:{seed}, FOLD:{cv_num} >>>> val_score:{score:.4f}")

        # sort as default
        va_idxes = np.concatenate(va_idxes)
        oof = np.concatenate(oof)
        order = np.argsort(va_idxes)
        oof = oof[order]
        oof_seeds.append(oof)
        scores_seeds.append(np.mean(scores))

    oof = np.mean(oof_seeds, axis=0)
    print(f"model:{name} score:{metrics(train_y, oof):.4f}\n")
    return oof, models
    
# oof予測値と学習済みモデルをGet!
oof, models = predict_cv(train_x, train_y, model,  # 説明変数, 目的変数, モデル
                         model_params,  # モデルパラメータ 
                         fit_params,  # fit時のパラメータ
                         cv=make_gkf,  # cv strategy 
                         folds=5,   # num folds
                         seeds=[0, 1, 2], # 3 seed average
                         metrics=root_mean_squared_error,  # 評価関数
                         name="lgbm"
                        )			
SEED:0, FOLD:1 >>>> val_score:0.5399
SEED:0, FOLD:2 >>>> val_score:0.5804
SEED:0, FOLD:3 >>>> val_score:0.5198
SEED:0, FOLD:4 >>>> val_score:0.5376
SEED:1, FOLD:0 >>>> val_score:0.5294
SEED:1, FOLD:1 >>>> val_score:0.5456
SEED:1, FOLD:2 >>>> val_score:0.5112
SEED:1, FOLD:3 >>>> val_score:0.5494
SEED:1, FOLD:4 >>>> val_score:0.5384
SEED:2, FOLD:0 >>>> val_score:0.5494
SEED:2, FOLD:1 >>>> val_score:0.5505
SEED:2, FOLD:2 >>>> val_score:0.5533
SEED:2, FOLD:3 >>>> val_score:0.5280
SEED:2, FOLD:4 >>>> val_score:0.5530
model:lgbm score:0.5353

結果は... CV:0.5353!!でした、どうなんでしょうか
次に学習済みモデルを用いてテストデータを予測してみます

# テストデータを学習済みモデルで予測
def predict_test(test_x, models):
    preds = []
    for name, est in models.items():
        print(f"{name} >>>> ")
        _pred = est.predict(test_x)
        preds.append(_pred)
    preds = np.mean(preds, axis=0)
    return preds
    
# 提出用の予測値をGet!
preds = predict_test(test_x, models)
lgbm_SEED0_FOLD0_model >>>> 
lgbm_SEED0_FOLD1_model >>>> 
lgbm_SEED0_FOLD2_model >>>> 
lgbm_SEED0_FOLD3_model >>>> 
lgbm_SEED0_FOLD4_model >>>> 
lgbm_SEED1_FOLD0_model >>>> 
lgbm_SEED1_FOLD1_model >>>> 
lgbm_SEED1_FOLD2_model >>>> 
lgbm_SEED1_FOLD3_model >>>> 
lgbm_SEED1_FOLD4_model >>>> 
lgbm_SEED2_FOLD0_model >>>> 
lgbm_SEED2_FOLD1_model >>>> 
lgbm_SEED2_FOLD2_model >>>> 
lgbm_SEED2_FOLD3_model >>>> 
lgbm_SEED2_FOLD4_model >>>> 

得られたoof予測値、テストデータ予測値、真の目的変数の値をプロットしてみます

def result_plot(train_y, oof, preds):
    name = "result"
    fig, ax = plt.subplots(figsize=(8, 8))
    sns.distplot(train_y, label='trian_y')
    sns.distplot(oof, label='oof')
    sns.distplot(preds, label='preds')
    ax.legend()
    ax.grid()
    ax.set_title(name)
    fig.tight_layout()
    plt.show()

result_plot(train_y, oof, preds)

result

真の値の分布とは少し違う感じですね

make submission

最後に提出用ファイルを拵えておしまいです。このときテストデータ予測値は対数をとった値になっているので直す必要があります

sample_sub["LandPrice"] = np.expm1(preds)  # RMSE → RMSLE
sample_sub.to_csv('../submission/baseline001.csv', index=False)

投稿!!
LB

summary

  • 以上が今回の記事の内容です! 何か参考になるところなんかあれば嬉しいです。
  • 大きく効いた特徴は、PlaceIDの場所は一体どんな場所なのかを示せるような特徴量す。つまりPlaceID単位で同じ値をとる特徴量ですね。時間変化の重みを場所ごとで同一な値につけて変化を加える拡張や、モデル自体をパネルデータを扱えるもに変えるなどもできそうです。
  • 面積はその土地を表す代表的な指標でもあるので、面積まわりで特徴量を考えてもいいのかもしれないです。
  • 学習用データを、欠損値補完をし完全なパネルデータの形にすることも考えられます。
  • クラスタリングで似たような地域をグルーピングすることで、新しい特徴量を作ることも可能だと思います。
  • CV-LB相関については、自分の提出結果をみると、多少ブレがあるのももの相関はとれていそうです。

参考記事

この記事ではあまり触れていなかったことは以下の記事で詳しい説明があります。
私自身もめちゃくちゃ参考にしました!ありがとうございます!