“ばらつき”を説明するためのデータ分析

今日の目標

メルボルンの住宅価格データを題材にして、データがもつ“ばらつき”という特徴について理解をする。また、データの可視化と機械学習のプログラミングを通して、ばらつきをどのように“説明する”のかを実際にみてみる。

# 必要なライブラリを読み込む
library(tidyverse) # データ集計・可視化に必要なツールが一通り含まれているライブラリ
library(rpart) # 決定木という機械学習アルゴリズムを使うためのライブラリ
library(rpart.plot) # 決定木の結果を可視化するためのライブラリ
library(randomForest) # ランダムフォレストという機械学習アルゴリズムを使うためのライブラリ

theme_set(theme_gray(base_family="HiraKakuProN-W3")) # 日本語フォントを使用する設定
options(scipen=100) # グラフ軸ラベルの数字が指数表示されるのを回避する設定

Part Ⅰ: データの集計・可視化

データの読み込みと基礎俯瞰

csvの読み込みにはread_csvを使用する。データには文字列型(character)や整数型(integer)などのデータ型がある。読み込み時に型が自動で設定された旨のメッセージが表示されるがcol_types = cols()で非表示にできる。

raw_df <- read_csv("../../data/house_prices.csv", col_types = cols())

今回のデータには次のような変数が含まれている。このうちPrice (住宅価格) に注目して分析をする。

変数	意味	備考
Suburb	地区
Address	住所
Rooms	部屋数
Type	住宅タイプ	br - bedroom(s); h - house,cottage,villa, semi,terrace; u - unit, duplex; t - townhouse;
Price	住宅価格
Method	売却方式	S - property sold; SP - property sold prior; PI - property passed in; PN - sold prior not disclosed; SN - sold not disclosed; VB - vendor bid; W - withdrawn prior to auction; SA - sold after auction; SS - sold after auction price not disclosed.
SellerG	売主
Date	売却年月日
Postcode	郵便番号
Regionname	地域	北部・東部・南西部など
Propertycount	地区内の物件数
Distance	市街中心部からの距離
CouncilArea	選挙区(？)

house_prices.csv

glimpseという関数でデータフレームの内容（レコード数、変数、データ型など）を確認できる。

glimpse(raw_df)

## Observations: 48,433
## Variables: 13
## $ Suburb        <chr> "Melbourne", "Melbourne", "Melbourne", "Melbourn...
## $ Address       <chr> "66/538 Little Lonsdale St", "814/70 Queens Rd",...
## $ Rooms         <int> 2, 2, 1, 2, 3, 2, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, ...
## $ Type          <chr> "u", "u", "u", "u", "u", "u", "u", "u", "u", "u"...
## $ Price         <int> 605000, 649000, 440000, 1320000, 753000, 452000,...
## $ Method        <chr> "S", "PI", "SP", "S", "S", "S", "SP", "S", "SP",...
## $ SellerG       <chr> "McGrath", "hockingstuart", "Harcourts", "Finder...
## $ Date          <chr> "1/04/2017", "1/04/2017", "1/04/2017", "1/09/201...
## $ Postcode      <int> 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, ...
## $ Regionname    <chr> "Northern Metropolitan", "Northern Metropolitan"...
## $ Propertycount <int> 17496, 17496, 17496, 17496, 17496, 17496, 17496,...
## $ Distance      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ CouncilArea   <chr> "Melbourne City Council", "Melbourne City Counci...

summaryという関数でデータフレームの要約統計量(平均値、中央値、最大値、最小値など)を確認できる。

summary(raw_df)

