python, R, vimでデータマイニング

python, R, vim で疑問に思ったことなどを

いろいろ試したい時のデータセット:タイタニック

分析データを見つけるのは難しい

何かを試したい時に適したデータセットを見つけるのはなかなか難しい。
複数のデータセットを目的に応じて使い分けるのも骨が折れる。
ある程度汎用的に使えるデータセットを見つけたい。

求めるデータセット

私が使いやすいデータセット

  1. 2値分類、多項分類、回帰など汎用的
  2. 多すぎず少なすぎず適度な件数
  3. 馴染みのあるデータ

iris?

有名なiris。

  • 植物に興味がないので興味のあるデータではない。でも何度も使用したことがあるので馴染みはある。
  • データ件数が150件で少ない。
  • そのままでは2値分類に使用できない

求めるデータセット候補:タイタニック

kaggleで提供しているtitanic。

データ件数も変数の種類も豊富。映画も見ているので馴染みがある。
求めるデータセットに近い。

ただ、問題がひとつ。Rにデータセットが見つからない。
datasetsパッケージのTitanicはクロス集計に加工されたデータでローデータではない。

Titanic
## , , Age = Child, Survived = No
## 
##       Sex
## Class  Male Female
##   1st     0      0
##   2nd     0      0
##   3rd    35     17
##   Crew    0      0
## 
## , , Age = Adult, Survived = No
## 
##       Sex
## Class  Male Female
##   1st   118      4
##   2nd   154     13
##   3rd   387     89
##   Crew  670      3
## 
## , , Age = Child, Survived = Yes
## 
##       Sex
## Class  Male Female
##   1st     5      1
##   2nd    11     13
##   3rd    13     14
##   Crew    0      0
## 
## , , Age = Adult, Survived = Yes
## 
##       Sex
## Class  Male Female
##   1st    57    140
##   2nd    14     80
##   3rd    75     76
##   Crew  192     20

kaggleのタイタニックデータセット

githubにパッケージを用意してくれていた。 https://github.com/paulhendricks/titanic

パッケージのインストール

githubからパッケージをインストール。

### パッケージのインストール
## library(devtools)
## install_github("paulhendricks/titanic")
library(titanic)

データ内容 titanic_train

library(knitr)
kable(head(titanic_train))
PassengerIdSurvivedPclassNameSexAgeSibSpParchTicketFareCabinEmbarked
1 0 3 Braund, Mr. Owen Harris male 22 1 0 A/5 21171 7.2500   S
2 1 1 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female 38 1 0 PC 17599 71.2833 C85 C
3 1 3 Heikkinen, Miss. Laina female 26 0 0 STON/O2. 3101282 7.9250   S
4 1 1 Futrelle, Mrs. Jacques Heath (Lily May Peel) female 35 1 0 113803 53.1000 C123 S
5 0 3 Allen, Mr. William Henry male 35 0 0 373450 8.0500   S
6 0 3 Moran, Mr. James male NA 0 0 330877 8.4583   Q
library(Hmisc)
describe(titanic_train)
## titanic_train 
## 
##  12  Variables      891  Observations
## ---------------------------------------------------------------------------
## PassengerId 
##       n missing  unique    Info    Mean     .05     .10     .25     .50 
##     891       0     891       1     446    45.5    90.0   223.5   446.0 
##     .75     .90     .95 
##   668.5   802.0   846.5 
## 
## lowest :   1   2   3   4   5, highest: 887 888 889 890 891 
## ---------------------------------------------------------------------------
## Survived 
##       n missing  unique    Info     Sum    Mean 
##     891       0       2    0.71     342  0.3838 
## ---------------------------------------------------------------------------
## Pclass 
##       n missing  unique    Info    Mean 
##     891       0       3    0.81   2.309 
## 
## 1 (216, 24%), 2 (184, 21%), 3 (491, 55%) 
## ---------------------------------------------------------------------------
## Name 
##       n missing  unique 
##     891       0     891 
## 
## lowest : Abbing, Mr. Anthony                   Abbott, Mr. Rossmore Edward           Abbott, Mrs. Stanton (Rosa Hunt)      Abelson, Mr. Samuel                   Abelson, Mrs. Samuel (Hannah Wizosky)
## highest: de Mulder, Mr. Theodore               de Pelsmaeker, Mr. Alfons             del Carlo, Mr. Sebastiano             van Billiard, Mr. Austin Blyler       van Melkebeke, Mr. Philemon           
## ---------------------------------------------------------------------------
## Sex 
##       n missing  unique 
##     891       0       2 
## 
## female (314, 35%), male (577, 65%) 
## ---------------------------------------------------------------------------
## Age 
##       n missing  unique    Info    Mean     .05     .10     .25     .50 
##     714     177      88       1    29.7    4.00   14.00   20.12   28.00 
##     .75     .90     .95 
##   38.00   50.00   56.00 
## 
## lowest :  0.42  0.67  0.75  0.83  0.92
## highest: 70.00 70.50 71.00 74.00 80.00 
## ---------------------------------------------------------------------------
## SibSp 
##       n missing  unique    Info    Mean 
##     891       0       7    0.67   0.523 
## 
##             0   1  2  3  4 5 8
## Frequency 608 209 28 16 18 5 7
## %          68  23  3  2  2 1 1
## ---------------------------------------------------------------------------
## Parch 
##       n missing  unique    Info    Mean 
##     891       0       7    0.56  0.3816 
## 
##             0   1  2 3 4 5 6
## Frequency 678 118 80 5 4 5 1
## %          76  13  9 1 0 1 0
## ---------------------------------------------------------------------------
## Ticket 
##       n missing  unique 
##     891       0     681 
## 
## lowest : 110152      110413      110465      110564      110813     
## highest: W./C. 6608  W./C. 6609  W.E.P. 5734 W/C 14208   WE/P 5735   
## ---------------------------------------------------------------------------
## Fare 
##       n missing  unique    Info    Mean     .05     .10     .25     .50 
##     891       0     248       1    32.2   7.225   7.550   7.910  14.454 
##     .75     .90     .95 
##  31.000  77.958 112.079 
## 
## lowest :   0.000   4.013   5.000   6.237   6.438
## highest: 227.525 247.521 262.375 263.000 512.329 
## ---------------------------------------------------------------------------
## Cabin 
##       n missing  unique 
##     204     687     147 
## 
## lowest : A10 A14 A16 A19 A20, highest: F33 F38 F4  G6  T   
## ---------------------------------------------------------------------------
## Embarked 
##       n missing  unique 
##     889       2       3 
## 
## C (168, 19%), Q (77, 9%), S (644, 72%) 
## ---------------------------------------------------------------------------

