Rで機械学習モデルを構築する方法

個人的にはpythonが得意なのですが、Rの復習も兼ねて簡単にデータから機械学習を構築するまでの流れを追ってみました。
機械学習初学者やpython使いだけどRを勉強したいという方の参考になれば幸いです。
(Kaggleのカーネルを参考にしています。)

データ確認

今回は不動産価格の予測を行います。
データは以下のgithubからcsvファイルをダウンロードしてください。
handson-ml/datasets/housing at master · ageron/handson-ml · GitHub
ちなみにこのデータ、私の大好きなオライリーの本で使用しているデータです。

scikit-learnとTensorFlowによる実践機械学習

それではデータを触っていきます。
まずはライブラリーの読み込みと、csvファイルからデータの読み込み。
データは地域ごとに存在する住宅の情報です。ここから、新しい地域の住宅情報を与えられた時に、住宅価格(median_house_value)を予測するモデルを構築します。住宅価格を目的変数と呼び、目的変数予測のために使う他の変数を説明変数と呼びます。

library(tidyverse)
library(reshape2)

housing <- read.csv('housing.csv')
head(housing)

次にデータの概要を確認します。

summary(housing)

ここからtotal_bedroomsにはNAが入っているということと、ocean_proximityがカテゴリーデータになっていることが分かります。

次にデータを可視化します。ちょっと軸の数字が重なっていますが気にしないでください。。

par(mfrow=c(2,5))
ggplot(data = melt(housing), mapping = aes(x = value)) + 
  geom_histogram(bins = 30) + facet_wrap(~variable, scales = 'free_x')

データクレンジング

次にデータクレンジングをします。
total_bedroomsにデータの欠損値が入っていることが確認できたのでここは同じカラムの値の中央値で埋めることにします。(王道のやり方なのであまり深く考えずに処理します。本当は、中央値で埋めていいかどうかは検討が必要。)

housing$total_bedrooms[is.na(housing$total_bedrooms)]<- median(housing$total_bedrooms , na.rm = TRUE)

次に世帯当たりのベッドルームの数と部屋の数を新しくデータとして作成します。
これは地域に存在する世帯数によって部屋の数が異なるのは当たり前のことで、世帯ごとの部屋の数のデータにしか価値がないからです。

housing$mean_bedrooms<- housing$total_bedrooms/housing$households
housing$mean_rooms<- housing$total_rooms/housing$households
drops<-c('total_bedrooms', 'total_rooms')
housing<-housing[ , !(names(housing) %in% drops)]
head(housing)

次にカテゴリのデータをone hot encodingします。
なんだそれという方、ググっていてください。すぐに分かります。

簡単に説明すると、文字列などでカテゴリに別れているデータはそのまま機械学習モデルに突っ込むことができないので、数値化してあげる必要があります。その時に、Aグループは1に変換、Bグループは2に変換、Cグループは3に変換、、、という具合に数値変換するとグループ間に重みの偏りが出てしまうので、Aグループかどうか(0 or 1)、Bグループかどうか(0 or 1)、という具合に各グループごとにデジタル値として変換してあげます。

pythonだとpandasにget_dummiesというメソッドがあるのでそのまま使えますが、Rはどうすればいいのか分からなかったのでチュートリアル通りに写経します。(あとで調査します。)
やっていることは簡単なのでコードを追ってみてください。
カテゴリーごとにカラムを生成してそこに0か1が入ったデータフレームが作成されます。

categories<-unique(housing$ocean_proximity)
cat_housing<-data.frame(housing$ocean_proximity)
colnames(cat_housing) <- 'ocean_proximity'

for(cat in categories){
  cat_housing[,cat]<-rep(0, times<-nrow(cat_housing))
}

for(i in 1:length(cat_housing$ocean_proximity)){
  cat<-as.character(cat_housing$ocean_proximity[i])
  cat_housing[,cat][i]<-1
}

head(cat_housing)

もともとのカテゴリ情報が入っていたocean_proximityは不必要になるので削除します。

cat_columns<-names(cat_housing)
keep_columns<-cat_columns[cat_columns != 'ocean_proximity']
cat_housing<-select(cat_housing,one_of(keep_columns))
head(cat_housing)

次に説明変数のスケーリング(標準化)を行います。スケーリングは外れ値や変動の大きい変数に予測結果が引っ張られないように、値を標準化します。

その前に説明変数として使用しないカラム(median_house_value目的変数,ocean_proximity使わない)をマスターデータフレームから削除します。その後にスケーリング。

drops<-c('ocean_proximity','median_house_value')
housing_num <-housing[ , !(names(housing) %in% drops)]

scaled_housing_num <- scale(housing_num)
head(scaled_housing_num)

スケーリングが完了したので、目的変数(median_house_value)とカテゴリーのone hotデータフレームを合体させます。

cleaned_housing <- cbind(cat_housing, scaled_housing_num, median_house_value=housing$median_house_value)
head(cleaned_housing)

モデリング

さてデータが整ったのでモデリング開始します。
機械学習モデルは通常学習データとテストデータに分けて、学習データを用いてモデルを作り、テストデータを使ってその精度を検証するというやり方をします。なので、ここでもまずデータを学習データとテストデータに分割します。

set.seed(1738)
sample<-sample.int(n = nrow(cleaned_housing), size<-floor(.8*nrow(cleaned_housing)), replace = F)
train<-cleaned_housing[sample, ] 
test<-cleaned_housing[-sample, ] 

それではモデリング＋検証を行います。
まずはGLMからです。コードではクロスバリデーションを使っています。詳しい説明はこちらを参考にしてください。

library('boot')
glm_house <- glm(median_house_value~median_income+mean_rooms+population, data=cleaned_housing)
k_fold_cv_error<-cv.glm(cleaned_housing , glm_house, K=5)
names(k_fold_cv_error)

>"call" "K" "delta" "seed"

クロスバリデーションの結果(RMSE)は以下のようになります。

glm_cv_rmse<-sqrt(k_fold_cv_error$delta)[1]
glm_cv_rmse

>83396.67

次にランダムフォレスト。先ほど分割した学習データとテストデータを使用します。
学習データを目的変数と説明変数に分けます。

library('randomForest')
names(train)
set.seed(1738)

train_y<-train[,'median_house_value']
train_x<-train[, names(train) !='median_house_value']

ランダムフォレストの場合クロスバリデーションをしなくても、ランダムに学習に使っていないデータ(out-of-bag (oob) )が残っているのでそのデータを用いいて精度検証できます。

rf_model<-randomForest(train_x, y = train_y , ntree = 500, importance = TRUE)
oob_prediction <- predict(rf_model)
train_mse <- mean(as.numeric((oob_prediction - train_y)^2))
oob_rmse <- sqrt(train_mse)
oob_rmse

>49148.28

最後にテストデータを使ってスコアリングをします。

test_y = test[,'median_house_value']
test_x = test[, names(test) !='median_house_value']
y_pred = predict(rf_model , test_x)
test_mse = mean(((y_pred - test_y)^2))
test_rmse = sqrt(test_mse)
test_rmse

>48273.99

この結果が良いか悪いかは求める精度次第ですので、ビジネスを始める時に目標精度定めておく必要があります。

簡単にモデリングまでの流れを追ってみました。
以上です。