Kaggle - House Prices挑戦記①：全ての作業を言語化する

0 はじめに

0.1 この記事の目的

今回はKaggleの有名なコンペティション、House Pricesに挑戦します。

このコンペ、あまりにも有名すぎて調べるといくらでも解説などが載っています。しかし、多くの解説が「これが上手くいきます」以上の説明がなく、その手前の「なぜその作業をしようと考えたのか」「どんなことを試してみたら上手く行かなかったのか」や「その失敗からなぜそのような作業をしたのか」という説明が少ないと思います。

なので、この記事では「データをのぞいてみる」ところから、可能な限り全ての作業を一つずつ「なぜそうしたのか」の説明もしながらHouse Pricesコンペティションを進めていきたいと思います。

ですから、この記事の目的は「Kaggleに勝つ」「ハイスコアを目指す」ではなく、「MLエンジニアとしての思考・試行を分解する」ことにあります。他人の記事を探せばいくらでもハイスコアを取る方法は書いてありますが、一つひとつの作業や試行を真の意味で理解をするために、それらはいっさい参考にしていません。

言うなればこの記事は「コンペの挑戦実況」を書き起こしたものです。ですから、とても冗長な記事なる上に、つど話数が分かれていくと思います。悪しからず。

0.2 筆者について

筆者は2025年4月1日から機械学習エンジニアにキャリアチェンジをした、完全なる初心者です。前職も全く関係のないことをしていました。ですから、この記事で書き起こしたコードや考え方には間違いや瑕疵、非効率的な部分が多分にあるかと思います。また、必ずしも業界的なベストプラクティスに沿ったものではないです。ご了承ください。

0.3 環境

OS: Linux (WSL - Ubuntu)
Jupyter Notebook (VSCode)
Python 3.12.3

1 セットアップ

まじでここから言語化していきます。
まずはプロジェクト用のディレクトリを作成し、移動。

mkdir -p kaggle-house-price/{data/{raw,processed,image},notebooks/{utils,constants},models}

cd kaggle-house-price

大体この階層にしています。dataは無加工と処理済み、それからmatplotlibなどで視覚化したデータで分けられるようにしています。

notebooksの中には、使いまわせる関数や定数を入れておくutilsモジュールとconstantsモジュールを作っています。constantsは主に保存、読み込み用のパスを保存したりしています。

また、学習済みモデルを保存したりそこから呼び出すためのmodelsも作っておきます。

次に、仮想環境を作成して仮想環境内に入ります。

python -m venv .venv

source .venv/bin/acitvate

ほぼ確定で使うパッケージをインストール。

pip install numpy pandas scikit-learn jupyter matplotlib seaborn

今回は、最近勉強した勾配ブースティングを使ってみようかなーと思ったので、XGBoostも追加しました。

VSCodeで開きます。

code .

準備は整いました。

2 EDA（Exploratory Data Analysis：探索的データ分析）

2.1 ノートブック設定

まずは探索用のノートブックを作成。

touch notebooks/{0000_eda.ipynb,0001_process.ipynb,0100_baseline.ipynb}

探索、前処理、ベースラインモデル（最初に作る、基準となるモデル）、各モデルで分けて作ることが多いです。また、探索と前処理を同時に開いておき、探索中に役立ったデータ整形を即時に前処理に反映させたりしています。
また、各ノートブックには自分なりのわかりやすい４桁コードをつけています。上2桁は用途・モデルごとに分け、下二桁を順番にしています。

上2桁	用途・モデル
00	データ探索・処理など
01	ベースライン
02	［以下、モデルごとにナンバリング］
$\vdots$	$\vdots$

0000_eda.ipynbを開いたら、まずは必要なモジュールをインポート。

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

plt.style.use('ggplot')
sns.set_palette('Set2')

下の2行で図表の色を指定していますが、これは好みの問題なので別になくてもいいです。

2.2 データの概観を確認

Kaggleのコンペティションページから取得したデータを開き、とりあえず中身を確認します。

train = pd.read_csv('../data/raw/train.csv', index_col=0)
train.head()

# 別のセルで
train.info()

# OUTPUT