データ内容 titanic_test

kable(head(titanic_test))
PassengerIdPclassNameSexAgeSibSpParchTicketFareCabinEmbarked
892 3 Kelly, Mr. James male 34.5 0 0 330911 7.8292   Q
893 3 Wilkes, Mrs. James (Ellen Needs) female 47.0 1 0 363272 7.0000   S
894 2 Myles, Mr. Thomas Francis male 62.0 0 0 240276 9.6875   Q
895 3 Wirz, Mr. Albert male 27.0 0 0 315154 8.6625   S
896 3 Hirvonen, Mrs. Alexander (Helga E Lindqvist) female 22.0 1 1 3101298 12.2875   S
897 3 Svensson, Mr. Johan Cervin male 14.0 0 0 7538 9.2250   S
describe(titanic_test)
## titanic_test 
## 
##  11  Variables      418  Observations
## ---------------------------------------------------------------------------
## PassengerId 
##       n missing  unique    Info    Mean     .05     .10     .25     .50 
##     418       0     418       1    1100   912.9   933.7   996.2  1100.5 
##     .75     .90     .95 
##  1204.8  1267.3  1288.2 
## 
## lowest :  892  893  894  895  896, highest: 1305 1306 1307 1308 1309 
## ---------------------------------------------------------------------------
## Pclass 
##       n missing  unique    Info    Mean 
##     418       0       3    0.83   2.266 
## 
## 1 (107, 26%), 2 (93, 22%), 3 (218, 52%) 
## ---------------------------------------------------------------------------
## Name 
##       n missing  unique 
##     418       0     418 
## 
## lowest : Abbott, Master. Eugene Joseph                 Abelseth, Miss. Karen Marie                   Abelseth, Mr. Olaus Jorgensen                 Abrahamsson, Mr. Abraham August Johannes      Abrahim, Mrs. Joseph (Sophie Halaut Easu)    
## highest: de Brito, Mr. Jose Joaquim                    de Messemaeker, Mr. Guillaume Joseph          del Carlo, Mrs. Sebastiano (Argenia Genovesi) van Billiard, Master. James William           van Billiard, Master. Walter John             
## ---------------------------------------------------------------------------
## Sex 
##       n missing  unique 
##     418       0       2 
## 
## female (152, 36%), male (266, 64%) 
## ---------------------------------------------------------------------------
## Age 
##       n missing  unique    Info    Mean     .05     .10     .25     .50 
##     332      86      79       1   30.27     8.0    16.1    21.0    27.0 
##     .75     .90     .95 
##    39.0    50.0    57.0 
## 
## lowest :  0.17  0.33  0.75  0.83  0.92
## highest: 62.00 63.00 64.00 67.00 76.00 
## ---------------------------------------------------------------------------
## SibSp 
##       n missing  unique    Info    Mean 
##     418       0       7    0.67  0.4474 
## 
##             0   1  2 3 4 5 8
## Frequency 283 110 14 4 4 1 2
## %          68  26  3 1 1 0 0
## ---------------------------------------------------------------------------
## Parch 
##       n missing  unique    Info    Mean 
##     418       0       8    0.53  0.3923 
## 
##             0  1  2 3 4 5 6 9
## Frequency 324 52 33 3 2 1 1 2
## %          78 12  8 1 0 0 0 0
## ---------------------------------------------------------------------------
## Ticket 
##       n missing  unique 
##     418       0     363 
## 
## lowest : 110469      110489      110813      111163      112051     
## highest: W./C. 14260 W./C. 14266 W./C. 6607  W./C. 6608  W.E.P. 5734 
## ---------------------------------------------------------------------------
## Fare 
##       n missing  unique    Info    Mean     .05     .10     .25     .50 
##     417       1     169       1   35.63   7.229   7.642   7.896  14.454 
##     .75     .90     .95 
##  31.500  79.200 151.550 
## 
## lowest :   0.000   3.171   6.438   6.496   6.950
## highest: 227.525 247.521 262.375 263.000 512.329 
## ---------------------------------------------------------------------------
## Cabin 
##       n missing  unique 
##      91     327      76 
## 
## lowest : A11   A18   A21   A29   A34  
## highest: F G63 F2    F33   F4    G6    
## ---------------------------------------------------------------------------
## Embarked 
##       n missing  unique 
##     418       0       3 
## 
## C (102, 24%), Q (46, 11%), S (270, 65%) 
## ---------------------------------------------------------------------------

