以下の発表で紹介したデータについて解析したフローを示します。
Tokyo.Rでの発表スライドはこちら 。

今回は前回の解析に、ご恵贈頂いたデータ分析プロセス3章に記載されている
前処理のフローを加えた内容となります（変数代入はやりきれませんでした…すいません）。

はじめにdevtools::install_github(“dichika/jaguchi”)で
jaguchiパッケージをインストールしておきます。

パッケージインストール・読み込み

library(jaguchi)
library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
library(mice)
## Loading required package: Rcpp
## mice 2.22 2014-06-10

データ読み込み

まずデータを読み込みます。

siero <-jaguchi("siero")
## Loading required package: RCurl
## Loading required package: bitops
## 
## Attaching package: 'RCurl'
## 
## The following object is masked from 'package:mice':
## 
##     complete
siero2 <- t(siero)

欠損値確認・データ前処理

今回の解析でNAになっているのは、
サンプル中の濃度が低すぎて測定器で拾うことができなかった値です。

この場合だと欠損値の有無が欠損値を持つ変数自身と関係しているので、
Missing Not At Random (MNAR) と考えていいんでしょうか。
機器分析では伝統的に機器が検出可能な値の1/2の値を入れたり、
検出可能な値の1/2を平均にした正規分布の乱数を入れたりしてますが…。

ともあれ欠損値の数から確認していくことにします。

siero <-jaguchi("siero")
siero2 <- t(siero)
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物質ということになります。

続けて欠損値処理をしてみます。 前回ではとりあえずNA取り除いたらなんとかなるだろ、
ということで、前処理を行いましたが、
この手法はリストワイズ法と言う手法のようです。

リストワイズ法

y <- siero$Food
y1 <- as.factor(y)

x <- data.matrix(siero[,-c(1,2)], rownames.force  = NA)
xbase <- t(x)
xlist <- na.omit(xbase)
xlist2 <- t(xlist)

比較として多重代入法を試みたのですが、
デフォルトでは上手く動きませんでした。
ランダムフォレストを使うなど、幾つか試しては見たのですが、
変数の数が多いためか猛烈に計算時間がかかってしまったので、
ちょっとこの部分については妥協しました、すいません。
もうちょっと向いてそうなデータで再挑戦します。

今回はNAにガンマ分布から得た乱数を入れることにします。 ガンマ分布にしたのは、測定値なので負の値がないためです。

乱数代入

x <- data.matrix(siero[,-c(1,2)], rownames.force  = NA)
n = 1000
shape = 1000
rate = 5
set.seed(123)
x[is.na(x)]<- rgamma(n, shape, rate, scale = 1/rate)

続いて前回の結果と今回の結果を比較してみましょう。

Enetリストワイズ

#elastic net
set.seed(123)
con<-trainControl(method ="LOOCV",
                  preProc = c("center", "scale"))

train_grid = expand.grid(alpha = 1:10 * 0.1, lambda = 10^{1:5 * -1})

enet_fit = train(xlist2,y1,  
                 method = "glmnet", 
                 tuneGrid = train_grid,
                 trControl=con)
## Loading required package: glmnet
## Loading required package: Matrix
## Loading required package: foreach
## Loaded glmnet 2.0-2
enet_fit$bestTune
##   alpha lambda
## 5   0.1    0.1
plot(varImp(enet_fit), top = 10)

Enet乱数代入

#elastic net
set.seed(123)
con<-trainControl(method ="LOOCV",
                  preProc = c("center", "scale"))

train_grid = expand.grid(alpha = 1:10 * 0.1, lambda = 10^{1:5 * -1})

enet_fit2 = train(x,y1,  
                 method = "glmnet", 
                 tuneGrid = train_grid,
                 trControl=con)

enet_fit2$bestTune
##   alpha lambda
## 5   0.1    0.1
plot(varImp(enet_fit2), top = 10)

少なくとも変数重要度から選ばれる物質が変わってくることが確認できました。
次回は本の4章を参考に、予測モデルの構築について詳細を詰めていきたいと思います。