#   Column         Non-Null Count  Dtype  
---  ------         --------------  -----  
 0   MSSubClass     1460 non-null   int64  
 1   MSZoning       1460 non-null   object 
 2   LotFrontage    1201 non-null   float64
 3   LotArea        1460 non-null   int64  
...
 78  SaleCondition  1460 non-null   object 
 79  SalePrice      1460 non-null   int64  
dtypes: float64(3), int64(34), object(43)
memory usage: 923.9+ KB

80の変数と1460のサンプルがあるそうです。ひえー。
しかも、欠損値が結構ありそう。また、データの概要を見ると、順序尺度・名義尺度も多くありそうです。

順序尺度や名義尺度は、訓練データにはそのバリエーションが全て揃っているけれどテストデータには一部しかない、といった場合もあり、特徴量エンジニアリングにおいて、例えばone-hotエンコードをする際に訓練とテストでカラムの数が異なる、なんてことにもなりかねません。

訓練データでfitしたモデルは、同じ特徴量の構成をしているデータでないと学習できなくて面倒な思いをすることになるので、特徴量を作るときはtrainとtestを合体させて一つのデータセットとして作っていこうと思います。

次に、目的変数である家の販売価格の分布表を見てみます。

sns.histplot(train['SalePrice'], kde=True)

• 出力
  
  全体的にだいぶ右に裾が長いですね。多くのモデルが、正規分布に従うデータの方が扱いやすいので、いったん対数変換してみます。

train['log_price'] = np.log(train['SalePrice'])
plt.figure(figsize=(12,8))
sns.histplot(train['SalePrice_log'], kde=True)
plt.show()

• 出力
  
  だいぶ正規分布っぽくなりました。歪度と尖度を調べてみます。

train['log_price'].skew(), train['log_price'].kurtosis()

# 出力
# (np.float64(0.12133506220520406), np.float64(0.8095319958036296))

だいぶいい感じなので、この値を予測した上で、提出時に指数変換しようと思います。

2.3 思いつく限り特徴量を作ってみる

次に説明変数を見ていきたいのですが、このコンペのデータセットは特徴量が80個くらいあってかなり面倒くさそうです。自分が家を選ぶときの条件を参考にいくつか見当をつけた上で、目的変数との関係を見てみようと思います。

コンペの特徴量の中で特に重要そうなのは、

• 立地：Neighborhood
• 家の状態：OverallQual, OverallCond, ExterQual, ExterCond, KitchenQual, KithcenCond
• 広さ：LotArea, 1stFlrSF, 2ndFlrSF, GrLivArea, TotalBsmtSF
• 古さ：YearBuilt, YearRemodAdd

これらと対数価格の関係を見てみたいのだけど、Neighborhoodが文字列データかつアメリカのどこかの地名なので土地勘もない。
そこで、立地ごとに価格の平均を求めて、序列をつけてみる。

立地の特徴量

neighbor_price = train.groupby('Neighborhood')['SalePrice'].mean().sort_values(ascending=False)

neighbor_price

• 出力

Neighborhood
NoRidge    335295.317073
NridgHt    316270.623377
StoneBr    310499.000000
Timber     242247.447368
.
.
.

この順番をこのままナンバリングして順序尺度化してもいいのだけど、NoRidgeとNridgHtの平均値の差がStoneBrとTimberの差と同じとは限らないので、min-max法で係数化してみる。
（正直、ターゲットエンコーディングはリークが気になるところではありますが・・・）

from sklearn.preprocessing import MinMaxScaler

mm_scaler = MinMaxScaler()

neighbor_price_mmsclaled = mm_scaler.fit_transform(neighbor_price.values.reshape(-1, 1))
neighbor_price_dict = dict(zip(neighbor_price.index, neighbor_price_mmsclaled))

neighbor_price_dict

• 出力

'NoRidge': array([1.]),
 'NridgHt': array([0.91963169]),
 'StoneBr': array([0.89524992]),
 'Timber': array([0.60692665]),
 .
 .
 .
 'BrkSide': array([0.11092307]),
 'BrDale': array([0.02499708]),
 'IDOTRR': array([0.0065365]),
 'MeadowV': array([0.])}

一応これをカラムとして保存。

train['neighbor_price_scaled'] = train['Neighborhood'].map(neighbor_price_dict).astype(float)

