中古マンション価格予測
田村 圭祐
R言語で機械学習実践
R言語のモダンなデータ分析パッケージである『tidymodels』と『tidyverse』を使って中古マンション価格予想を実践します。
機械学習の分析の流れ
以下がtidyverse、tidymodelsを使った分析の流れです。{}内はパッケージ名
1. データ読み込み {readr}
2. データ理解 {dplyr}{ggplot2}
3. データ分割 {rsample}
4. データ加工 {recipes}
5. 機械学習ルール設定 {parsnip}
6. ワークフロー {workflows}
7. ハイパーパラメータチューニング {tune}{dials}
8. 機械学習 {tune}
9. モデル評価
10. データ予測
11. アンサンブル
1.データ読み込み {readr}
データをread_csv関数を使い読み込みます。ファイル数が9個に分割されているのでbind_rows関数で縦にデータを結合をしていきます。
また、複数のファイル読み込み時に同じカラムでも異なる型が読み込まれる恐れがあるため、col_types引数で型を予め明示します。
#ライブラリの読み込み
library(tidyverse)
library(tidymodels)
library(zipangu)
library(withr)#データの読み込み
df <- tibble()
for (i in 1:9) {
x <- tibble(read_csv(paste("train/0",as.character(i),".csv",sep = ""),col_types = 'dcldccccccclllcccclllcddcccd'))
df <- bind_rows(df,x)
}
rm(x,i)
#確認
df## # A tibble: 45,193 × 28
## ID 種類 地域 市区…¹ 都道府…² 市区…³ 地区名 最寄…⁴ 最寄…⁵ 間取り
## <dbl> <chr> <lgl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1088170 中古マン… NA 1204 北海道 旭川市 ７条通 旭川 16 ２ＤＫ
## 2 1048405 中古マン… NA 1106 北海道 札幌市… 真駒… 真駒内 25 ３Ｌ…
## 3 1030233 中古マン… NA 1104 北海道 札幌市… 菊水… 菊水 6 ４Ｌ…
## 4 1000080 中古マン… NA 1101 北海道 札幌市… 大通西 西１１… 4 ３Ｌ…
## 5 1065071 中古マン… NA 1108 北海道 札幌市… 大谷… ひばり… 9 ３Ｌ…
## 6 1040407 中古マン… NA 1105 北海道 札幌市… 月寒… 福住 8 ３Ｌ…
## 7 1041362 中古マン… NA 1105 北海道 札幌市… 平岸… 澄川 8 ３Ｌ…
## 8 1044353 中古マン… NA 1105 北海道 札幌市… 西岡… 月寒中… 30分?6… ３Ｌ…
## 9 1000091 中古マン… NA 1101 北海道 札幌市… 大通西 西１１… 2 ３Ｌ…
## 10 1002272 中古マン… NA 1101 北海道 札幌市… 北５… さっぽ… 10 １Ｋ
## # … with 45,183 more rows, 18 more variables: `面積（㎡）` <chr>,
## # 土地の形状 <lgl>, 間口 <lgl>, `延床面積（㎡）` <lgl>, 建築年 <chr>,
## # 建物の構造 <chr>, 用途 <chr>, 今後の利用目的 <chr>, `前面道路：方位` <lgl>,
## # `前面道路：種類` <lgl>, `前面道路：幅員（ｍ）` <lgl>, 都市計画 <chr>,
## # `建ぺい率（％）` <dbl>, `容積率（％）` <dbl>, 取引時点 <chr>, 改装 <chr>,
## # 取引の事情等 <chr>, `取引価格（総額）_log` <dbl>, and abbreviated variable
## # names ¹市区町村コード, ²都道府県名, ³市区町村名, ⁴`最寄駅：名称`, …
2.データ理解{dplyr}{ggplot2}
データを操作する前に簡単な前処理を行い、基本統計量やグラフでデータ構造を理解します。
日本語名カラムだと不備が出やすいためカラム名はローマ字に変換します。
同時に説明変数に使えそうなカラムを絞り込みます。
(本当はさらに前処理をしつつcount関数を使って詳細にデータ理解を進めたいですが文章量が莫大になるので割愛してます。)
#カラム名を日本語からローマ字に変換
data <- df %>%
select(KAKAKU = "取引価格（総額）_log" ,TODOHUKEN = "都道府県名",KYORI = "最寄駅：距離（分）",MENSEKI = "面積（㎡）" ,KENTIKUNEN = "建築年", KOUZOU = "建物の構造",YOUSEKI = "容積率（％）",KAISOU = "改装")
次に基本統計量やグラフでデータの分布を確認したり、欠損値や文字列の確認を行います。
#データ構造確認
data## # A tibble: 45,193 × 8
## KAKAKU TODOHUKEN KYORI MENSEKI KENTIKUNEN KOUZOU YOUSEKI KAISOU
## <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
## 1 6.74 北海道 16 45 昭和64年 ＲＣ 300 未改装
## 2 6.88 北海道 25 70 平成4年 ＲＣ 200 改装済
## 3 7.52 北海道 6 105 平成20年 ＲＣ 200 未改装
## 4 7.79 北海道 4 100 平成17年 ＲＣ 600 未改装
## 5 7 北海道 9 75 平成19年 ＲＣ 300 未改装
## 6 7.15 北海道 8 70 平成8年 ＳＲＣ 200 <NA>
## 7 6.65 北海道 8 55 昭和56年 ＲＣ 80 未改装
## 8 7.11 北海道 30分?60分 75 平成3年 ＳＲＣ 200 改装済
## 9 7.84 北海道 2 90 平成24年 ＲＣ 600 <NA>
## 10 6.56 北海道 10 30 昭和56年 ＳＲＣ 400 未改装
## # … with 45,183 more rows#基本統計量算出及び欠損値確認
summary(data)## KAKAKU TODOHUKEN KYORI MENSEKI
## Min. :4.000 Length:45193 Length:45193 Length:45193
## 1st Qu.:6.806 Class :character Class :character Class :character
## Median :7.079 Mode :character Mode :character Mode :character
## Mean :7.012
## 3rd Qu.:7.279
## Max. :8.255
##
## KENTIKUNEN KOUZOU YOUSEKI KAISOU
## Length:45193 Length:45193 Min. : 60.0 Length:45193
## Class :character Class :character 1st Qu.: 200.0 Class :character
## Mode :character Mode :character Median : 200.0 Mode :character
## Mean : 285.8
## 3rd Qu.: 400.0
## Max. :1000.0
## NA's :1920#可視化
ggplot(data = data, mapping = aes(x = KAKAKU, y = YOUSEKI)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)
基本統計量やグラフを見ると『容積』が200～280あたりに集中しているようです。
もっと視覚的に分かりやすいようにヒストグラムで確認してみましょう。
ggplot(data = data, mapping = aes(x = KAKAKU)) +
geom_histogram(colour="white") +
xlim(5, 9)
ggplot(data = data, mapping = aes(x = YOUSEKI)) +
geom_histogram(colour="white",bins = 10)
どうやら『価格』は中央値・平均値ともに7の正規分布になっており、
『容積』の場合は200～280の極端な分布になっているのが分かります。
3.データ分割 {rsample}
ここからは機械学習に向けて『tidymodels』を使っていきます。今までの『tidyverse』と同じ操作感覚でコーディングができ、パイプラインも整頓されている便利なパッケージです。
まずはデータを訓練データとテストデータに分割していきます。
データ全体の内『訓練データ8割：テストデータ2割』で分割します。
#初期分割
initial_split =
initial_split(data, prop = 0.8, strata = "KAISOU") %>% #『改装』の割合に合わせる
with_seed(1234, .) #シード値固定
train = training(initial_split)
test = testing(initial_split)
次に交差検証(Cross Validation)で学習ができるように『stratified分割』をしていきます。
ここでは訓練データ全体の内『学習データ3:検証データ1』で4分割します。
#stratified分割
stratified_splits =
vfold_cv(train, v = 4, strata = "KAISOU") %>%
with_seed(1234, .)
4.データ加工 {recipes}
データ加工(特徴量エンジニアリング)では機械学習がデータを使って学習しやすいように前処理を行います。今回の場合は『Xgboost』というアルゴリズムを使っているため本当は大部分が変換不要で学習が行えますが、{recipes}の簡単なお手本として書き記します。
#{recipes}でデータ加工のレシピを作成
rec <-
recipe(KAKAKU ~ ., data = train) %>%
#現年度より築年数を算出
step_mutate( KENTIKUNEN = str_replace(KENTIKUNEN , "戦前", "昭和20年")
,KENTIKUNEN = convert_jyear(KENTIKUNEN) #和暦→西暦
,KENTIKUNEN = abs(KENTIKUNEN - as.integer(substr(Sys.Date(),1,4)))) %>%
#距離の文字列を数値に変換、型も数値型に変える。
step_mutate( KYORI = str_replace(KYORI , "2H\\?" , "120")
,KYORI = str_replace(KYORI , "1H\\?1H30" , "75")
,KYORI = str_replace(KYORI , "1H30\\?2H" , "105")
,KYORI = str_replace(KYORI , "30分\\?60分", "45")
,KYORI = as.integer(KYORI)) %>%
#面積の文字列を数値に変換、型も数値型に変える。
step_mutate( MENSEKI = str_replace(MENSEKI , "2000㎡以上" , "2000")
,MENSEKI = as.integer(MENSEKI)) %>%
#構造は3つの単語にある程度まとめ、欠損値をＲＣにする。
step_mutate( KOUZOU = str_replace(KOUZOU , "ＲＣ、ブロック造|ＲＣ、木造|ＲＣ、鉄骨造", "ＲＣ")
,KOUZOU = str_replace(KOUZOU ,"ＳＲＣ、ＲＣ|ＳＲＣ、ＲＣ、鉄骨造|ＳＲＣ、鉄骨造", "ＳＲＣ")
,KOUZOU = str_replace(KOUZOU , "軽量鉄骨造", "鉄骨造")
,KOUZOU = replace_na(KOUZOU , "ＲＣ")) %>%
#数値型の欠損値を中央値で埋める
step_impute_median(all_numeric_predictors()) %>%
#chr型をfct型にする
step_string2factor(all_nominal_predictors()) %>%
#改装の欠損値を決定木で埋める
step_impute_bag(
KAISOU,
impute_with = imp_vars(all_predictors())
) %>%
#標準化(イオ・ジョンソン変換)
step_normalize(all_numeric_predictors()) %>%
step_YeoJohnson(all_numeric_predictors()) %>%
#カテゴリ変数の変換(one-hot encoding)
step_dummy(all_nominal_predictors(), one_hot = TRUE) #%>%
#以下の処理は不要。データを確認したい場合に使用する。
#上のパイプラインのコメントアウトも外して使ってください。
# prep() %>% #処理内容の決定
# bake(new_data = train) #処理適用
# rec長文になりましたが以上でデータ加工は完了です。{recipes}には加工するために様々な関数が用意されており、『tidymodels : step_mutate』 = 『tidyverse : mutate』と同じ感覚で関数が使えるのも大きなメリットです。
5. 学習ルール設定 {parsnip}
機械学習がどのようにしてデータを学習していくのかルールを設定します。
今回はメジャーなxgboost(勾配ブースティング木)で学習をします。
ハイパーパラメーターを設定することもできますが、具体的には後で{tune}と{dials}を使って設定します。
#{parsnip}学習ルール設定
boost_tree_xgboost_spec <-
boost_tree(tree_depth = tune(), #ハイパラの設定
trees = tune(),
learn_rate = tune(),
min_n = tune(),
loss_reduction = tune(),
sample_size = tune(),
stop_iter = tune()
) %>%
set_engine('xgboost') %>%
set_mode('regression') #回帰分析
6.ワークフロー {workflows}
{recipes}で作った『データ加工のレシピ』と{parsnip}で作った『学習ルールの設定』を{workflows}を使って『ワークフロー』としてひとまとめにします。
ワークフローにまとめることで、後で機械学習をする時のコーディングの負担が大幅に減り、今後の改修も楽になります。
#{workflows}ワークフロー
wf = workflow() %>%
add_model(boost_tree_xgboost_spec) %>% #学習ルールの設定
add_recipe(rec) #データ加工のレシピ
7.ハイパーパラメータチューニング {tune}{dials}
{tune}と{dials}を使ってハイパーパラメータ(以後ハイパラ)を設定します。
まずは初期設定のままですがハイパラを確認をしましょう。
#ハイパラ確認
wf %>%
parameters() %>%
pull_dials_object("trees")## # Trees (quantitative)
## Range: [1, 2000]wf %>%
parameters() %>%
pull_dials_object("min_n")## Minimal Node Size (quantitative)
## Range: [2, 40]wf %>%
parameters() %>%
pull_dials_object("tree_depth")## Tree Depth (quantitative)
## Range: [1, 15]wf %>%
parameters() %>%
pull_dials_object("learn_rate")## Learning Rate (quantitative)
## Transformer: log-10 [1e-100, Inf]
## Range (transformed scale): [-10, -1]wf %>%
parameters() %>%
pull_dials_object("loss_reduction")## Minimum Loss Reduction (quantitative)
## Transformer: log-10 [1e-100, Inf]
## Range (transformed scale): [-10, 1.5]wf %>%
parameters() %>%
pull_dials_object("sample_size")## Proportion Observations Sampled (quantitative)
## Range: [0.1, 1]wf %>%
parameters() %>%
pull_dials_object("stop_iter")## # Iterations Before Stopping (quantitative)
## Range: [3, 20]
次にこれらハイパラをチューニングをしていきましょう。
探索範囲を更新し、グリッド(探索するときのルール)を作成します。
#探索範囲更新
range <-
wf %>%
parameters() %>%
update(trees = trees(c(5, 2500))) %>%
finalize(rec %>% prep() %>% bake(train) %>% select(-KAKAKU))
#グリッド作成
grid <-
range %>%
grid_latin_hypercube(size = 50) %>% #パラメータ空間全体をカバーする配置
with_seed(1234, .) #シード値固定ハイパラチューニングは以上です。
次はいよいよModeling(機械学習)をしていきます。
8.機械学習 {tune}
ここからようやく{tune}を使った機械学習をしていきます。作成した訓練データやワークフロー、グリッドをセットしましょう。また、損失関数(誤差関数)にはMAEを使用します。Xgboostの交差検証法なので学習にはかなりの時間が掛かります。
tune <-
wf %>% #ワークフロー
tune_grid(
resamples = stratified_splits, #訓練データ
grid = grid, #グリッド
control = control_grid( #出力結果のオプション のちのモデルや予想結果を出す
save_pred = TRUE,
extract = extract_model,
),
metrics = metric_set(mae) #損失関数(誤差関数)
)#確認
tune## # Tuning results
## # 4-fold cross-validation using stratification
## # A tibble: 4 × 6
## splits id .metrics .notes .extracts .predictions
## <list> <chr> <list> <list> <list> <list>
## 1 <split [27115/9039]> Fold1 <tibble [50 × 11]> <tibble> <tibble> <tibble>
## 2 <split [27115/9039]> Fold2 <tibble [50 × 11]> <tibble> <tibble> <tibble>
## 3 <split [27116/9038]> Fold3 <tibble [50 × 11]> <tibble> <tibble> <tibble>
## 4 <split [27116/9038]> Fold4 <tibble [50 × 11]> <tibble> <tibble> <tibble>これで学習は以上です。tibble形式でデータが格納されていることが分かります。
交差検証法のために4分割し、さらにグリッドで探索範囲を50パターンのハイパラで学習を行っております。
9.モデル評価
tuneに格納されているtibble形式のデータを確認していきます。
学習して得た予想結果と損失関数に基づいた評価指標を見ていきましょう。
#特定のFoldに対するassessmentの予測値
tune %>%
pluck(".predictions", 1) #Fold1の予測値
#全Foldの評価指標まとめ
tune %>%
collect_metrics()## # A tibble: 451,950 × 11
## .pred KAKAKU .row trees min_n tree_depth learn_rate loss_r…¹ sampl…² stop_…³
## <dbl> <dbl> <int> <int> <int> <int> <dbl> <dbl> <dbl> <int>
## 1 6.94 7.08 3 428 2 1 0.0206 2.93e-8 0.257 9
## 2 6.80 7.00 5 428 2 1 0.0206 2.93e-8 0.257 9
## 3 6.94 6.96 6 428 2 1 0.0206 2.93e-8 0.257 9
## 4 7.39 7.64 7 428 2 1 0.0206 2.93e-8 0.257 9
## 5 7.07 7.08 12 428 2 1 0.0206 2.93e-8 0.257 9
## 6 7.42 7.45 17 428 2 1 0.0206 2.93e-8 0.257 9
## 7 6.78 7.08 24 428 2 1 0.0206 2.93e-8 0.257 9
## 8 7.39 7.36 25 428 2 1 0.0206 2.93e-8 0.257 9
## 9 7.29 7.43 27 428 2 1 0.0206 2.93e-8 0.257 9
## 10 7.19 6.99 28 428 2 1 0.0206 2.93e-8 0.257 9
## # … with 451,940 more rows, 1 more variable: .config <chr>, and abbreviated
## # variable names ¹loss_reduction, ²sample_size, ³stop_iter## # A tibble: 50 × 13
## mean trees min_n tree_depth learn_…¹ loss_…² sampl…³ stop_…⁴ .metric .esti…⁵
## <dbl> <int> <int> <int> <dbl> <dbl> <dbl> <int> <chr> <chr>
## 1 0.139 428 2 1 2.06e- 2 2.93e-8 0.257 9 mae standa…
## 2 0.195 572 3 12 6.56e- 3 1.53e-7 0.268 17 mae standa…
## 3 6.24 1412 4 6 3.00e- 5 7.70e-8 0.935 16 mae standa…
## 4 6.51 197 5 2 5.08e- 8 3.51e-4 0.668 18 mae standa…
## 5 6.43 1365 5 4 8.83e- 6 9.75e-9 0.384 7 mae standa…
## 6 6.51 263 7 14 4.54e- 9 1.58e-9 0.341 6 mae standa…
## 7 6.51 1806 7 15 7.30e-10 4.48e-2 0.362 19 mae standa…
## 8 6.51 1666 8 5 1.46e- 8 1.16e-5 0.684 15 mae standa…
## 9 6.48 377 8 7 1.28e- 5 1.35e-6 0.606 6 mae standa…
## 10 6.51 1933 9 6 4.17e- 7 2.86e+1 0.504 8 mae standa…
## # … with 40 more rows, 3 more variables: n <int>, std_err <dbl>, .config <chr>,
## # and abbreviated variable names ¹learn_rate, ²loss_reduction, ³sample_size,
## # ⁴stop_iter, ⁵.estimator特定のFoldをハイパラごとに予測値(.pred)と実測値(KAKAKU)を見比べることができます。
次の評価指標まとめではmeanが評価指標となります。今回の評価指標はMAEを使っているので、小さいほど精度が高いです。
※実行環境によって値は僅かに変わります。
それでは評価指標の結果が良かった2つを見てみましょう。
#評価指標が良い結果を2つ抽出
tune %>%
show_best(
metric = "mae",
n = 2
) %>%
select(mean,everything())## # A tibble: 2 × 13
## mean trees min_n tree_depth learn_r…¹ loss_…² sampl…³ stop_…⁴ .metric .esti…⁵
## <dbl> <int> <int> <int> <dbl> <dbl> <dbl> <int> <chr> <chr>
## 1 0.118 929 17 10 0.00941 0.00122 0.985 5 mae standa…
## 2 0.127 1505 29 14 0.0732 0.00425 0.912 8 mae standa…
## # … with 3 more variables: n <int>, std_err <dbl>, .config <chr>, and
## # abbreviated variable names ¹learn_rate, ²loss_reduction, ³sample_size,
## # ⁴stop_iter, ⁵.estimator
次に評価指標の結果が良かったハイパラをモデルにセットします。
# 良さげなハイパラをワークフローにセット
good_hypara_1 =
tune %>%
show_best(
metric = "mae",
n = 2
) %>%
.[1,]
good_hypara_2 =
tune %>%
show_best(
metric = "mae",
n = 2
) %>%
.[2,]
final_wf_1 =
wf %>%
finalize_workflow(good_hypara_1)
final_wf_2 =
wf %>%
finalize_workflow(good_hypara_2)
final_model_1 =
final_wf_1 %>%
fit(train) %>%
with_seed(1234, .)
final_model_2 =
final_wf_2 %>%
fit(train) %>%
with_seed(1234, .)
10.データ予測
いよいよラストスパートです。モデルを使ってデータ予測をしていきます。全体の内の2割をテストデータとして使えるので、テストデータで予測を行います。
#評価データにモデルを適用し，予測値を算出
pred_1 =
predict(final_model_1, new_data = test) %>%
set_names(str_c(names(.), "_1"))
pred_2 =
predict(final_model_2, new_data = test) %>%
set_names(str_c(names(.), "_2"))
pred_1 %>%
bind_cols(pred_2) %>%
mutate(KAKAKU = test$KAKAKU)## # A tibble: 9,039 × 3
## .pred_1 .pred_2 KAKAKU
## <dbl> <dbl> <dbl>
## 1 7.08 7.07 7.15
## 2 6.99 6.99 6.68
## 3 7.32 7.33 7.40
## 4 7.46 7.43 7.34
## 5 6.86 7.03 6.53
## 6 7.17 7.17 7.20
## 7 7.37 7.15 7.18
## 8 6.94 7.01 7.08
## 9 6.78 6.80 7.04
## 10 7.28 7.36 7.28
## # … with 9,029 more rows2つのモデルで予測ができていることが分かります。
11.アンサンブル
最後に複数のモデルを合体(アンサンブル)させてみましょう。アンサンブルをすることで過学習を抑え、より柔軟な予測ができるようになり精度を上げることができます。
ensemble =
bind_cols(pred_1,pred_2) %>%
transmute(
.pred = (.pred_1 + .pred_2) / 2 ,
) %>%
mutate(KAKAKU = test$KAKAKU)
ensemble %>%
mae(truth = KAKAKU, estimate = .pred)## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 mae standard 0.116.estimateの数値が最終的な評価指標となります。アンサンブルをすることで以前より精度が少し高まりました。
それでは最後に機械学習が予測した予測値(.pred)と実測値(KAKAKU)を見比べてみましょう！
ensemble## # A tibble: 9,039 × 2
## .pred KAKAKU
## <dbl> <dbl>
## 1 7.08 7.15
## 2 6.99 6.68
## 3 7.32 7.40
## 4 7.44 7.34
## 5 6.94 6.53
## 6 7.17 7.20
## 7 7.26 7.18
## 8 6.97 7.08
## 9 6.79 7.04
## 10 7.32 7.28
## # … with 9,029 more rowsEnjoy!