##     Suburb            Address              Rooms            Type          
##  Length:48433       Length:48433       Min.   : 1.000   Length:48433      
##  Class :character   Class :character   1st Qu.: 2.000   Class :character  
##  Mode  :character   Mode  :character   Median : 3.000   Mode  :character  
##                                        Mean   : 3.072                     
##                                        3rd Qu.: 4.000                     
##                                        Max.   :31.000                     
##      Price             Method            SellerG         
##  Min.   :   85000   Length:48433       Length:48433      
##  1st Qu.:  620000   Class :character   Class :character  
##  Median :  830000   Mode  :character   Mode  :character  
##  Mean   :  997898                                        
##  3rd Qu.: 1220000                                        
##  Max.   :11200000                                        
##      Date              Postcode     Regionname        Propertycount  
##  Length:48433       Min.   :3000   Length:48433       Min.   :   39  
##  Class :character   1st Qu.:3051   Class :character   1st Qu.: 4280  
##  Mode  :character   Median :3103   Mode  :character   Median : 6567  
##                     Mean   :3123                      Mean   : 7566  
##                     3rd Qu.:3163                      3rd Qu.:10412  
##                     Max.   :3980                      Max.   :21650  
##     Distance    CouncilArea       
##  Min.   : 0.0   Length:48433      
##  1st Qu.: 7.0   Class :character  
##  Median :11.7   Mode  :character  
##  Mean   :12.7                     
##  3rd Qu.:16.7                     
##  Max.   :55.8

データの前処理

変数Priceは桁が大きくてみづらいため、百万ドル単位に直す。変数Dateは日付などの時間を表すPOSIXctという型に変換する。また、後に機械学習モデルの入力として使うため、character型の変数はfactor型に変換する。factor型はカテゴリ変数を数値計算で扱うためのR特有のデータ型。

df <- raw_df %>% 
  mutate(
    Price = Price / 1000000,
    Date = as.POSIXct(Date, format="%d/%m/%Y"),
    Suburb = as.factor(Suburb),
    Type = as.factor(Type),
    Method = as.factor(Method),
    SellerG = as.factor(SellerG),
    Regionname = as.factor(Regionname),
    CouncilArea = as.factor(CouncilArea)
    )

単変数の可視化

まず、それぞれの変数がどのように分布しているかを確認する。カテゴリ変数であれば棒グラフ、連続変数であればヒストグラムで可視化する場合が多い。分析の初期段階で興味のある変数について網羅的に分布を可視化しておくことが効率的に分析するコツ。後続の分析時にもこの図表を都度確認しながら進めることになる。

変数Method (売却方式)の分布

“S”(通常の売却?)が最も多く、一方“SA”(競売での売却?)はほとんど含まれていない。

df %>% 
  group_by(Method) %>% 
  summarise(n=n()) %>% 
  ggplot(aes(x=Method, y=n, fill=Method)) +
  geom_bar(stat="identity") +
  guides(fill=FALSE)

Regionname (地域)の分布

“Metropolitan”地域が多く、それ以外の“Victoria”州のデータは少ない。

df %>% 
  group_by(Regionname) %>% 
  summarise(n=n()) %>% 
  ggplot(aes(x=Regionname, y=n, fill=Regionname)) +
  geom_bar(stat="identity") +
  coord_flip() + # x軸とy軸を入れ替えて横棒グラフにする
  guides(fill=FALSE)

Distance (市街中心部からの距離)の分布

5km〜20kmくらいがボリュームゾーンで、最大50kmくらいまで分布している。

df %>% 
  ggplot(aes(x=Distance)) +
  geom_histogram(binwidth = 2, fill="orange", color="grey20") +
  coord_cartesian(xlim=c(0, 50))

Price (住宅価格)の分布

住宅価格はやや右に裾の長い滑らかな単峰の分布になっている。

df %>% 
  ggplot(aes(x=Price)) +
  geom_histogram(binwidth = 0.1, fill="blue", color="black", alpha=0.2) +
  geom_vline(xintercept = mean(df$Price), color="red") + # 平均値
  geom_vline(xintercept = median(df$Price), color="blue") + # 中央値
  coord_cartesian(xlim=c(0, 5))

多変数の可視化

2つ以上の変数を使ってデータを可視化する方法はたくさんあるが、ここでは代表的な可視化方法である箱ひげ図と散布図を描画してみる。

Type（住宅タイプ）とPrice (住宅価格)の箱ひげ図

箱ひげ図の箱の下端は価格が小さい順に25%の価格、上端は75%の価格を表す。中央の太線は中央値(50%の価格)。ひげの長さはデータの最大値と最小値を表す。ただし、箱の上(下)端からみて箱の長さの1.5倍よりも大きい(小さい)データ点は外れ値として点で表される。