min-maxスケーリングは住宅価格のような極端な外れ値があるデータに弱い（外れ値の影響を受けやすい）。
ただし、「超高級住宅」それ自体にも情報としての価値はあると考えたため、外れ値はそのまま扱いたい。
そこで、これとは別に、対数変換した価格の平均を立地ごとに集計し、こちらは標準化により係数化して、カラムとして保存。

from sklearn.preprocessing import StandardScaler

std_scaler = StandardScaler()

neighbor_log_price = train.groupby('Neighborhood')['LogPrice'].mean().sort_values(ascending=False)

neighbor_log_price_stdsclaled = std_scaler.fit_transform(neighbor_log_price.values.reshape(-1, 1))
neighbor_log_price_dict = dict(zip(neighbor_log_price.index, neighbor_log_price_stdsclaled))
train['neighbor_weight'] = train['Neighborhood'].map(neighbor_log_price_dict).astype(float)

中身はこんな感じ。

{'NoRidge': array([1.89799433]),
 'NridgHt': array([1.73171346]),
 'StoneBr': array([1.63202722]),
 .
 .
 .
 'BrDale': array([-1.41692402]),
 'MeadowV': array([-1.63243298]),
 'IDOTRR': array([-1.71363052])}

一応、目的変数との関係を可視化してみます。

sns.boxplot(x=train['neighbor_weight'] , y=train['log_price'])
plt.xticks(range(len(neighbor_log_price)), neighbor_log_price.index, rotation=90)
plt.show()

また、こうして特徴量を追加したときは、data/内にadded_deatures.mdなどを作って追加した特徴量とその内容をメモしておき、後かからどの特徴量がいつどんな処理をしたものかがすぐわかるようにしています。

# 0000_EDA

    - log_price: 'SalePrice'を自然対数変換したもの。

    - neighbor_price_scaled: 'Neighborhood'ごとに'SalePrice'の平均を集計し、min-max法により立地を係数化したもの。

    - neighbor_weight 'Neighborhood'ごとに'log_price'の平均を集計し、標準化により立地を係数化したもの。

この処理がかなり使えそうだったので、utils.pyに関数化して保存しようと思います。

# notbooks/utils/utils.py

def ordinal_to_weight(df: pd.DataFrame, col_name: str, target: str='SalePrice', log_transform=False, scaler_type: str='minmax'):
    """
    順序尺度のカテゴリをターゲットエンコーディングする関数

    Args:
        df (DataFrame): 対象データフレーム
        col_name (str): 順序尺度の列名
        target (str, optional): Defaults to 'SalePrice'.
        log_transform (bool, optional): 対数変換するかどうか. Defaults to False.
        scale_typy (str): スケーリング方法（'standard', 'minmax'）. Defaults to 'minmax'.
        
    Returns:
        DataFrame: 新しい特徴量が追加されたデータフレーム
        dict: 変換用の辞書（テストデータ用）
    """
    
    df_copy = df.copy()
    
    # 対数変換の処理
    if log_transform and target == 'SalePrice':
        target_col = f'log_{target}'
        if target_col not in df_copy.columns:
            df_copy[target_col] = np.log(df_copy[target])
    else:
        target_col = target
        
    # 各カテゴリの平均値を計算
    category_weight = df_copy.groupby(col_name)[target_col].mean()
    
    # スケーリング
    if scaler_type == 'standard':
        from sklearn.preprocessing import StandardScaler
        scaler = StandardScaler()
    elif scaler_type == 'minmax':
        from sklearn.preprocessing import MinMaxScaler
        scaler = MinMaxScaler()
    else:
        raise ValueError('scale_type must be "standard" or "minmax"')
    
    if scaler_type != 'none':
        scaled_values = scaler.fit_transform(category_weight.values.reshape(-1, 1)).flatten()
        category_weight_dict = dict(zip(category_weight.index, scaled_values))
        
    if log_transform and target == 'SalePrice':
        new_col_name = f'log_{col_name}_weight'
    else:
        new_col_name = f'{col_name}_weight'
        
    df_copy[new_col_name] = df_copy[col_name].map(category_weight_dict)
    
    return df_copy, category_weight_dict

この関数を使って、どんどん順序尺度を係数化していきます。

順序・名義尺度を係数化

他に、（勘で）ある程度価格と相関が強そうな順序尺度は、

