niszetの日記

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

(R) 隠れtidyverse modelr編 とりあえず使ってみる 2回目 fit_withと関連する関数

前回の続きです。

niszet.hatenablog.com

fit_with

ヘルプのパクリになってしまいますが、この関数は予測変数が複数ないとあまり意味がないので、使用例と同じmtcarsデータセットを素直に使用します。
用語として、説明変数なのか予測変数なのかちょっとよくわかっていませんが、moderlrでは式の右辺をpredictorと表現しているため予測変数と書いています。統計ガチ勢の方からの優しいマサカリをお待ちしております。。。

mtcarsデータセットの詳細については?mtcarsでヘルプを参照するか、こちらを参照ください。 R: Motor Trend Car Road Tests

でもリンク先見るのが面倒くさいので、一部抜粋して貼っておきます。
mtcarsデータセットは下記のようなデータになっています。

head(mtcars, n=3)
#                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
# Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
# Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
# Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1

ヘルプから、それぞれの列は下記を示しているそうです。

[, 1]   mpg Miles/(US) gallon
[, 2]  cyl Number of cylinders
[, 3]  disp    Displacement (cu.in.)
[, 4]  hp  Gross horsepower
[, 5]  drat    Rear axle ratio
[, 6]  wt  Weight (1000 lbs)
[, 7]  qsec    1/4 mile time
[, 8]  vs  V/S
[, 9]  am  Transmission (0 = automatic, 1 = manual)
[,10]  gear    Number of forward gears
[,11]  carb    Number of carburetors

さて、前回同様に、下記のライブラリをロードしておきます。

library(tidyverse)
library(broom)
library(modelr)

さて、ヘルプにあった例をそのまま実行してみます。…と思ったのですが、最初のモデルだけmodel1という名前に変えています。結果の差異を見たかったのでその名残…。

disp_fits <- mtcars %>% fit_with(lm, formulas(~disp,
model1 = ~drat + cyl,
interaction = ~drat * cyl,
full = add_predictors(interaction, ~am, ~vs)
))

色々と知らない関数が使われていて、わかりづらいので1つずつ見ていきます。

formulaについて

これは深入りすると結構重い話になるのでさらっと。
前回の例のように、yが応答変数xが予測変数であるような場合は

y~x
# y ~ x

のようにして書きますね。xyは未定義の変数で良いです。

str()で見ると、こんな感じですね。

str(y~x)
# Class 'formula'  language y ~ x
#   ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 

また、左辺と右辺を別々に、~y, ~x のように書いて使用することもあります(あとで使います) ~ x + z~ x:z, ~x * z 等々、色々あるのですが…ここでは割愛。

formulaについてはこちらの記事が良いと思います。私も、大変参考になりました。

ill-identified.hatenablog.com

また、そのリンク先のこちらの記事も参考になると思います。

m884.hateblo.jp

modelr::formulas

これもmodelrパッケージに含まれている関数です。modelr::formulaeという関数もありますが、中身は同じです。

# 確認。
identical(modelr::formulae, modelr::formulas)
# [1] TRUE

さて、この関数formulasformulaの右辺を複数受け取り、複数のformulaを作ることが出来ます。…なんのこっちゃ?
実際に先の例からこの部分を取り出して実行してみます。

tmp <- formulas(~disp,
         model1 = ~drat + cyl,
         interaction = ~drat * cyl,
         full = add_predictors(interaction, ~am, ~vs))

tmp
# $model1
# disp ~ drat + cyl
# 
# $interaction
# disp ~ drat * cyl
# 
# $full
# disp ~ drat * cyl + am + vs

…と、このようにリストが返却されます。formulasの引数として、要素の名前(モデルの名前)を指定して式の右辺を与えることが出来ます。
左辺の変数は固定のようです(ヘルプより)ので、応答変数が複数ある場合は別々に行う必要がありますね。

さて、この関数の引数にもまた関数が…

add_predictors