df %>% 
  ggplot(aes(x=Type, y=Price, color=Type)) +
  geom_boxplot() +
  coord_cartesian(ylim=c(0,3)) +
  guides(color=FALSE)

Distance（市街中心部からの距離）とPrice (住宅価格)の散布図

DistanceとPriceの相関関係をみるために散布図を描画する。Distanceが大きくなると住宅価格が安くなる傾向がみられるが、それほど顕著ではない。

df %>% 
  ggplot(aes(x=Distance, y=Price)) +
  geom_point(color="orange", alpha=0.2) +
  coord_cartesian(xlim=c(0,50), ylim=c(0,5))

ばらつきを表す統計量

上の散布図やヒストグラムをみるとPriceの値は「ばらつき」があることがわかる。このばらつきの大きさを表す指標として次のようなものがある。

偏差平方和(SSD: Sum of Squared Deviation)：平均値と各データとの差を二乗して全てのデータについて足し合わせたもの。 \[ SSD = \sum_{i=1}^{N}(y_i-\bar{y})^2 \]

分散：偏差平方和をデータ数\(N\)で割ったもの。 \[ \sigma^2 = \frac{1}{N}\sum_{i=1}^{N}(y_i-\bar{y})^2 \]

標準偏差：分散の平方根をとったもの。これは平均値と各データが平均してどのくらい離れているかを表す統計量。 \[ \sigma = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_i-\bar{y})^2} \]

Priceの統計量を計算する

Priceの標準偏差は$0.593Mとなっているので、適当に選んだ住宅の価格は平均値$0.997Mと大体$0.593M程度異なっていると期待される。

mean(df$Price) # 平均値

## [1] 0.9978982

sum((df$Price - mean(df$Price))^2) # 偏差平方和

## [1] 17059.73

var(df$Price) # 分散

## [1] 0.352241

sd(df$Price) # 標準偏差

## [1] 0.5934989

偏差平方和を計算する関数

よく使う計算は関数(function)化すると同じコードを繰り返し書く手間が省ける。この先、偏差平方和の計算を何度かすることになるのでssd()という関数を自作しておくことにする。

ssd <- function(x){
  s <- sum((x - mean(x))^2)
  return(s)
}

# Priceの偏差平方和
ssd(df$Price)

## [1] 17059.73

# Distanceの偏差平方和
ssd(df$Distance)

## [1] 2760767

Type (住宅タイプ)別の Price (住宅価格)の分布

df %>% 
  ggplot(aes(x=Price, fill=Type)) +
  geom_histogram(binwidth = 0.1) +
  coord_cartesian(xlim=c(0, 5)) +
  facet_grid(Type~.)

Type別にPriceの分布をみると、“h”という住宅タイプは低価格から高価格まで幅広く分布しているが、住宅タイプ“t”と“u”は低価格帯が多くばらつきは小さい。

df %>%
  group_by(Type) %>%
  summarise(
    平均 = mean(Price),
    標準偏差 = sd(Price),
    偏差平方和 = ssd(Price)
  )

## # A tibble: 3 x 4
##   Type   平均 標準偏差 偏差平方和
##   <fct> <dbl>    <dbl>      <dbl>
## 1 h     1.11     0.638     13900.
## 2 t     0.911    0.367       671.
## 3 u     0.630    0.286       760.

Type別のPriceの偏差平方和を全て足し合わせると\(13900 + 671 + 760 = 15331\)となり、Type別にしない場合の偏差平方和\(17059\)よりも小さくなっている。

このように、Typeという変数を含めてPriceをみるとTypeを含めずにみた場合よりもばらつきが小さくなる。

「なぜPriceにはこんなにばらつきがあるのだろうか？」という問いに対しては「価格が異なる複数の住宅タイプがデータに含まれているから」と言うことができる。これが「Typeという変数でPriceという変数のばらつきを説明する」ということ。

ここで「Typeという変数は“どの程度”Priceのばらつきを説明しているか？」と定量的に考えてみたい。それにはTypeを導入したことによって偏差平方和がどの程度減少したかを計算すればよい。Typeを含まないときの偏差平方和は\(17059\)、Typeを考慮したときの偏差平方和は\(15331\)で、減少分は\(17059 - 15331 = 1728\)となる。これは変数Typeの「説明力」と言うことができる。

