Kaggleデータ探索(Backorder Data)

前回「Kaggle dataをあさってみる」で、Human Resource Analysisのデータを扱ったのですが、http://ritsuan.com/blog/5840/　あの時はまだ人気6-7番目くらいだった今ではKaggleである種一番人気のある(hotnessでsort)データとなっています。

Google検索で、HumanResourceAnalysisと入力すると僕が書いた記事がトップ3に来るので、僕の記事がダウンロード数アップに貢献したのかもしれません笑。

引き続き、Kaggleのデータセット採掘作業頑張りたいと思います。
今回は、Kaggle datasetsのページで、Hotnessでsortして、上から34番目にある(2017/05/10)
「Can You Predict Product Backorders?」に興味をそそられたので、これを見てみることにします。

https://www.kaggle.com/datasets?sortBy=hottest&group=featured

Back orderとは「取り寄せ注文」のことを言います。
さて、データをダウンロードしてみますと、、
解凍したフォルダには、Kaggle_Training_Dataset_v2.csvファイルと、Kaggle_Test_Dataset_v2.csvファイルがが入っています。
それぞれ、僕のwork directory(Documentsディレクトリ)に移して、中身を見ていきたいと思います。
#readrパッケージ読み込み
library(readr) #もしinstall.packagesしてなければしてください。
#データの読み込み
d=read_csv(“Kaggle_Training_Dataset_v2.csv”)
これを実行すると、1687861行目においてparsing failuresを吐き出すので、この1行を消去するときれいになります。(実際この列はすべて欠損値です)
d=d[-1687861,]
d

