例のデータ解析3: データ分析プロセス4章を添えて
@siero5335
2015年7月13日
以下の発表で紹介したデータについて解析したフローを示します。
Tokyo.Rでの発表スライドはこちら 。
今回は前々回、前回の解析に、ご恵贈頂いたデータ分析プロセス４章に記載されている
予測モデルの構築の手法を盛り込んだ内容です。
はじめにdevtools::install_github(“dichika/jaguchi”)で
jaguchiパッケージをインストールしておきます。
パッケージ読み込み
library(jaguchi)
library(dplyr)
library(caret)
library(mice)
library(pROC)データ読み込み
データを読み込みます。
siero <-jaguchi("siero")## Loading required package: RCurl
## Loading required package: bitops
##
## Attaching package: 'RCurl'
##
## 以下のオブジェクトは 'package:mice' からマスクされています:
##
## completesiero2 <- t(siero)欠損値確認・データ前処理
今回の解析でNAになっているのは、
サンプル中の濃度が低すぎて測定器で拾うことができなかった値です。
はじめに欠損値の数から確認していくことにします。
md.pattern(siero2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## 261 1 1 1 1 1 1 1 1 1 1 1 1 0
## 4 1 1 1 1 1 1 1 1 1 1 0 1 1
## 2 1 1 1 1 1 0 1 1 1 1 1 1 1
## 2 1 1 1 0 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 0 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 0 1 1 1 1
## 3 1 1 1 1 1 1 1 0 1 1 1 1 1
## 1 1 1 0 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 0 1 1 1
## 1 1 0 1 1 1 1 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1 1 1 1 1 0 1
## 1 1 1 1 1 1 1 0 1 1 1 0 1 2
## 1 1 1 1 1 1 0 0 1 1 1 1 1 2
## 2 1 1 1 1 1 0 1 1 0 1 1 1 2
## 1 1 1 1 0 1 1 1 0 1 1 1 1 2
## 1 1 1 1 0 1 1 1 1 1 1 1 0 2
## 1 1 1 1 1 1 1 1 1 1 0 1 0 2
## 1 1 1 1 1 0 1 0 0 1 1 1 1 3
## 1 1 1 0 1 1 1 1 1 0 1 0 1 3
## 1 1 0 1 1 1 1 1 1 1 0 1 0 3
## 1 1 0 1 1 0 1 1 1 0 1 0 1 4
## 1 1 1 0 1 1 1 1 0 0 1 0 0 5
## 1 1 1 1 1 0 1 0 1 0 0 1 0 5
## 1 1 0 1 1 0 0 1 0 1 0 1 0 6
## 1 0 1 0 0 0 0 0 1 1 0 1 0 8
## 1 0 1 0 0 0 1 0 1 0 0 0 0 9
## 2 0 0 0 0 0 0 0 0 0 0 0 0 12
## 4 6 7 8 8 9 9 9 10 10 11 12 103すべての項目で欠損がないのは261件。 IDと食事の項目があるので実際は259物質です。
データ分割
次にデータを学習用、検証用に分割します。
set.seed(123)
samples=sort(sample(nrow(siero), nrow(siero)*0.5))
trainpre = siero[samples, -c(1:2)]
testpre = siero[-samples, -c(1:2)]
trainout = siero$Food[samples]
testout = siero$Food[-samples]
trainout2 = as.factor(trainout)
testout2 = as.factor(testout)
trainpre %>% dim## [1] 6 296testpre %>% dim## [1] 6 296乱数代入
続いて欠損値処理をします。
前回同様、NAにはガンマ分布から得た乱数を入れることにします。
ガンマ分布にしたのは物質の測定値なので、0はあっても負の値はないためです。
trainx <- data.matrix(trainpre, rownames.force = NA)
testx <- data.matrix(trainpre, rownames.force = NA)
n = 1000
rate = 5
set.seed(123)
trainx[is.na(trainx)]<- rgamma(n, rate, scale = 1/rate)
testx[is.na(testx)]<- rgamma(n, rate, scale = 1/rate)Enet学習用
Enetで学習用モデルを作ります。
#elastic net
set.seed(123)
con<-trainControl(method ="LOOCV",
preProc = c("center", "scale"),
summaryFunction = twoClassSummary,
classProbs = TRUE)
train_grid = expand.grid(alpha = 1:10 * 0.1, lambda = 10^{1:5 * -1})
enet_fit = train(trainx,trainout2,
method = "glmnet",
tuneGrid = train_grid,
trControl=con,
metric ="ROC")## Loading required package: glmnet
## Loading required package: Matrix
## Loading required package: foreach
## Loaded glmnet 2.0-2
##
##
## Attaching package: 'glmnet'
##
## 以下のオブジェクトは 'package:pROC' からマスクされています:
##
## aucenet_fit$bestTune## alpha lambda
## 5 0.1 0.1plot(varImp(enet_fit), top = 10)
Enet検証用
続いてこの学習用モデルを検証データに当てはめてみます。
pred <- predict(enet_fit, testx)
(conf.mat <- table(pred,testout2))## testout2
## pred After Before
## After 3 0
## Before 0 3どうやら上手く当てはまっているようですね。
ROC曲線
確認のためにROC曲線も書いておきましょう。
prob <- predict(enet_fit, testx, type = "prob")
rocCurve <- roc(response = testout2, predictor = prob$After)
plot(rocCurve, legacy.axes = TRUE)
##
## Call:
## roc.default(response = testout2, predictor = prob$After)
##
## Data: prob$After in 3 controls (testout2 After) > 3 cases (testout2 Before).
## Area under the curve: 1見事に直角ですね。 (2015/10/22 追記: 著者の福島さんからご連絡いただき、上記のコードを修正いたしました。)
ここまで、前処理からモデル構築をデータ分析プロセスに沿って解析してきました。
書籍にある基本的な流れについては触れられたかな、と思っていますが、
ここでは到底やりきれなかったたくさんの技が データ分析プロセス本編には載っています。
スタッキングなどの技術も次に書く論文では取り入れていければ…などと思いつつ、
一連の解析は一区切りにしようと思います。
Rで機械学習を使う方であれば、抑えておきたい本だと思います。