<補足>散布図でみるType別のPriceのばらつき

df %>% 
  ggplot(aes(x=Distance, y=Price, color=Type)) +
  geom_point() +
  geom_hline(yintercept = mean(df$Price)) +
  geom_hline(yintercept = mean(subset(df, Type=="h")$Price), color="red") +
  geom_hline(yintercept = mean(subset(df, Type=="t")$Price), color="green") +
  geom_hline(yintercept = mean(subset(df, Type=="u")$Price), color="blue") +
  coord_cartesian(xlim=c(0,50), ylim=c(0,5))

Distance (中心市街地からの距離)別の Price (住宅価格)の分布

前にみたようにDistanceが大きくなると住宅価格が安くなる傾向がある。そこで、Distanceが大きいデータと小さいデータに二等分して、それぞれの偏差平方和を計算してみる。

threshold <- 30 # Distanceで分割する閾値
temp1 <- df %>% filter(Distance<threshold)  # Distanceが小さい住宅のデータ
temp2 <- df %>% filter(Distance>=threshold) # Distanceが大きい住宅のデータ

ggplot(mapping = aes(x=Distance, y=Price)) +
  geom_point(data=temp1, color="lightblue", alpha=0.3) +
  geom_point(data=temp2, color="pink", alpha=0.3) +
  geom_segment(aes(x=0, xend=threshold, y=mean(temp1$Price)), yend=mean(temp1$Price), color="blue") +
  geom_segment(aes(x=threshold, xend=50, y=mean(temp2$Price)), yend=mean(temp2$Price), color="red") + 
  geom_vline(xintercept = threshold, linetype='dashed') +
  coord_cartesian(xlim=c(0,50), ylim=c(0,5))

Distanceを考慮しない場合のPriceの偏差平方和

ssd(df$Price)

## [1] 17059.73

Distanceの大きさで分割した場合のPriceの偏差平方和

ssd(temp1$Price) + ssd(temp2$Price)

## [1] 16831.68

このように、Distanceという連続変数を使うことでもPriceのばらつきを説明することができる。

Part II: 機械学習の基礎〜決定木とランダムフォレスト〜

ここまでは丁寧にヒストグラムや散布図を描きながら、変数のばらつきを別の変数を使って説明できるということをみてきた。しかし、変数の数が多くなると単純な集計によってばらつきの説明を続けるのが難しくなってくる。機械学習を使うと多変数の複雑な関係性も効率的に扱うことができるようになる。今回はシンプルな機械学習アルゴリズムの一つである決定木を実行するプログラムを書いてみる。

決定木の学習

{rpart}というライブラリで決定木の学習をする。「rpart」はRecursive Partitioning and Regression Treesのこと。~の左側に予測したい変数(目的変数)、右側に予測に使う変数(説明変数)を記述する。まずはTypeだけを説明変数にして決定木を学習する。下記のコード一行で決定木の学習が実行される。

tree_model = rpart(Price ~ Type, data = df)

学習済みの決定木の情報を次のように表示できる。しかし、これでは何がどうなったのかわからない。

tree_model

## n= 48433 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
## 1) root 48433 17059.7300 0.9978982  
##   2) Type=u 9292   760.4334 0.6301053 *
##   3) Type=h,t 39141 14743.9600 1.0852120  
##     6) Type=t 4980   671.0723 0.9111480 *
##     7) Type=h 34161 13900.0100 1.1105870 *

決定木の構造を可視化

学習済みの決定木の構造を見やすく可視化する。rpart.plot()という関数にモデルを渡すと自動的に決定木を描画してくれる。

一番上のノードに「100%」と書かれているのはここに全てのデータサンプルが含まれていることを示す。「1」と書かれているのはこのノードに含まれるデータサンプルの住宅価格Priceの平均値を表している(小数点以下が丸められているので紛らわしいが実際は0.998という値が入っている)。

