【CodeIQ】Rで解くデータサイエンティスト問題の解説（R Advent Calendar2012）

先日より、リクルート様のITエンジニアのための実務スキル評価サービス「CodeIQ」で、データサイエンティストに関する問題を出題させて頂いております（問題集はこちら）。先日12/12のおしゃスタ＠リクルートでも少し解説しましたが、Rでの解答例をお見せする時間がなかったので、この機会にブログで公開します（おしゃスタに関するCodeIQ様のブログはこちら）。去年に引き続き勢いだけで参加したR Advent Calendar 2012でしたが、ちゃんとネタが見つかって良かった！！！でも無計画に参加したらクリスマスイブの日に当たってしまったので、、、日付が変わるくらいにさっさと書いてしまいたいと思います！！！爆

【データサイエンティスト初級問題】

【前提】
とある転職サイトから、「とりあえずデータがあるんだけど、、、」と言われてデータを受け取りました。先方は何をして欲しいかまだはっきりと決まってない様子。

受け取ったデータは、
「応募した人の属性データ(oubo_zokusei.csv)」
「応募した時間の記録データ(oubo_kiroku.csv)」
の2種のCSVファイルです。

【課題】
幸いにもデータサイズはあまり大きくなく、Rで読み込める程度。まずはRで読み込んでどんなデータか把握してください。
解答はテキストファイル(.txt)で、

1. Rコード
2. データの概要：どんなデータかまとめた文章
3. 分析提案：もし応募数の予測モデルを作るとしたらあった方がよさそうな変数とその理由

の3つを、この順番で書いて下さい。

データサイエンティストによくある状況なんですが、「データがあるからとりあえず何かやってよ」と言われているパターンです。「仮説がはっきりしてないと分析できませんよ」と言いたくなる気持ちを抑えつつも、何か結果を出す必要があります。受け取ったときはデータの定義ファイルなんかもありません。データを読み込んで始めて分かる状況です。解答例としては、例えば以下のようになります。

【解答例】

1. Rコード

# データチェック
setwd("......")

ID   <- read.csv(file("oubo_zokusei.csv", encoding="cp932"), as.is=T, header=F)
oubo <- read.csv("oubo_kiroku.csv", as.is=T, header=F)

head(ID)
head(oubo)

dim(ID)
dim(oubo)

names(ID)   <- c("ID", "join_time", "age", "gakureki", "shokureki")
names(oubo) <- c("ID", "time")
oubo$month  <- months(as.Date(oubo$time))

hist(ID$age)
summary(ID$age)
table(ID$gakureki, useNA="always")
table(ID$shokureki, useNA="always")

table(oubo$month)


# 全期間での応募数
ouboNum  <- aggregate(oubo$ID, list(oubo$ID), length)
names(ouboNum) <- c("ID", "oubo_num")

# 各月の応募数
ouboNumMonth <- aggregate(oubo$ID, list(oubo$ID, oubo$month), length)
names(ouboNumMonth) <- c("ID", "month", "oubo_num")
ouboNumMonth$month <- ifelse(ouboNumMonth$month == "5月", "mon5", 
                      ifelse(ouboNumMonth$month == "6月", "mon6", "mon7"))

library(reshape)
ouboNumMonth1 <- cast(ouboNumMonth, ID~month, value="oubo_num")
ouboNumMonth1 <- as.data.frame(ouboNumMonth1)
ouboNumMonth1[is.na(ouboNumMonth1)] <- 0

# 入会データとマージ
ouboID <- merge(ID, ouboNum, by="ID", all.x=T)
ouboID <- merge(ouboID, ouboNumMonth1, by="ID", all.x=T)
ouboID[is.na(ouboID)] <- 0

# 可視化
cor(ouboID[, 7:9])
plot(ouboID[, 7:9])

# モデル
ouboLM <- lm(mon7~age+gakureki+shokureki+mon5+mon6, data=ouboID)
summary(ouboLM)

cor(predict(ouboLM), ouboID$mon7)
plot(predict(ouboLM), ouboID$mon7, xlim=c(0, 100), ylim=c(0, 100))
abline(a=0, b=1)


library(randomForest)
ouboRF <- randomForest(mon7~mon5+mon6, data=ouboID)
summary(ouboRF)

cor(predict(ouboRF), ouboID$mon7)
plot(predict(ouboRF), ouboID$mon7, xlim=c(0, 100), ylim=c(0, 100))
abline(a=0, b=1)

2. データの概要レポート

zokuseiデータ：恐らくID・登録日・年齢・学歴・勤務状況が記録されている
ouboデータ：恐らく誰が・いつ応募したのか記録されている
20代〜30代、大卒が中心に利用。離職中が在職中の約半分。

とりいそぎ過去の応募数から次月の応募数を予測するモデルを作ってみたところ、単純な線形回帰モデルでR二乗が約0.33であった。ランダムフォレストも同じくらいの精度であった。次月の応募数のうち約３割が説明可能という解釈ができるが、プロットを確認すると応募数が多い外れ値に引きずられているようにも見える。