• MSSubClass: 住居のタイプ（数字なので、文字列に直す必要あり）
• MSZoning: 住居の立地のタイプ（居住用地なのか、など）
• OverallQual: 住宅の材質と仕上げの質（1-10）
• OverallCond: 住宅の全体的な状態（1-10）
• ExterQual: 外装の材質の質（Ex, Gd, TA, Fa, Po）
• ExterCond: 外装の現在の状態（Ex, Gd, TA, Fa, Po）
• BsmtQual: 地下室の高さの評価（Ex, Gd, TA, Fa, Po, NA）
• BsmtCond: 地下室の一般的な状態（Ex, Gd, TA, Fa, Po, NA）
• KitchenQual: キッチンの品質（Ex, Gd, TA, Fa, Po）
• FireplaceQu: 暖炉の品質（Ex, Gd, TA, Fa, Po, NA）
• GarageQual: ガレージの品質（Ex, Gd, TA, Fa, Po, NA）
• GarageCond: ガレージの状態（Ex, Gd, TA, Fa, Po, NA）
• PoolQC: プールの品質（Ex, Gd, TA, Fa, NA）
• Functional: 住宅の機能性（Typ, Min1, Min2, Mod, Maj1, Maj2, Sev, Sal）
• Utilities: 住宅の電気・ガス・上下水道（AllPub, NoSewr, NoSeWa, ELO）

面倒なので、ループで全部係数化します。意味がないものもあるかもしれませんが、特徴量を作るだけ作ってからモデルに突っ込んだり主成分分析で次元を削減したり特徴量重要度を見て取捨選択すればいいやの精神です。

NaN値があると上手く係数化できないので、データを確認して欠損値を適宜補填してから、係数を作っていきます。
まずは、欠損値があるカラムを確認。

ordinal_cols = ['OverallQual', 'OverallCond', 'ExterQual', 'ExterCond', 
                'BsmtQual', 'BsmtCond','KitchenQual', 'FireplaceQu',
                'GarageQual', 'GarageCond', 'PoolQC', 'Functional']

train[ordinal_cols].isna().sum()

• 出力

MSSubClass 0
MSZoning 0
OverallQual 0
OverallCond 0
ExterQual 0
ExterCond 0
BsmtQual 37
BsmtCond 37
KitchenQual 0
FireplaceQu 690
GarageQual 81
GarageCond 81
PoolQC 1453
Functional 0
Utilities 0
dtype: int64

地下室、暖炉、ガレージ、プールに欠損値があるようです。これらの欠損値は、そもそもそれらの設備がないのでは？と考え、それぞれの設備の有無を確認します。

has_bsmt = (train['TotalBsmtSF'] > 0).astype(int)

has_fireplace = (train['Fireplaces'] > 0).astype(int)

has_garage = (train['GarageArea'] > 0).astype(int)

has_pool = (train['PoolArea'] > 0).astype(int)

has_bsmt.value_counts(), has_fireplace.value_counts(), has_garage.value_counts(), has_pool.value_counts()

• 出力
  

思ったとおり、欠損値はそもそもその設備がないことを意味していました。欠損値は「設備なし」という情報を持つので、全て'NONE'文字列で変換していきます。

# 欠損値を'NA'に変換
for col in ordinal_cols:
    if col in train.columns:
        train[col] = train[col].fillna('NONE')

欠損値を埋めたところで、係数化していきます。

weight_dicts = {}

from utils import ordinal_to_weight

# 対数価格についての標準化係数を取得
for col in ordinal_cols:
    if col in train.columns:
        train, weight_dict = ordinal_to_weight(
            df=train,
            col_name=col,
            target='log_price',
            scaler_type='standard'
            )
        weight_dicts[col] = weight_dict
        print(f"{col} を重み付けしました。")
        
# 確認
train.head()

ここまでの作業で、

• 立地：Neighborhood
• 家・設備の状態

ここまで処理が終わりました。次は、以下を考えていきます。

• 広さ：LotArea, 1stFlrSF, 2ndFlrSF, GrLivAreaなど
• 古さ：YearBuilt, YearRemodAdd

広さの特徴量作成