dplyr:group_by に嵌る

問題

dplyrのgroup_byを使用してグループ別の集計処理をしたかったのだが
グループ処理が働かずに嵌った。

問題例

library(dplyr)
library(plyr)

mtcars %>% group_by(vs, am) %>% summarise(max = max(mpg))
##    max
## 1 33.9

グループ別の最大値を求めたいのに一つの値しか得られない。

検索 “dplyr group_by not working”

検索ワード「dplyr group_by not working」で検索して情報を得た。

http://stackoverflow.com/questions/26923862/why-are-my-dplyr-group-by-summarize-not-working-properly-name-collision-with

I believe you’ve loaded plyr after dplyr, which is why you are getting an overall summary instead of a grouped summary.

dplyrの後にplyr をロードするとだめらしい。

問題解決例

library(dplyr)
library(plyr)
detach(package:plyr)

mtcars %>% group_by(vs, am) %>% summarise(max = max(mpg))
## Source: local data frame [4 x 3]
## Groups: vs
## 
##   vs am  max
## 1  0  0 19.2
## 2  0  1 26.0
## 3  1  0 24.4
## 4  1  1 33.9

xgboostを試してみる

xgboostで精度の高いモデルを構築できるらしい。
それから、fevalでモデル精度指標を指定できるところも良さそう。
xgboostを試してみたい。

まずは試してみる。

caretを使用すると色々なモデルが共通の文法で書けるので便利。

library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
model1 <- train(
  Species ~ .,
  iris,
  method="xgbTree"
)
## Loading required package: xgboost
## Loading required package: plyr
model1
## eXtreme Gradient Boosting 
## 
## 150 samples
##   4 predictor
##   3 classes: 'setosa', 'versicolor', 'virginica' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 150, 150, 150, 150, 150, 150, ... 
## Resampling results across tuning parameters:
## 
##   max_depth  nrounds  Accuracy   Kappa      Accuracy SD  Kappa SD  
##   1           50      0.9485004  0.9220085  0.02628702   0.03960171
##   1          100      0.9463580  0.9187624  0.02515099   0.03790241
##   1          150      0.9455416  0.9174981  0.02573184   0.03884031
##   2           50      0.9468626  0.9194865  0.02621117   0.03952144
##   2          100      0.9475898  0.9206234  0.02740915   0.04130599
##   2          150      0.9475898  0.9206234  0.02740915   0.04130599
##   3           50      0.9506307  0.9251477  0.02432345   0.03674940
##   3          100      0.9513580  0.9262899  0.02443380   0.03685994
##   3          150      0.9512716  0.9261698  0.02402944   0.03625242
## 
## Tuning parameter 'eta' was held constant at a value of 0.3
## Accuracy was used to select the optimal model using  the largest value.
## The final values used for the model were nrounds = 100, max_depth = 3
##  and eta = 0.3.
summary(model1)
##             Length Class              Mode       
## handle          1  xgb.Booster.handle externalptr
## raw         90996  -none-             raw        
## xNames          4  -none-             character  
## problemType     1  -none-             character  
## tuneValue       3  data.frame         list       
## obsLevels       3  -none-             character