これもmodelr(ry

これは上に結果があるので良いですかね?右辺に変数を追加してくれます。 予測変数を足していって結果を見る際に便利ですね。例でも、交互作用のあるモデルinteractionamvsを追加していますね。

interaction <- ~drat * cyl
interaction
# ~drat * cyl
full <-  add_predictors(interaction, ~am, ~vs)
full
# ~drat * cyl + am + vs

モデルから変数を取り除く方はできないと思います(あれば教えてください。)

さて、最初の実行結果はリストだったので、要素にアクセスすれば、

disp_fits$model1
# 
# Call:
# .f(formula = disp ~ drat + cyl, data = data)
# 
# Coefficients:
# (Intercept)         drat          cyl  
#       18.72       -35.83        55.09  

となりますね。例によってtidyで整形すれば、

> broom::tidy(disp_fits$model1)
         term  estimate  std.error  statistic      p.value
1 (Intercept)  18.71579 127.836537  0.1464041 8.846154e-01
2        drat -35.83058  25.150896 -1.4246243 1.649392e-01
3         cyl  55.09059   7.529809  7.3163322 4.652654e-08

このように。 listの処理なので、purrr::mapの出番ですね。

> disp_fits %>%  map(tidy)
$model1
         term  estimate  std.error  statistic      p.value
1 (Intercept)  18.71579 127.836537  0.1464041 8.846154e-01
2        drat -35.83058  25.150896 -1.4246243 1.649392e-01
3         cyl  55.09059   7.529809  7.3163322 4.652654e-08

$interaction
         term    estimate std.error  statistic    p.value
1 (Intercept) -193.487202 335.44306 -0.5768109 0.56867604
2        drat   20.389801  85.87195  0.2374442 0.81404100
3         cyl   88.674520  49.59003  1.7881522 0.08457908
4    drat:cyl   -9.154652  13.35808 -0.6853269 0.49877264

$full
         term    estimate std.error  statistic   p.value
1 (Intercept) -164.832719 385.89073 -0.4271487 0.6727877
2        drat   27.577726  90.83136  0.3036146 0.7638379
3         cyl   75.843245  56.04986  1.3531390 0.1876606
4          am  -38.217318  31.42983 -1.2159567 0.2349265
5          vs  -15.529580  40.06479 -0.3876117 0.7014585
6    drat:cyl   -6.969767  14.00552 -0.4976442 0.6229177

こんな感じに出力されました。便利ですね。mapではなくmap_dfrに与えることで下記のようになります。

disp_fits %>%  map_dfr(tidy)
#           term    estimate  std.error  statistic      p.value
# 1  (Intercept)   18.715788 127.836537  0.1464041 8.846154e-01
# 2         drat  -35.830577  25.150896 -1.4246243 1.649392e-01
# 3          cyl   55.090585   7.529809  7.3163322 4.652654e-08
# 4  (Intercept) -193.487202 335.443063 -0.5768109 5.686760e-01
# 5         drat   20.389801  85.871954  0.2374442 8.140410e-01
# 6          cyl   88.674520  49.590030  1.7881522 8.457908e-02
# 7     drat:cyl   -9.154652  13.358081 -0.6853269 4.987726e-01
# 8  (Intercept) -164.832719 385.890732 -0.4271487 6.727877e-01
# 9         drat   27.577726  90.831360  0.3036146 7.638379e-01
# 10         cyl   75.843245  56.049855  1.3531390 1.876606e-01
# 11          am  -38.217318  31.429835 -1.2159567 2.349265e-01
# 12          vs  -15.529580  40.064790 -0.3876117 7.014585e-01
# 13    drat:cyl   -6.969767  14.005520 -0.4976442 6.229177e-01

出力が1つのデータフレームになりますが、このままではどのモデルだったのかがわからなくなってしまいますね。もう少し工夫が必要そうです。

なお、broom::tidy以外にも、broom::glance, broom::augmentに食わせてあげることもできます。

disp_fits %>%  map(glance) 
# $model1
#   r.squared adj.r.squared   sigma statistic      p.value df    logLik      AIC      BIC deviance df.residual
# 1 0.8258511     0.8138408 53.4748  68.76206 9.848099e-12  3 -171.1657 350.3315 356.1944 82927.08          29
# 
# $interaction
#   r.squared adj.r.squared    sigma statistic      p.value df    logLik      AIC      BIC deviance df.residual
# 1  0.828724      0.810373 53.97056  45.15962 7.426765e-11  4 -170.8996 351.7991 359.1278 81559.01          28
# 
# $full
#   r.squared adj.r.squared    sigma statistic      p.value df    logLik      AIC      BIC deviance df.residual
# 1 0.8381416     0.8070149 54.44635  26.92684 1.658079e-09  6 -169.9947 353.9894 364.2496 77074.53          26

augmentは結果が長くなるので1つだけ(かつhead()に食わせて)、下記のようになりますね。

# disp_fits$model1 %>%  augment() %>%  head()
#           .rownames disp drat cyl  .fitted  .se.fit     .resid       .hat   .sigma      .cooksd .std.resid
# 1         Mazda RX4  160 3.90   6 209.5200 11.59803 -49.520050 0.04704028 53.57031 0.0148068727 -0.9486251
# 2     Mazda RX4 Wag  160 3.90   6 209.5200 11.59803 -49.520050 0.04704028 53.57031 0.0148068727 -0.9486251
# 3        Datsun 710  108 3.85   4 101.1304 15.94755   6.869592 0.08893844 54.40433 0.0005894354  0.1345885
# 4    Hornet 4 Drive  258 3.08   6 238.9011 16.90633  19.098877 0.09995406 54.28819 0.0052464695  0.3764665
# 5 Hornet Sportabout  360 3.15   8 346.5742 13.68171  13.425848 0.06546100 54.35801 0.0015748949  0.2597131
# 6           Valiant  225 2.76   6 250.3669 23.99234 -25.366907 0.20130145 54.15633 0.0236699055 -0.5307951

なお、このbroom::augmentは使用すると

Warning messages:
1: Deprecated: please use `purrr::possibly()` instead 

のようにwarningが出てきます。しかしこれはbroom::augmentが内部で使用しているbroom:::augment_columns()が出しているようなので、たぶん無視してよいと思います(いずれ動かなくなるのだろうか…?)

長くなってしまったのでこの辺で(まだまだかかりますね…)