niszetの日記

10年目エンジニアが最近勉強したことを忘れないようにメモをする

(R) 隠れtidyverse modelr編 とりあえず使ってみる 1回目

久々のRの日記な気がする…

これから数回にかけて、modelrというパッケージに入っている関数を使ってみる、という趣旨で日記を書いていきます。

知らないうちに、インストールされてる…

modelrパッケージについて、日本語で書かれた情報が全然見当たりませんが、実はtidyverseをインストールすると一緒に入っています。ただし、library(tidyverse)ではロードされないので、意識して読み込む必要があります。

# インストールされてなければ
# install.packages("tidyverse")

# tidyverseをロード
library(tidyverse)

# modelrをロード
library(modelr)

# broomをロード
library(broom)

確認してみましょう。

 tidyverse::tidyverse_packages()
#  [1] "broom"     "dplyr"     "forcats"   "ggplot2"   "haven"     "httr"      "hms"       "jsonlite"  "lubridate"
# [10] "magrittr"  "modelr"    "purrr"     "readr"     "readxl"    "stringr"   "tibble"    "rvest"     "tidyr"    
# [19] "xml2"      "tidyverse"

お。11番目にありますね。

なお、broomとmodelrそれぞれにbootstrapという関数が含まれているため、このようなメッセージが出てきます。

 次のパッケージを付け加えます: ‘broom’ 

 以下のオブジェクトは ‘package:modelr’ からマスクされています: 

     bootstrap 

このbootstrapのそれぞれのパッケージでの差についてはいずれ…(まだちゃんと見てない)とりあえず気にせずに進めましょう。

broomと仲が良い。

以前、broom::tidy、broom::tidyMCMCについて書きました。

niszet.hatenablog.com

このときはstanの結果についてみてみましたが、tidyはlmやglmの結果をtidyにしてくれます(実際はtidyでtidyMCMCと同じようにstanの結果も見れる)

modelrはlmやglmなどの解析をtidyにしてくれます。また、ブートストラップやらCVなどに使える関数も色々入っており、broomとあわせることで解析が楽になるのではないかと思います。

実際に使ってみる

とりあえず、実際に使って結果を見る方が早いので。 modelrの例文だとmtcarsを使っていることが多いのですが、簡単のためにirisのデータを使ってみます。

plot(iris$Sepal.Length ~ iris$Petal.Length)

こんな図になりますね。特に面白みのない図です。 f:id:niszet:20170823225820p:plain

さて、lmに入れてみます。

lm(iris$Sepal.Length ~ iris$Petal.Length)
# 
# Call:
# lm(formula = iris$Sepal.Length ~ iris$Petal.Length)
# 
# Coefficients:
#       (Intercept)  iris$Petal.Length  
#            4.3066             0.4089  

せっかくなのでbroom::tidyで整形します。

broom::tidy(lm(iris$Sepal.Length ~ iris$Petal.Length))
#                term  estimate  std.error statistic       p.value
# 1       (Intercept) 4.3066034 0.07838896  54.93890 2.426713e-100
# 2 iris$Petal.Length 0.4089223 0.01889134  21.64602  1.038667e-47

printで見れないパラメータはsummaryで見れる値です。std.error、statistic、p.valueですね。

summary(lm(iris$Sepal.Length ~ iris$Petal.Length))
# 
# Call:
# lm(formula = iris$Sepal.Length ~ iris$Petal.Length)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.24675 -0.29657 -0.01515  0.27676  1.00269 
# 
# Coefficients:
#                   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)        4.30660    0.07839   54.94   <2e-16 ***
# iris$Petal.Length  0.40892    0.01889   21.65   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4071 on 148 degrees of freedom
# Multiple R-squared:   0.76,  Adjusted R-squared:  0.7583 
# F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

summaryから、statistic = t value、std.error = Std. Error、p.value = Pr(>|t|)が対応していますね。tidyMCMCと違って、オプションは特にないのだと思いますが、見つけたら教えて下さい…。

add_predictions

broomの話しかしてなかった。add_precictionsは、データとモデルの式を与えると、もとのデータ列にモデルでの予測値を追加したものを返します。

head(modelr::add_predictions(iris, lm(iris$Sepal.Length ~ iris$Petal.Length)),n=3)
#   Sepal.Length Sepal.Width Petal.Length Petal.Width Species     pred
# 1          5.1         3.5          1.4         0.2  setosa 4.879095
# 2          4.9         3.0          1.4         0.2  setosa 4.879095
# 3          4.7         3.2          1.3         0.2  setosa 4.838202

先の結果から、Sepal.Length = 4.3066034 + 0.4089223 * Petal.Length が得られていますので、1行目の1.4を与えて計算すると4.879095が得られ、一致しますね。

add_residuals

今度は残差を足す方です。

head(modelr::add_residuals(iris, lm(iris$Sepal.Length ~ iris$Petal.Length)),n=3)
#   Sepal.Length Sepal.Width Petal.Length Petal.Width Species      resid
# 1          5.1         3.5          1.4         0.2  setosa  0.2209054
# 2          4.9         3.0          1.4         0.2  setosa  0.0209054
# 3          4.7         3.2          1.3         0.2  setosa -0.1382024

一行目、予測値4.879095に対して実際の値は5.1なので、5.1 - 4.879095 = 0.220905 となり一致します(桁数の問題で完全一致ではないですが)
予測値、残差ともに加える場合、pipeフレンドリーな関数なので、このように書くことが出来ますね。もちろん、先の例も全部pipeで受け渡すようにも書けます。

modeleq <- lm(iris$Sepal.Length ~ iris$Petal.Length)

iris %>%  add_predictions(modeleq) %>%  add_residuals(modeleq) %>%  head(n=3)
#   Sepal.Length Sepal.Width Petal.Length Petal.Width Species     pred      resid
# 1          5.1         3.5          1.4         0.2  setosa 4.879095  0.2209054
# 2          4.9         3.0          1.4         0.2  setosa 4.879095  0.0209054
# 3          4.7         3.2          1.3         0.2  setosa 4.838202 -0.1382024

とりあえず今回はここまで。
含まれている関数が結構多いので全部やると時間がかかりそうですが…。

tidyverseに含まれているだけあって、使っていると便利さがわかってきました。
分析を本業にされている方はもう使っているのかもしれませんが…

まぁ、R4DSにも書かれていますけどね。modelr, broomは23章あたりから。

r4ds.had.co.nz

読んでからそろそろ半年なのでいろいろ忘れてきましたが…

本業はRからも分析・機械学習その他からも遠いところにあるので、このあたりは勉強しながら・・・です。コツコツいきましょう。。。

modelrパッケージもぜひぜひ使ってみてください。

Enjoy!