思ったよりも簡単にできた。

feval

fevalで独自の評価関数を指定したい。

マニュアルを確認する。
feval custimized evaluation function.
Returns list(metric=‘metric-name’, value=‘metric-value’)
with given prediction and dtrain.

caret:train で試してみる。
評価関数はrmspeを自作する。

rmspe <- function(prediction, dtrain) {
  .observed <- getinfo(dtrain, "label") 
  list(
    metric = "RMSPE",
      value = sqrt(mean(((.observed- prediction) / .observed) ^ 2))
  )
}

model2 <- train(
  Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
  iris,
  method = "xgbTree",
  feval = rmspe
)
model2
## eXtreme Gradient Boosting 
## 
## 150 samples
##   4 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 150, 150, 150, 150, 150, 150, ... 
## Resampling results across tuning parameters:
## 
##   max_depth  nrounds  RMSE       Rsquared   RMSE SD     Rsquared SD
##   1           50      0.3583983  0.8086825  0.02845745  0.04300446 
##   1          100      0.3592624  0.8085968  0.02550611  0.04302866 
##   1          150      0.3633599  0.8053671  0.02583143  0.04393736 
##   2           50      0.3638931  0.8022369  0.03145324  0.04766807 
##   2          100      0.3780569  0.7881979  0.03225296  0.05362911 
##   2          150      0.3892047  0.7772812  0.03294389  0.05639148 
##   3           50      0.3745153  0.7923352  0.02861350  0.04923131 
##   3          100      0.3925762  0.7754515  0.03001914  0.05058433 
##   3          150      0.4014436  0.7671566  0.03145388  0.05187344 
## 
## Tuning parameter 'eta' was held constant at a value of 0.3
## RMSE was used to select the optimal model using  the smallest value.
## The final values used for the model were nrounds = 50, max_depth = 1
##  and eta = 0.3.

だめ。RMSEが評価指標として採用されている。 私の知識ではcaret:trainでfevalを指定できない。

feval xgboost:xgb.train でfevalを指定

caretを使用しないと新しい文法をマニュアルなどで調べないといけないので嫌だが、とにかく試してみる。

watchlist

目的変数やクロスバリデーション用のデータはwatchlistに指定するらしい。
また、xgb.DMatrix 形式にする必要があるらしい?

.feature <- data.matrix(iris[, 2:4])
.watchlist <- list(
  val = xgb.DMatrix(data = .feature, label = iris$Sepal.Length),
  train = xgb.DMatrix(data = .feature, label = iris$Sepal.Length)
)

xgb.train で試してみる。

model.3 <- xgb.train(
  data = xgb.DMatrix(data = .feature, label = iris$Sepal.Length),
  nround = 10,
  watchlist = .watchlist,
  feval=rmspe
)
## [0]  val-RMSPE:0.641076263296985 train-RMSPE:0.641076263296985
## [1]  val-RMSPE:0.450932635857561 train-RMSPE:0.450932635857561
## [2]  val-RMSPE:0.318373538117156 train-RMSPE:0.318373538117156
## [3]  val-RMSPE:0.227032609122525 train-RMSPE:0.227032609122525
## [4]  val-RMSPE:0.164214553125964 train-RMSPE:0.164214553125964
## [5]  val-RMSPE:0.120398490746834 train-RMSPE:0.120398490746834
## [6]  val-RMSPE:0.0907635245058501    train-RMSPE:0.0907635245058501
## [7]  val-RMSPE:0.0706479694183348    train-RMSPE:0.0706479694183348
## [8]  val-RMSPE:0.0571393392305592    train-RMSPE:0.0571393392305592
## [9]  val-RMSPE:0.0484302684656021    train-RMSPE:0.0484302684656021
summary(model.3)
##        Length Class              Mode       
## handle    1   xgb.Booster.handle externalptr
## raw    7312   -none-             raw