その下の部分にType=uと書かれている。このノードのサンプルのうち、Typeが“u”のもの(“yes”)は左側のノードに、それ以外(“no”)は右側のノードに分割される。左側のノードには全体の19%のサンプルが入り、右側中段のノードには81%のサンプルが入ることになる。

ここで決定木がやっていることは前章の<補足>の図と同じで、Type別に価格の平均値を計算してそれを予測値としている。

rpart.plot(tree_model)

残差平方和と決定係数

偏差平方和(再掲)：変数の平均値とデータとの差の二乗を全てのデータについて足し合わせたもの。

\[ SSD = \sum_{i=1}^{N}(y_i-\bar{y})^2 \] 住宅価格の偏差平方和

ssd(df$Price)

## [1] 17059.73

残差平方和：モデルの予測値とデータとの差の二乗を全てのデータについて足し合わせたもの。残差平方和はモデルが説明しきれなかったばらつきの大きさとして解釈することができる。

\[ RSD = \sum_{i=1}^{N}(y_i-\hat{y_i})^2 \]

Rでは学習済み機械学習モデルの残差を関数residual()で取得することができる。これを二乗して総和をとると残差平方和になる。

sum(residuals(tree_model)^2)

## [1] 15331.51

決定木を使うことで残差平方和は前の章で集計をしていたときよりもだいぶ小さくなり、ばらつきの多くが説明できるようになってきた。しかし、この残差平方和の値そのものは大きいとも小さいとも解釈のしようがない。

ここで、下式で表される決定係数\(R^2\)という指標を導入する。決定係数は元々変数が持っているばらつきのうちどのくらいの割合をモデルが説明できているかを表す。すべてのサンプルに対して平均値で予測するモデルの決定係数は0である。一方、全てのばらつきを説明できるような完璧なモデルであれば決定係数は1になる。

\[ \large R^2 = 1 - \frac{モデルが説明できなかったばらつき(=残差平方和)}{変数のばらつき(=偏差平方和)} \]

上記の決定木モデルの決定係数は次のように計算できる。このモデルが説明できたばらつきは全体の10.1%程度である。

SSD <- ssd(df$Price) # 偏差平方和
RSD <- sum(residuals(tree_model)^2) # 残差平方和
R2 <- 1 - RSD/SSD # 決定係数
R2

## [1] 0.1013041

複雑な決定木

次にType、Rooms、Distance、Method、Regionnameの5変数を使って決定木をつくる。説明変数が複数ある場合は+で区切って記述する。rpart.plot()で構造を可視化すると、先ほどの決定木よりも多数の分岐がある複雑なモデルになっていることがわかる。決定木は数値計算によって「最もデータのばらつきを減らすことができる変数と値の組み合わせ」を探し出してノードを分岐させる。

tree_model = rpart(Price ~ Type + Rooms + Distance + Method + Regionname, data = df)
rpart.plot(tree_model)

変数を増やした分、説明できるばらつきの割合が増えて決定係数が0.558まで大きくなる。

SSD <- ssd(df$Price) # 偏差平方和
RSD <- sum(residuals(tree_model)^2) # 残差平方和
R2 <- 1 - RSD/SSD # 決定係数
R2

## [1] 0.558215

決定木のアンサンブル(ランダムフォレスト)

多変数を同時に考慮するモデルができたとはいえ、決定木では末端ノードの数(上記のモデルでは11個)しか予測値を返すことができない。より精密な予測をするために、多数の決定木を組み合わせるアンサンブルというテクニックを使う。

まず、3つの決定木をつくる。3つの決定木は使用している変数も学習データも異なっている。sample(df, 10000, replace = TRUE)はdfの中から10000個のサンプルをランダムに抽出している(replace = TRUEは一度抽出したサンプルをデータに戻すように指定している。FALSEの場合は一度抽出すると二度と選ばれない)。

set.seed(2018) # 乱数シードを固定（ランダムな抽出を何度実行しても同じ結果が出るようにする）
model1 = rpart(Price ~ Method + CouncilArea, data = sample(df, 10000, replace = TRUE))
model2 = rpart(Price ~ Type + Distance, data = sample(df, 10000, replace = TRUE))
model3 = rpart(Price ~ Distance + Rooms, data = sample(df, 10000, replace = TRUE))

