前回の続きです。
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
のようにして書きますね。x
とy
は未定義の変数で良いです。
str()
で見ると、こんな感じですね。
str(y~x) # Class 'formula' language y ~ x # ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
また、左辺と右辺を別々に、~y
, ~x
のように書いて使用することもあります(あとで使います)
~ x + z
や ~ x:z
, ~x * z
等々、色々あるのですが…ここでは割愛。
formula
についてはこちらの記事が良いと思います。私も、大変参考になりました。
また、そのリンク先のこちらの記事も参考になると思います。
modelr::formulas
これもmodelr
パッケージに含まれている関数です。modelr::formulae
という関数もありますが、中身は同じです。
# 確認。 identical(modelr::formulae, modelr::formulas) # [1] TRUE
さて、この関数formulas
はformula
の右辺を複数受け取り、複数の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
これは上に結果があるので良いですかね?右辺に変数を追加してくれます。
予測変数を足していって結果を見る際に便利ですね。例でも、交互作用のあるモデルinteraction
にam
とvs
を追加していますね。
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()
が出しているようなので、たぶん無視してよいと思います(いずれ動かなくなるのだろうか…?)
長くなってしまったのでこの辺で(まだまだかかりますね…)