なんかできた。

dplyr, mutateを用いたデータハンドリング, data wrangling

モデル構築前の前処理

モデル構築する前には前処理が必要。
むしろ、前処理の方がが時間がかかる。
モデル構築ではトレーニング用とテスト用に分割するので
2つのデータセットに同じ前処理を施す。

この前処理群の管理や新変数名の設定方法など
煩雑になりがちで困っていたが自分なりの実践方法を整理しておきたい。

mutate_eachの気になる点

以前のブログで触れたがmutate_eachの気になる点として、
適用する関数が一つの場合に元変数を上書きしてしまう。

一つの関数で元変数が上書きされる例

library(dplyr)  
## 
## Attaching package: 'dplyr'
## 
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag 
## 
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union
iris <- tbl_df(iris)  
mutate_each(iris, funs(log), -Species)  
## Source: local data frame [150 x 5]
## 
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           (dbl)       (dbl)        (dbl)       (dbl)  (fctr)
## 1      1.629241    1.252763    0.3364722  -1.6094379  setosa
## 2      1.589235    1.098612    0.3364722  -1.6094379  setosa
## 3      1.547563    1.163151    0.2623643  -1.6094379  setosa
## 4      1.526056    1.131402    0.4054651  -1.6094379  setosa
## 5      1.609438    1.280934    0.3364722  -1.6094379  setosa
## 6      1.686399    1.360977    0.5306283  -0.9162907  setosa
## 7      1.526056    1.223775    0.3364722  -1.2039728  setosa
## 8      1.609438    1.223775    0.4054651  -1.6094379  setosa
## 9      1.481605    1.064711    0.3364722  -1.6094379  setosa
## 10     1.589235    1.131402    0.4054651  -2.3025851  setosa
## ..          ...         ...          ...         ...     ...

解決策

あまり綺麗な方法ではないがeval関数を使用することで解決する。

mutate_each(iris, funs(eval, log), -Species)  
## Source: local data frame [150 x 13]
## 
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           (dbl)       (dbl)        (dbl)       (dbl)  (fctr)
## 1           5.1         3.5          1.4         0.2  setosa
## 2           4.9         3.0          1.4         0.2  setosa
## 3           4.7         3.2          1.3         0.2  setosa
## 4           4.6         3.1          1.5         0.2  setosa
## 5           5.0         3.6          1.4         0.2  setosa
## 6           5.4         3.9          1.7         0.4  setosa
## 7           4.6         3.4          1.4         0.3  setosa
## 8           5.0         3.4          1.5         0.2  setosa
## 9           4.4         2.9          1.4         0.2  setosa
## 10          4.9         3.1          1.5         0.1  setosa
## ..          ...         ...          ...         ...     ...
## Variables not shown: Sepal.Length_eval (dbl), Sepal.Width_eval (dbl),
##   Petal.Length_eval (dbl), Petal.Width_eval (dbl), Sepal.Length_log (dbl),
##   Sepal.Width_log (dbl), Petal.Length_log (dbl), Petal.Width_log (dbl)

余分に作成される変数 *_eval はあとで削除する。

mutate_each(iris, funs(eval, log), -Species) %>%  
  select(-ends_with("_eval"))  
## Source: local data frame [150 x 9]
## 
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           (dbl)       (dbl)        (dbl)       (dbl)  (fctr)
## 1           5.1         3.5          1.4         0.2  setosa
## 2           4.9         3.0          1.4         0.2  setosa
## 3           4.7         3.2          1.3         0.2  setosa
## 4           4.6         3.1          1.5         0.2  setosa
## 5           5.0         3.6          1.4         0.2  setosa
## 6           5.4         3.9          1.7         0.4  setosa
## 7           4.6         3.4          1.4         0.3  setosa
## 8           5.0         3.4          1.5         0.2  setosa
## 9           4.4         2.9          1.4         0.2  setosa
## 10          4.9         3.1          1.5         0.1  setosa
## ..          ...         ...          ...         ...     ...
## Variables not shown: Sepal.Length_log (dbl), Sepal.Width_log (dbl),
##   Petal.Length_log (dbl), Petal.Width_log (dbl)

前処理の関数化

前処理を関数化することで異なるデータセットに処理を適用するようにする。

  1. Species を文字列に型変換
  2. 数値型の4変数を対数変換
  3. 数値型の値が5以上の場合は5に丸める。