3つのモデルに予測させる住宅データを一つ抽出する。

n <- 2500
df[n,]

## # A tibble: 1 x 13
##   Suburb Address Rooms Type  Price Method SellerG Date               
##   <fct>  <chr>   <int> <fct> <dbl> <fct>  <fct>   <dttm>             
## 1 Abbot… 18 Val…     3 h      1.37 S      Jellis  2018-03-17 00:00:00
## # ... with 5 more variables: Postcode <int>, Regionname <fct>,
## #   Propertycount <int>, Distance <dbl>, CouncilArea <fct>

3つの決定木にそれぞれ予測させてみる。予測値を出力するにはpredictという関数を使う。全く同じデータを入力したにもかかわらず異なる価格を予測していることがわかる。

pred1 <- predict(model1, newdata=df[n,])
pred2 <- predict(model2, newdata=df[n,])
pred3 <- predict(model3, newdata=df[n,])

print(paste(pred1, pred2, pred3))

## [1] "1.11541012302345 1.53371316816709 1.25445517793031"

これら全ての出力値を考慮し、3つのモデルの平均を最終的な予測値として採用することにする。この予測値は一つの決定木から得られたものよりも信頼性が高いと思われる。

pred <- (pred1 + pred2 + pred3) / 3
print(paste(pred))

## [1] "1.30119282304028"

ランダムフォレストはこのように多数の異なる決定木の出力を組み合わせて高精度の予測値をつくる機械学習モデル。考え方は比較的シンプルだが、数ある機械学習モデルの中でも特に高い性能を発揮する万能モデルである。今回は3つの決定木で実験したが、一般的には数百本程度の決定木を組み合わせることが多い。

訓練用データ・テスト用データの分割

データをモデルの訓練用とテスト用に分割する。ひとまずデータの最初の3万件を訓練用、残り(18,433件)をテスト用とする。

df_train <- df[1:30000, ]
df_test <- df[30001:nrow(df), ]

訓練データとテストデータでPriceの分布を比較する。明らかに分布が異なり、訓練データの方が高価格側に偏った分布になっていることが分かる。

ggplot(mapping = aes(x=Price)) +
  geom_histogram(data=df_train, binwidth = 0.1, fill="blue", alpha=0.5) +
  geom_histogram(data=df_test, binwidth = 0.1, fill="red", alpha=0.5) +
  coord_cartesian(xlim=c(0,5))

Priceの分布に差異があるのは元のデータがDistanceの小さい順に並んでいることが原因。実際にDistanceの分布を比較してみると下図のように訓練データにはDistanceが小さいサンプル、テストデータにはDistanceが大きいサンプルしか含まれていないことがわかる。

ggplot(mapping=aes(x=Distance)) +
  geom_histogram(data=df_train, binwidth = 2, fill="blue", alpha=0.5) +
  geom_histogram(data=df_test, binwidth = 2, fill="red", alpha=0.5) +
  coord_cartesian(xlim=c(0,50))

このように、訓練データとテストデータの作り方は機械学習をする上で注意が必要なポイントである。

今回はデータからランダムに3万件をピックアップしてきて訓練データとし、残りをテストデータとする(※)。
※ランダムに分ける方法が常に適切であるとは限らない。

set.seed(2018)
train_index <- sample(nrow(df), 30000) #学習用データをサンプリング
df_train <- df[train_index,]
df_test <- df[-train_index,] # 「-train_index」とするとtrain_index以外のサンプルを抽出する。

ランダムフォレストモデルの構築

Rの{randomForest}ライブラリを使ってランダムフォレストのモデルを構築する。今回はRooms、Distance, Type、Regionname、CouncilArea、Method、Propertycountの7変数を使用する。決定木の数はntreeで指定する。ランダムフォレストの中の各決定木はmtryで指定した数の変数を7つの変数の中からランダムに選んで使用する。

forest_model = randomForest(Price ~ Rooms + Distance + Type + Regionname + CouncilArea + Method + Propertycount,
                            data = df_train,
                            ntree = 30,
                            mtry = 2,
                            importance=TRUE
                            )

