★データ解析備忘録★

ゆる〜い技術メモ

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)

この例でこの比較に意味があるのかはさておき、この2つのモデルに対して回帰係数をそれぞれ見たいとき、

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 パッケージ群です。tidyversetidymodels のようなパッケージ「群」です。インストールをして library(easystats) を実行すれば複数のパッケージがまとめて読み込まれます。つまり tidyversedplyr パッケージや tidyr パッケージを組み合わせていたように、複数のパッケージをを組み合わせることが想定されているということです。細かいパッケージの使い方は別記事に譲るとして、本記事では library(easystats) を実行したときに読み込まれるパッケージの紹介と、「例えばこんなことができるよ」的なものを紹介できればと思います。

easystatsパッケージ

f:id:songcunyouzai:20191201125950p:plain f:id:songcunyouzai:20191201130300p:plain ↑ロゴがファンキー

インストールと読み込み

個別のパッケージは 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

easystats.github.io

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

easystats.github.io

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)

f:id:songcunyouzai:20191201131915p:plain

performance

easystats.github.io

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)

f:id:songcunyouzai:20191201133012p:plain

parameters

easystats.github.io f:id:songcunyouzai:20191201133319p:plain

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

easystats.github.io

see パッケージは、easystats の他のパッケージと組み合わせて、いい感じの図を提供してくれるパッケージです。ggplot2 と組み合わせて使えるカラーテーマや、 plot() で簡単に図を出せるようにしてくれます。例えば、先程の parametes パッケージと組み合わせて、

plot(model5_params)

f:id:songcunyouzai:20191201134959p:plain

また、bayestestR パッケージなどと組み合わせる例もGitHubで紹介されています。

github.com

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!

「白金鉱業 Meetup Vol.12(白金鉱業 x IBM社内勉強会dsn 合同企画回)」に参加しました。

白金鉱業 Meetup Vol.12(白金鉱業 x IBM社内勉強会dsn 合同企画回)に参加しました。 白金鉱業 Meetupの参加自体始初めてだったのですが、非常に楽しめました。

brainpad-meetup.connpass.com

続きを読む

「データアーキテクト(データ整備人)を”前向きに”考える会」に参加しました

はじめに

11/27に開催された「データアーキテクト(データ整備人)を”前向きに”考える会」にブログ枠として参加させて頂いたのでイベントの様子をレポートとして書きます。

analytics-and-intelligence.connpass.com

会場: 株式会社オプト

続きを読む

RからGoの関数をつかう→はやい

これはなに

こんな話があります(だいぶ昔

qiita.com

qiita.com

で、RからRustを使おうというのが ゆたにサンの記事。*1

notchained.hatenablog.com

それで、この記事はRからGoを使おうって話です。

速度の比較をやりたいので、本記事ではフィボナッチ数列の関数を書いて比較します。

※この記事は別に「みんなGoを使おうぜ」みたいな趣旨ではないです。現状だとC++使うほうが応用範囲は広いと思うので

*1:余談ですが、僕はこの記事に関する発表をTokyo.Rで見た覚えがあるんですが、もう3年以上前のことなのですね。。。光陰矢の如し。

続きを読む

Rおじさん、Lisperへの道 (Common Lisp環境設定編)

はじめに

Twitter上でこんなやりとりがあって、Lispをちゃんと勉強しようかなと思った次第。これまでEmacsを使って設定ファイルも(汚いけど)それなりに書いてきたので基礎を知ってるくらいなので、ここは一つCommon Lispをちゃんとやっておきたいなというのが主な動機です。

ところで、このブログは「データ解析周りを中心に、プログラミング全般について書く」と明言してるわけですがこれからやるCommon Lispがどのくらい役に立つかは未知数。

続きを読む