広さについては、敷地面積（LotArea)、階ごとの面積（1stFlrSF, 2ndFlrSF, TotalBsmtSF）、居住可能面積と、かなりかぶっている感があります。上手く圧縮できそうです。
まずは、居住面積に関わる変数同士の相関行列（ヒートマップ）と散布図を見ていきます。

area_features = [
    '1stFlrSF', '2ndFlrSF', 'TotalBsmtSF',
    'GrLivArea', 'LotArea', 'GarageArea'
]

# 相関行列のヒートマップ
plt.figure(figsize=(10, 8))
area_corr = train[area_features].corr()
sns.heatmap(area_corr, annot=True, cmap='coolwarm')
plt.title('Correlation Matrix: Area Features')
plt.show()

# 相互散布図
sns.pairplot(train[area_features], diag_kind='kde')
plt.show()

こんな感じでした。

ヒートマップとペアプロット図からは、次のようなことが読み取れそうです。

• 1階と地下室の強い相関：地下室は1階の真下にあることが多いから、当然か。
• 1階と2階は逆相関：2階建ての家は各階が小さめで平屋は1階が広めっていう傾向がありそう。
• 総居住面積への寄与度：総居住面積（GrLivArea）に対して、２階の面積の相関（0.69）の方が１階のそれ（0.57）より高いのは面白い。2階の有無が総面積の変動に大きく効いてそう。
• 総敷地面積（LotArea）は他とあまり相関しない：敷地面積は建物面積とあまり強くは相関していない（0.3前後）。高級住宅地では敷地が広くても家は小さかったり、逆に敷地が狭くても縦に大きく建てたりがあるから当然か。

これらを踏まえて、床面積因子のような特徴量を作成。また、２階の有無、地下室の有無に対してフラグを設定する。

# 床面積因子を作成
train['total_sf'] = train['1stFlrSF'] + train['2ndFlrSF'].fillna(0) + train['TotalBsmtSF'].fillna(0)

# ２階、地下室、ガレージの有無のフラグを設定
train['has_2nd'] = (train['2ndFlrSF'] > 0).astype(int)

train['has_bsmt'] = (train['TotalBsmtSF'] > 0).astype(int)

train['has_garage'] = (train['GarageArea'] > 0).astype(int)

また、床面積因子に関しても、かなり左に寄った分布なので、対数変換してみる。

同様に、総敷地面積も左によっていたので対数化。

sns.distplot(train['LotArea'])
plt.show() # ものすごく左寄り

train['log_lot_area'] = np.log(train['LotArea'])

他に重要そうな情報としては、ガレージの許容量（車を何台止められるか）、家からでしょうか。こちらはそのまま使用できそうです。

フラグに関して、設備のフラグもここでまとめて作っておきます。

# プールの有無
train['has_pool'] = (train['PoolArea'] > 0).astype(int)

# 全館空調の有無
train['has_central_ac'] = (train['CentralAir'] == 'Y').astype(int)

# 暖炉の有無
train['has_fireplace'] = (train['Fireplaces'] > 0).astype(int)

古さの特徴量作成

次に、建物の古さの特徴量を作っていきます。

まずは、販売時の築年数と、販売時のリフォームからの経過年の特徴量を作成し、目的変数との関係を見てみます。

# 販売年 - 建築年
train['house_age'] = train['YrSold'] - train['YearBuilt']

# 販売年 - 改築年
train['remod_age'] = train['YrSold'] - train['YearRemodAdd']

# 可視化
sns.scatterplot(x='remod_age', y='log_price', data=train, hue='house_age', alpha=0.8)
sns.regplot(x='remod_age', y='log_price', data=train, scatter=False, ci=90, line_kws={'color': 'green', 'lw': 1})
plt.show()

最終改築年が古ければ古いほど、価格が下がるのは思った通りで、また最終改築年が同じでも、築年数が古いものは価格が安くなる傾向にありそうなのも、直感的に理解できます。

次に、築年代ごとのカテゴリを作成しようと思いますが、カテゴリの範囲を調べるために度数分布表を作ります。
度数分布表を作成する関数をutils内に作り（参考：「Pythonで度数分布表を一発で自動生成する」 @TakuTaku36）、年代ごとの度数分布を確認します。

