久々の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について書きました。
このときはstanの結果についてみてみましたが、tidyはlmやglmの結果をtidyにしてくれます（実際はtidyでtidyMCMCと同じようにstanの結果も見れる）
modelrはlmやglmなどの解析をtidyにしてくれます。また、ブートストラップやらCVなどに使える関数も色々入っており、broomとあわせることで解析が楽になるのではないかと思います。
実際に使ってみる
とりあえず、実際に使って結果を見る方が早いので。 modelrの例文だとmtcarsを使っていることが多いのですが、簡単のためにirisのデータを使ってみます。
plot(iris$Sepal.Length ~ iris$Petal.Length)
こんな図になりますね。特に面白みのない図です。
さて、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章あたりから。
読んでからそろそろ半年なのでいろいろ忘れてきましたが…
本業はRからも分析・機械学習その他からも遠いところにあるので、このあたりは勉強しながら・・・です。コツコツいきましょう。。。
modelrパッケージもぜひぜひ使ってみてください。
Enjoy!