3. 分析提案

詳細な分析を行うためには、属性データの詳細が必要。またサイトのアクセス履歴も分析することで、ユーザーの行動特性が見つかるかもしれない（ただ、アクセス解析はRではなく専用ソフトが向いているかも）。
将来の応募数を予測することで、営業支援や売上予測につなけることも可能である。そのためには現在のモデルより精度を上げなければならないため、新たな変数が必要であるし、モデルも吟味する必要がある。

【解説】

だいたいここまで分析時間と文章に起こす時間で、１〜２時間といったところです（エラーが出たりしたらもっとかかりますが）。むしろここまでは割とさくっとできるのですが、ここから先が大変です。仮説をヒアリングしたり、この結果を丁寧に説明したり。丁寧に説明しなくても分かる人が増えてくれれば良いんですがね。マッキンゼーが「データサイエンティストをマネジメントする人材が150万人足りなくなる」って言っているのはこのことだと感じてます（詳しくはこの記事など）。あと実務でこういう事をやろうとすると、データを作るまでに結構な時間がかかってしまったりします。分析に慣れている方がDBを操作していると、簡単に作ってくれたりすることもあります。

しかし、出題した後に、結構レベルが高かったかな〜と思ったり。。まだみなさんの回答を見れてないのですが、むしろこんな出題の仕方で回答してもらえて恐縮でした。。。初級の次に中級・上級を予定してたのですが、それは断念しました。

その反省点を活かして、次の問題からは課題をはっきりとさせました。データサイエンティストで必要な学問知識のうち重要な「統計学」と「機械学習」に分けて出題しています。Codingして結果を出せるように、Rを使って回答してもらう、という形式にしています。以下問題と、私の解答例を紹介します。

データサイエンティスト〜統計学編1〜

【問題】
問1. Rを使い、5人分の身長に関し、下記の6つの統計量を求めてください。

　# 5人分の身長データ
　height1 <- c(168, 173, 152, 181, 175)

　(1) 5人分の身長データの平均値
　(2) 5人分の身長データの中央値
　(3) 5人分の身長データの標本分散
　(4) 5人分の身長データの不偏分散
　(5) 5人分の身長データの標準偏差（不偏分散を使う）
　(6) 5人分の身長データの標準誤差（不偏分散を使う）

問2. 下記のrnorm関数を使うと5人分の身長のデータを擬似生成できます。擬似生成データを使った以下の質問に答えてください。

　# 5人分の身長を擬似生成するコード
　height2 <- 170 + 10*rnorm(5)

　(1) 擬似生成した5人分の身長データを使って
　　　- 標本分散
　　　- 不偏分散
　　　- 標準偏差（不偏分散を使う）
　　　- 標準誤差（不偏分散を使う）
　　の4つの統計量を計算してください。
　　計算に使ったRのコードと計算結果を提出してください。

　(2) 「100人分の身長データ」、「1,000人分の身長データ」を擬似生成し、上記4つの統計量（標本分散、不偏分散、標準偏差、標準誤差）を計算してみてください。人数（データ件数）が増えたとき、この4つの関係性はどうなるか議論してください。

【解答例】

# --1.
height1 <- c(168, 173, 152, 181, 175)
mean(height1)
median(height1)
sum( (height1 - mean(height1)) ** 2)/length(height1)
var(height1)
sd(height1)
sd(height1)/sqrt(length(height1))

# --2.
stat <- function(vec){
  a <- sum( (vec - mean(vec))**2)/length(vec)
  b <- var(vec)
  c <- sd(vec)
  d <- sd(vec)/sqrt(length(vec))
  c(a, b, c, d)
}

height2 <- 170 + 10*rnorm(5)
stat(height2)

height3 <- 170 + 10*rnorm(100)
stat(height3)

height4 <- 170 + 10*rnorm(1000)
stat(height4)

# 標本分散と不偏分散はほとんど等しくなっていく
# 標準偏差は10付近をばらついている
# 標準誤差はゼロに近づいていく

データサイエンティスト〜統計学編2〜

【問題】
問1. DataScience_stat2.csvのデータは100万人分の身長を模擬的に作成したものです。
この100万人分の身長データから、1000人分のデータをRを使ってサンプリングしてください。

問2. 問1でランダムサンプリングしたデータを使って以下の値を求めてください。
　　 2-1. 平均値
　　 2-2. 標準誤差（不偏分散を使う）
　　 2-3. 平均値の95%信頼区間

問3. 問2-3の信頼区間は、もとの100万人の集団の平均身長（真値）を推測するものです。サンプリングを無限に繰り返すと、理論的には95%の割合で信頼区間の範囲に真値が含まれます。ランダムサンプリングを10000回繰り返し、この事を確認してください。計算に使ったRのコードと簡単な説明を提出してください。

※DataScience_stat2.csvは以下のRプログラムによって作成しています。
set.seed(1)
height <- 170 + 10*rnorm(1000000)

【解答例】

# --1.
height <- read.csv("DataScience_stat2.csv", as.is=T)