パッと見だと2列目のlead_time列が一つ置きに欠損値のようです。。(下の方は違うみたいですが)
str(d)
#168万列ありますね。。
このデータの説明ページ(https://www.kaggle.com/tiredgeek/predict-bo-trial)
にある各変数の説明を和訳してみます。

sku – 製品のランダムID
national_inv – 部品の現在の在庫レベル
lead_time – 製品の移動時間（利用可能な場合）
in_transit_qty – ソースからの輸送中の商品の量
forecast_3_month – 今後3ヶ月間の売上予測
forecast_6_month – 今後6ヶ月間の売上予測
forecast_9_month – 今後9ヶ月間の売上予測
sales_1_month – 前1か月の販売数量
sales_3_month – 前3か月の販売数量
sales_6_month – 前6か月の販売数量
sales_9_month – 前9ヶ月の販売数量
min_bank – 在庫の最小推奨量
potential_issue – 特定された部分の商品発送元の問題
pieces_past_due – 商品発送元から延滞されているパーツ
perf_6_month_avg – 過去6か月間の商品発送元の成果
perf_12_month_avg – 過去12か月間の商品発送元の成果
local_bo_qty – 延滞されている在庫発注量
deck_risk – パーツリスクフラグ
oe_constraint – パーツリスクフラグ
ppap_risk – パーツリスクフラグ
stop_auto_buy – パーツリスクフラグ
rev_stop – パーツリスクフラグ
went_on_backorder – 実際に取り寄せ注文になりつづけた商品。これが目的変数。

それでは各変数の要約統計量を見てみましょう。
summary(d)
lead_timeに欠損値がありますね。。
次に質的変数の要約統計量を見てみましょう。
このままだと扱いにくいので、data.frame型にします。
#データフレーム型への変換
t=data.frame(d)
#文字列型をFactor型へ変換
for(i in c(13,18,19,20,21,22,23)){t[,i]=as.factor(t[,i])}
#結果確認
str(t)

#統計量の算出
summary(t)
ということで質的変数の7列は、すべてYesかNoの二択であることがわかります。
さて次に3列目のlead_time列にある欠損値の補完を行いましょう。10万個くらい欠損値がありますね。
miceパッケージで欠損値の補完を試してみたのですが、10万個の欠損値補完はmiceパッケージではきつかったみたいで、methodとして、fastpmmとrf(randomforest)を試してみたのですが、前者ではエラー、後者では計算が20分経っても終わらなかったので断念しました。
なので、欠損値を目的変数にして、機械学習で欠損値を補完しようと思います。
まずデータを分割します。
a=complete.cases(t)
head(a)

#学習用
train=t[a, ]
#欠損値の入ったデータを切り分ける
library(dplyr)
test=filter(t, !(sku %in% train$sku))
#中身確認
str(train)
#中身確認
str(test)
これでうまく分けられましたね。
ちなみにですが、NAが入っているデータと入っていないデータを切り分けるのは実は大変です。思いつく限りのいろんな方法で試してみれば、一筋縄でいかないことがわかるかと思います。上のやり方は数少ないうまくいった方法です。
僕が大好きなrandomforestで学習器を作成します。
#ランダムフォレストパッケージのロード
library(randomForest)
#CPUコアをすべて使用する命令
library(doParallel)
registerDoParallel(detectCores())
#普通にやったらメモリが5.9Gb必要でできなかったのでorzデータを1/100サンプリング
ind=sample(nrow(train),floor(nrow(train)/100))
train_new=train[ind,]
dim(train_new)
 #説明変数の作成
train_x=train_new[,c(2,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23)]
#目的変数の作成
train_y=train_new[,3]
#ランダムフォレストの実行(数十分かかるので、C++で書かれているrangerパッケージやRboristパッケージを使った方がいいかもです。Rboristパッケージなら、rf_result=Rborist(train_x,train_y)でできますが、rangerパッケージだとformulaしか受けつけずひと手間かかります。)
#Rborist(https://cran.r-project.org/web/packages/Rborist/Rborist.pdf)
#ranger(https://cran.r-project.org/web/packages/ranger/ranger.pdf)
#今回はrandomForestパッケージで実行します。
rf_result=randomForest(train_x,train_y)
#モデル保存(一応保存しておきます)
saveRDS(rf_result,”rf_result.rds”)
#結果
rf_result

#テスト用の説明変数の作成
test_x=test[,c(2,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23)]
#欠損値の予測
prediction=predict(rf_result,test_x)
str(prediction)

#欠損値の補完
test[,3]=prediction

欠損値の補完ができたので次は、今までのtrainデータとtestデータを結合させましょう。
tableの結合には、dplyrパッケージのfull_join関数を使用します。
full_data=full_join(train,test)
str(full_data)
このデータ数1687860個は、はじめ読み込んだ個数と同じです。うまくいっていますね。
では次に欠損値がないか見てみましょう。
summary(full_data)
欠損値はないですね。
欠損値補完はうまくいったようです。
#データ保存(一応保存しておきます)
write.csv(full_data,”full_data.csv”)

さて、もともとこのデータは、目的変数が、23列目のwent_on_backorder列です。
これを目的変数(応答変数)として、分類モデルを構築することがこのデータのそもそもの目標です。
のでやっていきましょう。
せっかくなので、Rboristパッケージを使用して、ランダムフォレストによる分類を行っていきたいと思います。。。と言いたいところなんですが(実際に試したのですが)メモリが少ないためか(4GB)10回くらいRが落ちたので今回は断念します(スペックの高いPCを手に入れたら再挑戦します)。
学習用データ:テスト用データ=80:20でやっていこうと思います。
id列は使わないので、省きます。
#id列消去
full_data1=full_data[,2:23]

また目的変数を見てみると、Noが1676567個、Yesが11293個であり、二桁違うことがわかります。
不均衡データへの対応は,現在の実行環境がメモリが少なく、ダウンサンプリングで実行してみます(SMOTE,ROSE,upsamplingのいづれもPCが悲鳴をあげたため)。

#データ分割
set.seed(123)
ind=sample(nrow(full_data1),floor(nrow(full_data1)*0.8))
train1=full_data1[ind,]
test1=full_data1[-ind,]
dim(train1)
dim(test1)

#ダウンサンプリング
library(caret)
set.seed(1357)
down_train = downSample(x=train1[,-22],y=train1[,22])
table(down_train[,22])

#↑これでNoとYesが同じ数になり、不均衡が改善しました！
#ランダムフォレストパッケージによるランダムフォレスト
library(randomForest)
rf=randomForest(down_train[,-22],down_train[,22],importance=T)
rf_pred=predict(rf,test1[,-22])

#混合行列
a=table(rf_pred,test1[,22])
a

#Accuracy
accuracy=(a[1,1]+a[2,2])/(a[1,1]+a[1,2]+a[2,1]+a[2,2])
accuracy
Accuracyが0.88あればとりあえずよいでしょう。一応モデル構築はできました。
これで一応今回の目標は達成しました。
一応説明変数の重要度を見てみましょう。
varImpPlot(rf)
とりあえずnational_inv (現在の在庫レベル)が重要なのはわかりました。データの意味を考えると、まあ、確かにな。。といった感じです。
左右の表の順序があまり一致していないのでまだ安定したモデルとは言えないようですね。ハイパーパラメータチューニングしていなかったりdownsamplingで多い方のラベルのデータを切り捨てたことが原因でしょう。

ちなみに最後に、データをダウンロードした際にフォルダに入っていたもう一つのデータ、
Kaggle_Test_Dataset_v2.csv　の中身を見て終わりにしましょう。
c=read.csv(“Kaggle_Test_Dataset_v2.csv”,header=T)
str(c)
どうやらデータ数が一桁少ないだけで、目的変数のwent_on_backorderもあり、これをテストデータにすべきでしたね。。モデル保存した学習器を使って欠損値を埋めてから、予測モデルを適用すればよいようです。今回はやりませんが。

ということで、内容が面白そうなので手を出したBackorderデータでしたが、内容にはほとんど踏み込むことなく(内容の詳細解析以外のところで多大な労力がかかった)、データ数が100万件以上あり欠損値が10万個以上ということもあり、
“ビックデータ”、”欠損値対応”、”不均衡データ対応”、”corei3+メモリ4GBでのロースペックでの解析”がテーマとなる、まさかのたいへんな解析例になってしまったのですが、データ解析を実際にやっている方には、見ごたえがあった方もいらっしゃったかなと思います。
今回はここで、このデータの解析を終えたいと思います。

鈴木瑞人
東京大学大学院新領域創成科学研究科　メディカル情報生命専攻　博士課程
東京大学機械学習勉強会　代表
NPO法人Bizjapan テクノロジー部門BizXチームリーダー