.wrangling <- function(.data) {  
  .data %>%  
    mutate(Species_as.character = as.character(Species)) %>%  
    mutate_each(  
      funs(  
        eval,   
        log,  
        round5 = ifelse(. > 5, 5, .)  
      ),  
      -starts_with("Species")  
    ) %>%  
    select(-ends_with("_eval"))  
}  
iris.wrangling <- .wrangling(iris)  
summary(iris.wrangling)  
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species   Species_as.character Sepal.Length_log Sepal.Width_log 
##  setosa    :50   Length:150           Min.   :1.459    Min.   :0.6931  
##  versicolor:50   Class :character     1st Qu.:1.629    1st Qu.:1.0296  
##  virginica :50   Mode  :character     Median :1.758    Median :1.0986  
##                                       Mean   :1.755    Mean   :1.1074  
##                                       3rd Qu.:1.856    3rd Qu.:1.1939  
##                                       Max.   :2.067    Max.   :1.4816  
##  Petal.Length_log Petal.Width_log   Sepal.Length_round5 Sepal.Width_round5
##  Min.   :0.000    Min.   :-2.3026   Min.   :4.300       Min.   :2.000     
##  1st Qu.:0.470    1st Qu.:-1.2040   1st Qu.:5.000       1st Qu.:2.800     
##  Median :1.470    Median : 0.2624   Median :5.000       Median :3.000     
##  Mean   :1.175    Mean   :-0.1723   Mean   :4.955       Mean   :3.057     
##  3rd Qu.:1.629    3rd Qu.: 0.5878   3rd Qu.:5.000       3rd Qu.:3.300     
##  Max.   :1.932    Max.   : 0.9163   Max.   :5.000       Max.   :4.400     
##  Petal.Length_round5 Petal.Width_round5
##  Min.   :1.000       Min.   :0.100     
##  1st Qu.:1.600       1st Qu.:0.300     
##  Median :4.350       Median :1.300     
##  Mean   :3.565       Mean   :1.199     
##  3rd Qu.:5.000       3rd Qu.:1.800     
##  Max.   :5.000       Max.   :2.500

caret:trainに慣れる

パッケージ caret

The caret Package

色々なアルゴリズムを個別のパッケージで対応してきた。
でも、それぞれの使い方を調べながら対応するのが面倒。
caretは多くのアルゴリズムを一つのパッケージにまとめてくれている。
また、モデル構築で必要なツールを提供してくれている。

caret:train

caret:trainでモデルを構築する。

オプション:method

method = でアルゴリズムを指定。多くのアルゴリズムを指定できる。

指定できるアルゴリズム

今回は“lm”, “rf”, “gbm”を使用する。

  • “lm”は一般的な回帰分析の例として。
  • “rf”はよく利用していたランダムフォレスト。
  • gbm”は今後利用したいアルゴリズム。精度が高くランダムフォレストより高速に処理できる。

シンプルなモデル構築

library(caret)
library(dplyr)

iris <- tbl_df(iris)

model.lm <- train(
    data = iris,
    Sepal.Length ~ .,
    method = "lm"
)

model.rf <- train(
    data = iris,
    Sepal.Length ~ .,
    method = "rf"
)

