前回の続きです。
niszet.hatenablog.com
fit_with
ヘルプのパクリになってしまいますが、この関数は予測変数が複数ないとあまり意味がないので、使用例と同じmtcarsデータセットを素直に使用します。
用語として、説明変数なのか予測変数なのかちょっとよくわかっていませんが、moderlrでは式の右辺をpredictorと表現しているため予測変数と書いています。統計ガチ勢の方からの優しいマサカリをお待ちしております。。。
mtcarsデータセットの詳細については?mtcarsでヘルプを参照するか、こちらを参照ください。
R: Motor Trend Car Road Tests
でもリンク先見るのが面倒くさいので、一部抜粋して貼っておきます。
mtcarsデータセットは下記のようなデータになっています。
head(mtcars, n=3)
ヘルプから、それぞれの列は下記を示しているそうです。
[, 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)
))
色々と知らない関数が使われていて、わかりづらいので１つずつ見ていきます。
formulaについて
これは深入りすると結構重い話になるのでさらっと。
前回の例のように、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についてはこちらの記事が良いと思います。私も、大変参考になりました。
ill-identified.hatenablog.com
また、そのリンク先のこちらの記事も参考になると思います。
m884.hateblo.jp
modelr::formulas
これもmodelrパッケージに含まれている関数です。modelr::formulaeという関数もありますが、中身は同じです。
identical(modelr::formulae, modelr::formulas)
さて、この関数formulasはformulaの右辺を複数受け取り、複数のformulaを作ることが出来ます。…なんのこっちゃ？
実際に先の例からこの部分を取り出して実行してみます。
tmp <- formulas(~disp,
model1 = ~drat + cyl,
interaction = ~drat * cyl,
full = add_predictors(interaction, ~am, ~vs))
tmp
…と、このようにリストが返却されます。formulasの引数として、要素の名前（モデルの名前）を指定して式の右辺を与えることが出来ます。
左辺の変数は固定のようです（ヘルプより）ので、応答変数が複数ある場合は別々に行う必要がありますね。
さて、この関数の引数にもまた関数が…
add_predictors
これもmodelr(ry
これは上に結果があるので良いですかね？右辺に変数を追加してくれます。
予測変数を足していって結果を見る際に便利ですね。例でも、交互作用のあるモデルinteractionにamとvsを追加していますね。
interaction <- ~drat * cyl
interaction
full <- add_predictors(interaction, ~am, ~vs)
full
モデルから変数を取り除く方はできないと思います(あれば教えてください。)
さて、最初の実行結果はリストだったので、要素にアクセスすれば、
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)
出力が1つのデータフレームになりますが、このままではどのモデルだったのかがわからなくなってしまいますね。もう少し工夫が必要そうです。
なお、broom::tidy以外にも、broom::glance, broom::augmentに食わせてあげることもできます。
disp_fits %>% map(glance)
augmentは結果が長くなるので1つだけ（かつhead()に食わせて）、下記のようになりますね。
なお、このbroom::augmentは使用すると
Warning messages:
1: Deprecated: please use `purrr::possibly()` instead
のようにwarningが出てきます。しかしこれはbroom::augmentが内部で使用しているbroom:::augment_columns()が出しているようなので、たぶん無視してよいと思います（いずれ動かなくなるのだろうか…？）
長くなってしまったのでこの辺で（まだまだかかりますね…）