ランダムフォレストモデルによる予測と評価

テストデータから一つのサンプルを選んで予測する。

n <- 300
df_test[n,]

## # A tibble: 1 x 13
##   Suburb Address Rooms Type  Price Method SellerG Date               
##   <fct>  <chr>   <int> <fct> <dbl> <fct>  <fct>   <dttm>             
## 1 South… 345 Fe…     3 h       1.2 S      Greg    2018-10-13 00:00:00
## # ... with 5 more variables: Postcode <int>, Regionname <fct>,
## #   Propertycount <int>, Distance <dbl>, CouncilArea <fct>

決定木のときと同様、predict関数を使って予測する。

print(as.double(predict(forest_model, df_test[n,])))

## [1] 1.81593

全テストデータに対してまとめて予測して、結果をデータフレームに格納する。

pred <- predict(forest_model, df_test)
df_test <- df_test %>% 
  mutate(
    Prediction = pred, # 予測結果
    Error = Price - Prediction, # 誤差
    ErrorRate = Error / Price # 誤差率
    )

正解と予測結果の対応をプロットする。

ggplot(df_test, aes(x=Price, y=Prediction, color=Distance)) +
  geom_point()

誤差率の分布をプロットする。誤差率0%付近を中心にして山形の分布になっている。誤差の大きいものは±50%以上のものもある。

df_test %>% 
  ggplot(aes(x=ErrorRate)) +
  geom_histogram(binwidth=0.05, color="grey20", fill="red", alpha=0.7) +
  coord_cartesian(xlim=c(-1.5,1.5))

決定係数を計算すると0.728となった。ランダムフォレストモデルによって住宅価格のばらつきのうち72.8%を説明できたということができる。

SSD <- ssd(df_test$Price) # 偏差平方和
RSD <- sum((df_test$Price - df_test$Prediction)^2) # 残差平方和
R2 <- 1 - RSD/SSD
print(R2)

## [1] 0.7255503

特徴量重要度

importance()という関数でランダムフォレストの変数の重要度を計算することができる。この重要度は各変数がばらつきを減らすのに寄与した度合いで測っている。

importance(forest_model, type = 1)

##                 %IncMSE
## Rooms         23.728743
## Distance      16.759674
## Type          30.687427
## Regionname    10.289144
## CouncilArea   16.723824
## Method         5.109817
## Propertycount 26.501942

Type、Rooms、Distanceなどの変数の重要度が大きく、RegionnameやMethodは重要度が小さいことがわかる。

importance_df <- as.data.frame(importance(forest_model))
importance_df %>% 
  ggplot(aes(x=rownames(importance_df), y=`%IncMSE`)) +
  geom_bar(stat="identity")

おわりに

さまざまな変数を考慮したランダムフォレストのモデルを用いることによりばらつきの大部分を説明することができた。それでもなお残っているばらつきは誤差と呼ばれる。この誤差をさらに小さくするべくモデルを精密にチューニングしていくことも一つの進め方ではあるが、ここでは「なぜこの誤差があるのか？」と考えてみる。

仮に地域も売却方式も部屋の数も住宅タイプも全く同じ住宅が複数存在したとする。これらの住宅の価格は果たして同じだろうか？ある住宅の近くにはスーパーマーケットがあり、他の住宅の近くには墓地があるかもしれない。また、住宅の価格を決める不動産会社の担当者が異なれば違う価格をつけるかもしれない。このようにデータからは見えない要因は無数に挙げることができ、それらの要因によるばらつきは原理的に説明不可能な本質的ばらつきである。このような本質的ばらつきまで無理に説明しようとすると、欠陥のあるモデルが出来上がることになってしまう（今回は触れなかったが過学習と呼ばれる）。

製造現場の異常検知では正常信号のばらつきの範囲に入らないような異常信号を見つけようとするし、書籍のレコメンドでは好きなジャンルのばらつきが小さくなるような変数(年齢・性別・購買履歴など)を選んでユーザーをクラスタリングする。このように、一見全く別のことをしているように思える機械学習の様々な応用も、データのばらつきを説明するための手法を目的に応じて使い分けているとみることができる。