# utils.py
def freq_dist(data, class_width=None):
    """ 度数分布表を作成する関数 """
    data = np.asarray(data)
    if class_width is None:
        class_size = int(np.log2(data.size).round()) + 1
        class_width = round((data.max() - data.min()) / class_size)

    bins = np.arange(0, data.max()+class_width+1, class_width)
    hist = np.histogram(data, bins)[0]
    cumsum = hist.cumsum()

    return pd.DataFrame({'階級値': (bins[1:] + bins[:-1]) / 2,
                         '度数': hist,
                         '累積度数': cumsum,
                         '相対度数': hist / cumsum[-1],
                         '累積相対度数': cumsum / cumsum[-1]},
                        index=pd.Index([f'{bins[i]}以上{bins[i+1]}未満'
                                        for i in range(hist.size)],
                                       name='階級'))

# notebook内
year_freq = freq_dist(train['YearBuilt'], class_width=10)

year_freq.iloc[187:] # 1870年以降しかデータがないため

• 出力
  

大体1900年以前、1900~1940年、以降20年区切り、2000年以降をまとめればよさそう。

train['built_era'] = pd.cut(
    train['YearBuilt'],
    bins=[1800, 1900, 1940, 1960, 1980, 2000, 2020],
    labels=['pre_1900', '1900-1940', '40-60', '60-80', '80-2000', 'post-2000']
)

# log_priceとの関係を可視化
sns.boxplot(x='built_era', y='log_price', data=train)
plt.show()

1900年以前の古い建物が1900~40年代の建物よりも高いのは、感覚としては「保存価値のある古くて良い建物」が多いからだろうか。

住宅の古さと価格の関係も重要な情報なので、年代別の平均価格を作る。平滑化するために、ローリング平均をここでは使用する。

# 建築年ごとの平均価格
year_price = train.groupby('YearBuilt')['SalePrice'].mean()

# ローリング平均
year_price = year_price.rolling(window=5, center=True).mean()

# 欠損値を前後の値で補完
year_price = year_price.interpolate()

# DataFrameに追加
train['year_price'] = train['YearBuilt'].map(year_price)

# 可視化
plt.figure(figsize=(16,8))
sns.barplot(x='YearBuilt', y='year_price', data=train)
plt.xticks(rotation=90, fontsize=8)
plt.show()

どの作業の結果のデータなのかわかるように接頭辞をつけて保存します。

train.to_csv('./data/raw/0000_train.csv')

まとめ

ここまでで、価格に特に影響を与えそうな以下の特徴量を作成しました。

• 立地
• 家の状態
• 広さ
• 古さ

次回は、これらの特徴量をもとにベースラインモデルを作り、貢献度の評価を見てさらに特徴量を追加または削除していこうと思います。

追加した特徴量一覧

    - log_price: 'SalePrice'を自然対数変換したもの。

    - neighbor_price_scaled: 'Neighborhood'ごとに'SalePrice'の平均を集計し、min-max法により立地を係数化したもの。

    - neighbor_weight: 'Neighborhood'ごとに'log_price'の平均を集計し、標準化により立地を係数化したもの。

    - total_sf = 1stFlrSF + 2ndFlrSF (nan = 0) + TotalBsmtSF (nan = 0)

    - log_total_sf: total_sfを対数変換したもの

    - log_lot_area: LotAreaを対数化したもの

    - 下記のそれぞれを重みづけ: 各カテゴリごとにlog_priceの平均を集計。接尾辞に_weight
                'MSSubClass', 'MSZoning',
                'OverallQual', 'OverallCond', 'ExterQual', 'ExterCond', 
                'BsmtQual', 'BsmtCond', 'KitchenQual', 'Utilities,
                'GarageQual', 'GarageCond', 'PoolQC', 'Functional'

    - has_2nd: ２階の有無フラグ

    - has_bsmt: 地下室の有無フラグ

    - has_garage: ガレージの有無フラグ

    - has_pool: プールの有無フラグ
    
    - has_central_ac: 全館空調有無のフラグ

    - house_age: 販売年 - 建築年

    - remod_age: 販売年 - 改築年

    - built_era: 1900年以前、1900~1940年、以降20年区切り、2000年以降で建築年代をカテゴライズ

    - year_price:　建築年ごとの住宅価格のローリング平均

    # そのまま使う
    - GarageCars: ガレージの車両収容台数