model.gbm <- train(
    data = iris,
    Sepal.Length ~ .,
    method = "gbm"
)
model.lm
## Linear Regression 
## 
## 150 samples
##   4 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 150, 150, 150, 150, 150, 150, ... 
## Resampling results
## 
##   RMSE       Rsquared   RMSE SD     Rsquared SD
##   0.3248528  0.8457439  0.01828276  0.02130072 
## 
## 
summary(model.lm)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.79424 -0.21874  0.00899  0.20255  0.73103 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.17127    0.27979   7.760 1.43e-12 ***
## Sepal.Width        0.49589    0.08607   5.761 4.87e-08 ***
## Petal.Length       0.82924    0.06853  12.101  < 2e-16 ***
## Petal.Width       -0.31516    0.15120  -2.084  0.03889 *  
## Speciesversicolor -0.72356    0.24017  -3.013  0.00306 ** 
## Speciesvirginica  -1.02350    0.33373  -3.067  0.00258 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3068 on 144 degrees of freedom
## Multiple R-squared:  0.8673, Adjusted R-squared:  0.8627 
## F-statistic: 188.3 on 5 and 144 DF,  p-value: < 2.2e-16
model.rf
## Random Forest 
## 
## 150 samples
##   4 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 150, 150, 150, 150, 150, 150, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE       Rsquared   RMSE SD     Rsquared SD
##   2     0.3749241  0.8155266  0.03152655  0.03701402 
##   3     0.3629939  0.8235616  0.02747483  0.03561889 
##   5     0.3627847  0.8207195  0.02842011  0.03995335 
## 
## RMSE was used to select the optimal model using  the smallest value.
## The final value used for the model was mtry = 5.
summary(model.rf)
##                 Length Class      Mode     
## call              4    -none-     call     
## type              1    -none-     character
## predicted       150    -none-     numeric  
## mse             500    -none-     numeric  
## rsq             500    -none-     numeric  
## oob.times       150    -none-     numeric  
## importance        5    -none-     numeric  
## importanceSD      0    -none-     NULL     
## localImportance   0    -none-     NULL     
## proximity         0    -none-     NULL     
## ntree             1    -none-     numeric  
## mtry              1    -none-     numeric  
## forest           11    -none-     list     
## coefs             0    -none-     NULL     
## y               150    -none-     numeric  
## test              0    -none-     NULL     
## inbag             0    -none-     NULL     
## xNames            5    -none-     character
## problemType       1    -none-     character
## tuneValue         1    data.frame list     
## obsLevels         1    -none-     logical
model.gbm
## Stochastic Gradient Boosting 
## 
## 150 samples
##   4 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 150, 150, 150, 150, 150, 150, ... 
## Resampling results across tuning parameters:
## 
##   interaction.depth  n.trees  RMSE       Rsquared   RMSE SD   
##   1                   50      0.3913974  0.7822964  0.03051559
##   1                  100      0.3737524  0.7994068  0.02906737
##   1                  150      0.3706075  0.8027339  0.02712242
##   2                   50      0.3741219  0.7996173  0.02745996
##   2                  100      0.3688091  0.8048558  0.02672018
##   2                  150      0.3675264  0.8055333  0.02610464
##   3                   50      0.3712947  0.8009253  0.02792836
##   3                  100      0.3684993  0.8043281  0.02672187
##   3                  150      0.3710270  0.8023887  0.02685090
##   Rsquared SD
##   0.03417764 
##   0.03301104 
##   0.03097777 
##   0.02871627 
##   0.02685198 
##   0.02633916 
##   0.02546331 
##   0.02405717 
##   0.02456531 
## 
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
## 
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## RMSE was used to select the optimal model using  the smallest value.
## The final values used for the model were n.trees = 150,
##  interaction.depth = 2, shrinkage = 0.1 and n.minobsinnode = 10.
summary(model.gbm)

##                                 var    rel.inf
## Petal.Length           Petal.Length 79.9216869
## Sepal.Width             Sepal.Width  9.0176394
## Petal.Width             Petal.Width  5.4888793
## Speciesversicolor Speciesversicolor  5.3919231
## Speciesvirginica   Speciesvirginica  0.1798713

rfとgbmはパラメータチューニングのため複数モデルを構築し性能が良いモデルのパラメータを採用している。
rfはmtryを、gbmはinteraction.depthとn.treesをチューニングしている。
モデルの評価はRMSEを採用。
rfはmtry=5のモデル、gbmはinteraction.depth = 2, n.trees = 150 のモデルを採用。
最終結果として

  1. lm : RMSE = 0.324
  2. rf : RMSE = 0.362
  3. gbm : RMSE = 0.367

lmのモデルが最も性能が良い。

後でパラーメータチューニング、リサンプリングの指定方法を整理したいが今日はここまで。

欠損値置換や因子ベクトル化をmutate_eachで

 Rでのデータ加工がいつも捗らなくて困っていた。
dplyrをより活用することでもっと効率的にしたい。
そのなかでもmutate_eachは使う場面が多そうなので慣れておきたい。

全変数を因子ベクトル化

 例えば、カテゴリ型の変数名を指名して対象変数を因子ベクトルにしたいことは多い。
自作のコードが煩雑になり困っていた。
そこでdplyr:mutate_eachを使用してみる。
irisの全変数をfactorで因子ベクトル化する。

それから関数はmutate_eachとmutate_each_がある。 (違いは他を当たってください)

