easystatsについて①: パッケージ群の紹介
これはR Advent Calendar2019の第1日目の記事です。
はじめに
R言語の特徴として
- 統計解析向けの手法がたくさん実装されている
- CRANやGitHubに誰でもパッケージを公開できる
というものがあるかと思います。他にも tidyverse パッケージ群の登場によってデータハンドリング、可視化周りが強くなったり、shiny パッケージでwebアプリが作れたりと最近は色々できるようになっていますが、上記の特徴は大きな特徴の一つと僕は思います。
さて、いろんな手法を別の人が実装した結果、各パッケージのアウトプットが異なるという問題がおこります。（手法が違えば出力が違うのはしょうがないのですが。）例えば、同じデータでいくつかの手法で試してみてその結果を比較したいとき。
model1 <- lm(Sepal.Length ~ Petal.Length, data=iris) model2 <- rstanarm::stan_glm(Sepal.Length ~ Petal.Length, data=iris)
この例でこの比較に意味があるのかはさておき、この２つのモデルに対して回帰係数をそれぞれ見たいとき、
model1$coefficients #> (Intercept) Petal.Length #> 4.3066034 0.4089223 model2$coefficients #> (Intercept) Petal.Length #> 4.3076763 0.4085251
のように見ることもできますが、rstanarm::stan_glm() の出力のリストに coefficients という項目があることは事前知識があるか、model2の中を実際に見てみないとわかりません。
あるいは lme4 パッケージや brms パッケージでランダム効果を入れたモデルを作成したときに、そのランダム回帰係数を見たいと言われたら、 model_brms <- brmss(...) で作ったmodel_brms中を探しに行かなければなりません。
このような統計モデル系のパッケージの間を埋めてくれるのが easystats パッケージ群です。
そのあたりを解決してくれるのが easystats パッケージ群です。tidyverse や tidymodels のようなパッケージ「群」です。インストールをして library(easystats) を実行すれば複数のパッケージがまとめて読み込まれます。つまり tidyverse で dplyr パッケージや tidyr パッケージを組み合わせていたように、複数のパッケージをを組み合わせることが想定されているということです。細かいパッケージの使い方は別記事に譲るとして、本記事では library(easystats) を実行したときに読み込まれるパッケージの紹介と、「例えばこんなことができるよ」的なものを紹介できればと思います。
easystatsパッケージ
↑ロゴがファンキー
インストールと読み込み
個別のパッケージは CRAN にあるものもありますが、まとめてインストールするeasystatsのほうはまだCRANには上がっていないようなのでGitHubからインストールします。余談ですが、GitHubに公開されているRパッケージをインストールするときに devtools::install_github() 関数を使うように記述してあるものが多いですが、現在は remotes::install_github() 関数のほうがおすすめです（devtoolsは依存が多いしそもそもが開発者向けのパッケージなので）。
remotes::install_github("easystats/easystats") library(easystats) #> # Attaching packages #> ✔ insight 0.7.1 ✔ bayestestR 0.4.0.1 #> ✔ performance 0.4.0.1 ✔ parameters 0.3.0.1 #> ✔ see 0.3.0.1 ✔ effectsize 0.0.1
各パッケージの簡単な紹介
insight
insight パッケージは作成したモデルのオブジェクトからモデル式、パラメーターなど各種値をとってきてくれる関数を揃えたパッケージです。例えば、
insight::find_parameters(model2) #> $conditional #> [1] "(Intercept)" "Petal.Length" insight::find_predictors(model2) #> $conditional #> [1] "Petal.Length"
のような感じです。
対応しているクラスは以下の通り。
supported_models() #> [1] "aareg" "aov" "aovlist" #> [4] "bamlss" "bamlss.frame" "bayesx" #> [7] "BBmm" "BBreg" "betareg" #> [10] "BFBayesFactor" "bigglm" "biglm" #> [13] "blavaan" "bracl" "brmsfit" #> [16] "brmultinom" "censReg" "clm" #> [19] "clm2" "clmm" "clmm2" #> [22] "complmrob" "coxme" "coxph" #> [25] "crch" "crq" "feis" #> [28] "felm" "fixest" "flexsurvreg" #> [31] "gam" "Gam" "gamlss" #> [34] "gamm" "gamm4" "gbm" #> [37] "gee" "geeglm" "glimML" #> [40] "glm" "glmmPQL" "glmmTMB" #> [43] "glmrob" "glmRob" "gls" #> [46] "gmnl" "htest" "hurdle" #> [49] "iv_robust" "ivreg" "lavaan" #> [52] "lm" "lm_robust" "lme" #> [55] "lmrob" "lmRob" "logistf" #> [58] "LORgee" "lrm" "mclogit" #> [61] "MCMCglmm" "merMod" "mixed" #> [64] "MixMod" "mlm" "mlogit" #> [67] "mmlogit" "multinom" "ols" #> [70] "plm" "polr" "psm" #> [73] "rlm" "rlmerMod" "rma" #> [76] "rma.uni" "rq" "rqss" #> [79] "speedglm" "speedlm" "stanmvreg" #> [82] "stanreg" "survfit" "survreg" #> [85] "svyglm" "svyolr" "tobit" #> [88] "truncreg" "vgam" "vglm" #> [91] "wbgee" "wblm" "wbm" #> [94] "zeroinfl" "zerotrunc"
bayestestR
bayestestR パッケージはrstanarmやbrmsといったベイズ推定を簡単にやれるパッケージの出力から、事後分布の可視化や信用区間の計算などをサクッとやってくれるパッケージです。
例えば、 insight::get_parameters() で事後分布のサンプルを取得したあとに事後平均、事後中央値、MAP推定量（ bayestestR::map_estimate() ）を可視化するみたいなことが簡単にできるようになります。
ggplot2::ggplot(posteriors, ggplot2::aes(x = Petal.Length)) + ggplot2::geom_density(fill = "orange") + # 青: 事後平均 ggplot2::geom_vline(xintercept = mean(posteriors$Petal.Length), color = "blue", size = 0.5) + # 赤: 事後中央値 ggplot2::geom_vline(xintercept = median(posteriors$Petal.Length), color = "red", size = 0.5) + # 紫: MAP推定値 ggplot2::geom_vline(xintercept = bayestestR::map_estimate(posteriors$Petal.Length), color ="purple", size = 0.5)
performance
performance パッケージは文字通り、モデルのパフォーマンス (r2, MSE, Loglossなど) を簡単に算出できるパッケージです。
performance::r2(model1) #> # R2 for Linear Regression #> #> R2: 0.760 #> adj. R2: 0.758 performance::performance_rmse(model2) #> [1] 0.4043614
あるいは多重共線性を調べることなんかもできます。
data(mtcars) model3 <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars) colli <- performance::check_collinearity(model3) colli #> # Check for Multicollinearity #> #> Low Correlation #> #> Parameter VIF Increased SE #> gear 1.53 1.24 #> #> Moderate Correlation #> #> Parameter VIF Increased SE #> wt 5.05 2.25 #> cyl 5.41 2.33 #> disp 9.97 3.16 plot(colli)
parameters
parameters パッケージの説明は上記の図が一番わかりやすいのではないでしょうか。
いろんなパッケージの出力に対して、統一的な表を作成してくれます。
model4 <- lm(Sepal.Width ~ Petal.Length * Species + Petal.Width, data = iris) parameters::model_parameters(model4) #> Parameter | Coefficient | SE | 95% CI | t | df | p #> ------------------------------------------------------------------------------------------------ #> (Intercept) | 2.89 | 0.36 | [ 2.18, 3.60] | 8.01 | 143 | < .001 #> Petal.Length | 0.26 | 0.25 | [-0.22, 0.75] | 1.07 | 143 | 0.287 #> Species [versicolor] | -1.66 | 0.53 | [-2.71, -0.62] | -3.14 | 143 | 0.002 #> Species [virginica] | -1.92 | 0.59 | [-3.08, -0.76] | -3.28 | 143 | 0.001 #> Petal.Width | 0.62 | 0.14 | [ 0.34, 0.89] | 4.41 | 143 | < .001 #> Petal.Length * Species [versicolor] | -0.09 | 0.26 | [-0.61, 0.42] | -0.36 | 143 | 0.721 #> Petal.Length * Species [virginica] | -0.13 | 0.26 | [-0.64, 0.38] | -0.50 | 143 | 0.618
parameters::parameters_selection() 変数選択なんかもできます。
library(dplyr) model5_params <- rstanarm::stan_glm(mpg ~ ., data = mtcars, refresh = 0) %>% parameters::parameters_selection() %>% parameters::model_parameters() model5_params #> Registered S3 method overwritten by 'xts': #> method from #> as.zoo.xts zoo #> Parameter | Median | 89% CI | pd | % in ROPE | Rhat | ESS | Prior #> ----------------------------------------------------------------------------------------------- #> (Intercept) | 20.16 | [-2.87, 42.17] | 92.60% | 1.20% | 1.001 | 2320 | Normal (0 +- 60.27) #> wt | -3.93 | [-5.93, -1.96] | 99.78% | 0.57% | 1.001 | 2571 | Normal (0 +- 15.40) #> cyl | -0.49 | [-1.77, 0.85] | 73.17% | 46.85% | 1.001 | 2610 | Normal (0 +- 8.44) #> hp | -0.02 | [-0.04, 0.00] | 89.80% | 100% | 1.000 | 3283 | Normal (0 +- 0.22) #> am | 2.92 | [ 0.09, 5.88] | 94.65% | 7.50% | 1.000 | 2836 | Normal (0 +- 15.07) #> qsec | 0.79 | [-0.14, 1.80] | 90.80% | 35.62% | 1.001 | 2363 | Normal (0 +- 8.43) #> disp | 0.01 | [-0.01, 0.03] | 87.45% | 100% | 1.002 | 3019 | Normal (0 +- 0.12)
see
see パッケージは、easystats の他のパッケージと組み合わせて、いい感じの図を提供してくれるパッケージです。ggplot2 と組み合わせて使えるカラーテーマや、 plot() で簡単に図を出せるようにしてくれます。例えば、先程の parametes パッケージと組み合わせて、
plot(model5_params)
また、bayestestR パッケージなどと組み合わせる例もGitHubで紹介されています。
effectsize
effectsize パッケージは、文字通り効果量を計算してくれるパッケージです。
## [1] -4.21 hedges_g(iris$Sepal.Length, iris$Sepal.Width) ## [1] -4.2 glass_delta(iris$Sepal.Length, iris$Sepal.Width) ## [1] -3.36
標準化回帰係数を出すこともできます。
std_model <- effectsize::standardize(lm(Sepal.Length ~ Species, data = iris)) coef(std_model) #> (Intercept) Speciesversicolor Speciesvirginica #> -1.011191 1.123099 1.910475
まとめ
本記事では、 library(easystats) を実行したときに読み込まれる以下のパッケージについて簡単に紹介しました。
- insight
- bayestestR
- parameters
- performance
- see
- effectsize
各パッケージでできることはまだまだたくさんありますし、この他にも easystats には report , correlation , estimate パッケージがあります。
まだドキュメントが十分には感じないのですが、統計モデルのパッケージの種類に関わらず、よく使われる指標を簡単に算出できるのは便利に感じました。
Enjoy!