mean(height$x)
sd(height$x)
sd(height$x)/sqrt(length(height$x))


# --2.
set.seed(1)
height_sample <- sample(height$x, 1000)

MEAN <- mean(height_sample)
SE   <- sd(height_sample)/sqrt(length(height_sample))
LCL  <- MEAN - 1.96*SE
UCL  <- MEAN + 1.96*SE


# --3.
CI <- data.frame(LCL=NULL, UCL=NULL)
for(i in 1:10000){
  set.seed(i)
  height_sample <- sample(height$x, 1000)
  
  MEAN <- mean(height_sample)
  SE   <- sd(height_sample)/sqrt(length(height_sample))
  LCL  <- MEAN - 1.96*SE
  UCL  <- MEAN + 1.96*SE
  
  CI <- rbind(CI, c(LCL, UCL))
  print(i)
}
names(CI) <- c("LCL", "UCL")

# ----真値が信頼区間に含まれている割合
CI$BetweenCI <- ifelse(CI$LCL <= mean(height$x) & mean(height$x) <= CI$UCL, 1, 0)
table(CI$BetweenCI)
prop.table(table(CI$BetweenCI))

# ----グラフでチェックする
plot(c(168, 172), c(1, 100), type="n")
for(i in 1:100){
  lines(c(CI$LCL[i], CI$UCL[i]), c(i, i), col=2-CI$BetweenCI[i])
}
abline(v=mean(height$x), lty=2)

データサイエンティスト〜機械学習編1〜

【問題】
問1. Rを使い、DataScience_ML1.csvを読み込み以下の線形回帰モデルを作成してください。
　　 y=x1+x2

問2. Rを使い、問1で作ったモデルに対して以下の回帰診断を行なってください。

　　(1) ローデータの散布図
　　(2) 調整済みR二乗
　　(3) 残差プロット
　　(4) キャリブレーションプロット

※DataScience_ML1.csvは以下のRプログラムによって作成しています。
set.seed(1)
x1 <- rnorm(10000)
set.seed(2)
x2 <- rnorm(10000)
set.seed(3)
y <- 2*x1 + x2**2 + rnorm(10000)
Data <- data.frame(x1 = x1, x2 = x2, y = y)

【解答例】

# --1.
Data <- read.csv(DataScience_ML1.csv, as.is=T)
LM1  <- lm(y ~ x1 + x2, data=Data)

# --2.
plot(Data)
summary(LM1)
plot(LM1)
plot(predict(LM1), Data$y)

# --3.
Data2 <- data.frame(Data, x3 = Data$x2**2)
LM2   <- lm(y ~ x1 + x2 + x3, data=Data2)
summary(LM2)
plot(predict(LM2), Data2$y)

データサイエンティスト〜機械学習編2〜

【問題】
問1. Rを使い、DataScience_ML1.csvを読み込み、以下のモデルを作成してください。結果変数はyとします。

・線形カーネルのSVRモデル
・3次多項式カーネルのSVRモデル
・ガウシアンカーネルのSVRモデル（radial basis）
・シグモイドカーネルのSVRモデル

問2. 作ったモデルに対して予測診断を行うために、各モデルのキャリブレーションプロットとR二乗をチェックしてください。

【解答例】

# --1.
library(e1071)
Data <- read.csv(DataScience_ML1.csv, as.is=T)
SVR1 <- svm(y~., data=Data, type="eps-regression", kernel="linear")
SVR2 <- svm(y~., data=Data, type="eps-regression", kernel="polynomial", degree=3)
SVR3 <- svm(y~., data=Data, type="eps-regression", kernel="radial")
SVR4 <- svm(y~., data=Data, type="eps-regression", kernel="sigmoid")


# --2.
plot(predict(SVR1), Data$y)
plot(predict(SVR2), Data$y)
plot(predict(SVR3), Data$y)
plot(predict(SVR4), Data$y)

cor(predict(SVR1), Data$y) ** 2
cor(predict(SVR2), Data$y) ** 2
cor(predict(SVR3), Data$y) ** 2
cor(predict(SVR4), Data$y) ** 2

以上が解答例になります。解説はまた今度の機会に行おうと思っています。CodeIQのサイトからデータのダウンロードはできなくなっていますが、統計学編と機械学習編は各自の環境でデータを作成できます。ぜひやってみて下さい。今回紹介した5問は分析では入り口ですが、これらができるだけでも実務ではとても役に立ちます。これくらい出来ていれば、他の手法も理解しやすくなると思います。おしゃスタでの発表資料でも問題の解説を少ししていますので、補助的に見て下さい。

おしゃスタ@リクルート from Issei Kurahashi

今後も分析力を身に着けてレベルアップしていけるように、問題をステップアップできればと思っています。最終的には私の修士論文とか博士論文のネタに辿り着ければな〜と思ったり（ここに置いてます）。

それでは、今回の記事はこのへんで失礼致します。

みなさまよいクリスマスと年末年始をお過ごし下さい！！！！！！！！！！