library(dplyr)
## 
## Attaching package: 'dplyr'
## 
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag 
## 
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union
packageVersion("dplyr")
## [1] '0.4.3'
iris <- tbl_df(iris)
iris.factor <- mutate_each(iris, funs(factor))
str(iris.factor)
## Classes 'tbl_df', 'tbl' and 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: Factor w/ 35 levels "4.3","4.4","4.5",..: 9 7 5 4 8 12 4 8 2 7 ...
##  $ Sepal.Width : Factor w/ 23 levels "2","2.2","2.3",..: 15 10 12 11 16 19 14 14 9 11 ...
##  $ Petal.Length: Factor w/ 43 levels "1","1.1","1.2",..: 5 5 4 6 5 8 5 6 5 6 ...
##  $ Petal.Width : Factor w/ 22 levels "0.1","0.2","0.3",..: 2 2 2 2 2 4 3 2 2 1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

変数を指定する。

 変数名はlistなどで指定できる。その他にも便利な指定方法があるようなので他を当たってみてください。

iris.factor2 <- mutate_each_(iris, funs(factor), list("Sepal.Length", "Species"))
str(iris.factor2)
## Classes 'tbl_df', 'tbl' and 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: Factor w/ 35 levels "4.3","4.4","4.5",..: 9 7 5 4 8 12 4 8 2 7 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
iris.factor2 <- mutate_each(iris, funs(factor), Species)
iris.factor2 <- mutate_each_(iris, funs(factor), list(quote(Species), quote(Sepal.Length)))
iris.factor2 <- mutate_each(iris, funs(factor), -Species)
iris.factor2 <- mutate_each_(iris, funs(factor), "-Sepal.Length")

新規追加する変数名

 一つ困っていることがある。
funsで関数名を指定する。指定した関数がひとつの場合は元の変数を上書きする。
指定した関数が2つ以上の場合は元の変数を上書きせず、関数名をsuffixとした新規変数を作成する。
適用したい関数がひとつの場合でも元の変数を残したいのだがその方法がわからない。

str(mutate_each_(iris, funs(half = ./2), list("Sepal.Length", "Sepal.Width")))
## Classes 'tbl_df', 'tbl' and 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  2.55 2.45 2.35 2.3 2.5 2.7 2.3 2.5 2.2 2.45 ...
##  $ Sepal.Width : num  1.75 1.5 1.6 1.55 1.8 1.95 1.7 1.7 1.45 1.55 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
str(mutate_each_(iris, funs(half = ./2, log), list("Sepal.Length", "Sepal.Width")))
## Classes 'tbl_df', 'tbl' and 'data.frame':    150 obs. of  9 variables:
##  $ Sepal.Length     : num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width      : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length     : num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width      : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species          : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Sepal.Length_half: num  2.55 2.45 2.35 2.3 2.5 2.7 2.3 2.5 2.2 2.45 ...
##  $ Sepal.Width_half : num  1.75 1.5 1.6 1.55 1.8 1.95 1.7 1.7 1.45 1.55 ...
##  $ Sepal.Length_log : num  1.63 1.59 1.55 1.53 1.61 ...
##  $ Sepal.Width_log  : num  1.25 1.1 1.16 1.13 1.28 ...

経過月数の計算 library(mondate)

 KAGGLEの分析で経過月数を求めたいと思った。

過去に書き散らかして放置していたブログでも同様の記事を投稿したことがあった。月日が立ち、そのときに利用していたライブラリ名すら忘れていた。

 「r months between」で検索してみるとStack Overflow に質問を投げている人がいる。でも、回答は自作関数ばかりで該当するパッケージのヒントは得られない。

 忘れていた過去ブログを探し出してパッケージ名を確認した。

「mondate」

 このブログは過去に放置した記事など整理も含めて、Rを使用して疑問に思ったことをまとめて記録するブログにしていきたい。

 

使用例 

mondate.ymd

mondateクラスのオブジェクトを作成するにはmondate.ymd関数を使用する。
引数dを省略するとその月の最後の日で定義する。

if (!require("mondate")) {
  install.packages("mondate")
  library("mondate")
}
## Loading required package: mondate
## 
## Attaching package: 'mondate'
## 
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      as.difftime
from <- mondate.ymd(2015, 10)
to <- mondate.ymd(2016, 1)
MonthsBetween(from, to)
## [1] 3
## attr(,"timeunits")
## [1] "months"
from2 <- as.mondate(Sys.Date())
MonthsBetween(from2, to)
## [1] 3.516129
## attr(,"timeunits")
## [1] "months"
from3 <- mondate.ymd(
  year(Sys.Date()),
  month(Sys.Date())
)

MonthsBetween(from3, to)
## [1] 3
## attr(,"timeunits")
## [1] "months"