glmやglmerによるパラメータ推定値と信頼区間をggplot2で描画する (easystats)

# 解析結果の視覚的な確認
![glmer2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/535720/5a08e0a8-8110-fe81-0d6f-f3af46a8adfe.png)



一般化線形モデルや一般化線形混合モデルのパラメータ推定をRで行う場合、よく用いられるのはglmやglmer（lmer）だと思います。

これらの関数を実行して得られるもっとも主要な結果はモデルにおけるパラメータの最尤推定値です。パラメータの推定値は`print`や`summary`で確認することが出来ますが、数字だけではイメージを掴みづらい場合があります。そのようなときには視覚的に推定値を確認するとわかりやすいです。
しかし、推定値をプロットしようとするとオブジェクトからのデータ抽出や95%信頼区間の算出が結構面倒です。またそれらの結果をggplot2で描画しようとするとコードが煩雑になります。

easystatsパッケージ群に含まれるparametersパッケージとggplot2と組み合わせて使うと簡潔なコードで推定値を描画できます。ここではその一例を紹介したいと思います。

例により、詳しい統計学的な説明の方はしない（できない）ので使い方のメモとして見てください。レポートなどに用いる際にはパッケージのヘルプやレファレンス、ソースコード等を見て導出過程や計算方法を確認するのが良いと思います。
**また信頼区間の算出方法については注意点があるので後述の説明も参考にしてください。**


# インストール
CRANとGitHubどちらからでもインストールできます。GitHubは更新が速いので安定版を使いたいのであればCRANのほうがおすすめです。

- CRAN

以下のコマンドでインストールできます。

```R
install.packages("parameters")
```

- GitHub

easystatsごとインストールすることでインストールできます。

```R
remotes::install_github("easystats/easystats")
```

すでにeasystatsに含まれるパッケージを部分的にインストールしている場合には、依存関係関連でエラーや警告が出たりします。そのため上記コマンドを実行する前にいったんそれらをアンインストールしたほうが良いです。


# `parameters::model_parameters`について
今回はパラメータの推定値および信頼区間の取得に`parameters::model_parameters`を使用します。

この関数を`glm`等の解析結果のオブジェクトに適用すると以下のような出力が得られます：

```R
> model <- glm(vs ~ wt + mpg, family = "binomial", data = mtcars)
> (model_parameters <- model_parameters(model))
Parameter   | Coefficient |   SE |         95% CI |     z | df |     p
----------------------------------------------------------------------
(Intercept) |      -12.54 | 8.47 | [-31.91, 2.92] | -1.48 | 29 | 0.139
wt          |        0.58 | 1.18 | [ -1.92, 2.94] |  0.49 | 29 | 0.623
mpg         |        0.52 | 0.26 | [  0.09, 1.17] |  2.01 | 29 | 0.044
```

コンソール出力では上記のようになっていますが、実際の中身は以下のようなデータフレームです：

```R
> as.data.frame(model_parameters)
    Parameter Coefficient        SE       CI_low  CI_high         z df_error         p
1 (Intercept) -12.5412218 8.4660329 -31.90842474 2.921046 -1.481358       29 0.1385113
2          wt   0.5828598 1.1844650  -1.91815266 2.936813  0.492087       29 0.6226579
3         mpg   0.5240640 0.2604188   0.09326514 1.165436  2.012389       29 0.0441789
```

データがこのような構造になっているためggplot2とうまく組み合わせて使うことが出来ます。

### 注意：信頼区間の値について
`parameters::model_parameters`を使って算出される信頼区間の値には注意が必要です。

デフォルトの設定では、対象が`glm`の結果（クラス`glm`のオブジェクト）の場合には**`ststs::confint`**で、対象が`glmer`の結果（クラス`merMod`のオブジェクト）の場合には**`parameters::ci_wald(dof = Inf)`**で信頼区間が算出されます。
`ststs::confint`についてはSEを使わない方法で信頼区間を算出しているようです（[参考](https://github.com/cran/MASS/blob/master/R/confint.R)）。よって、「Estimate ± 1.959964... x SE」というよく見る計算式による計算結果とは必ずしも一致しません。対して`parameters::ci_wald(dof = Inf)`による計算結果は「Estimate ± 1.959964... x SE」による計算結果と一致します（[参考](https://github.com/easystats/parameters/blob/master/R/ci_wald.R)。厳密にはt分布が使用されていますが自由度無限大なので正規分布のquantileと一致しています）。

ちなみにz値は`glm`や`glmer`で計算された値がそのまま使われるようです。p値の方は`glm`に対しては`glm`で計算された値が、`glmer`に対しては`parameters::p_value_wald`で計算された値（実際は`glmer`で計算された値と一致する？）が使用されるようです。
（これらはヘルプには書いていなかったので[GitHub上のソースコード](https://github.com/easystats/parameters/tree/master/R)で確認しました。間違ってはいないと思いますがご自分でチェックしたい方はソースコードを見てみてください）

算出方法が納得いかなかったり、満足いかない方は別の方法で信頼区間を出す方が良いと思います。また一方で、最尤推定値のSEや信頼区間の計算方法および意味を詳しく理解していないのであれば専門家が作成した関数を使ってしまうのも手かもしれません。
ここらへんは自分の用途や統計手法の使い方に対する態度によると思います。

ここでは説明していませんが信頼区間の算出方法にブーツストラップ法も指定できるようです。そちらを勉強してそっちを使ってしまうのもありかもしれません。


# 推定値の視覚化
ここでは以下のパッケージを使用します。

```R
library(parameters)
library(tidyverse)
library(lme4)
library(glmmTMB)
```

テストデータ、テストコードについては各パッケージのヘルプを参考にしました。一部収束エラーが出ていますが例示としては問題ないのでそのままにしてあります。

基本的には、パラメータ推定→`parameters::model_parameters`の実行→ggplot2で描画、という流れで推定値と信頼区間をプロットします。
推定値と信頼区間は`ggplot2::geom_pointrange`というとても便利な関数を使って描画します（[こちらのページ](https://peripatopsis.jimdofree.com/2020/02/26/geom-pointrange%E3%81%A7%E6%8E%A8%E5%AE%9A%E5%80%A4%E3%81%A895-%E4%BF%A1%E9%A0%BC%E5%8C%BA%E9%96%93%E3%81%AE%E5%9B%B3%E7%A4%BA/)で紹介されています）。またエラーバーを描画するgeom_*関数は[過去の記事](https://qiita.com/ocean_f/items/4d28b8707ff5c7b977e0)で列挙してあります。


## glm
`glm`の結果を描画すると以下のようになります。

```R
> glm_a <- glm(vs ~ wt + mpg, family = "binomial", data = mtcars)
> 
> (glm_a_param <- model_parameters(glm_a))
Parameter   | Coefficient |   SE |         95% CI |     z | df |     p
----------------------------------------------------------------------
(Intercept) |      -12.54 | 8.47 | [-31.91, 2.92] | -1.48 | 29 | 0.139
wt          |        0.58 | 1.18 | [ -1.92, 2.94] |  0.49 | 29 | 0.623
mpg         |        0.52 | 0.26 | [  0.09, 1.17] |  2.01 | 29 | 0.044
> 
> ggplot(data = glm_a_param, aes(x = Parameter)) + 
+   geom_pointrange(aes(y = Coefficient, ymin = CI_low, ymax = CI_high)) + 
+   geom_hline(yintercept = 0, lty = 2) + ylab("Estimate")
```

![glm_a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/535720/6ab9ae8c-485f-dd7b-37d6-0bf30ebfd074.png)


説明変数を正規化すると以下のようになります。軸を入れ替えて見た目を少し変えています。

```R
> glm_b <- glm(vs ~ scale(wt) + scale(mpg), family = "binomial", data = mtcars)
> 
> (glm_b_param <- model_parameters(glm_b))
Parameter   | Coefficient |   SE |        95% CI |     z | df |     p
---------------------------------------------------------------------
(Intercept) |       -0.14 | 0.51 | [-1.15, 0.90] | -0.27 | 29 | 0.787
wt          |        0.57 | 1.16 | [-1.88, 2.87] |  0.49 | 29 | 0.623
mpg         |        3.16 | 1.57 | [ 0.56, 7.02] |  2.01 | 29 | 0.044
> 
> ggplot(data = glm_b_param, aes(x = reorder(Parameter, desc(Parameter)))) + 
+   geom_pointrange(aes(y = Coefficient, ymin = CI_low, ymax = CI_high)) + 
+   geom_hline(yintercept = 0, lty = 2) + coord_flip() + xlab("Parameter") + ylab("Estimate")
```

![glm_b.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/535720/b65ebcdc-9402-f296-a2de-28a3dad9727e.png)

描画の参考：
[ggplot2で縦軸と横軸をひっくり返したい](https://qiita.com/kazutan/items/5c4ee243c48a64b44b2d)
[Reverse the order of a categorical axis in ggplot2](https://gist.github.com/jennybc/6f3fa527b915b920fdd5)


また以下のようにすることで`parameters::ci_wald`による信頼区間を描画することも出来ます。`dof`を指定しなかった場合には（自由度） = （残差自由度）であるt分布のquantileが使用されます（[参考](https://github.com/easystats/parameters/blob/master/R/ci_wald.R)）。
赤が`ststs::confint`による信頼区間、青が`parameters::ci_wald`による信頼区間です。

```R
> glm_c <- glm(vs ~ scale(wt) + scale(mpg), family = "binomial", data = mtcars)
> 
> (glm_c_param <- model_parameters(glm_c))
Parameter   | Coefficient |   SE |        95% CI |     z | df |     p
---------------------------------------------------------------------
(Intercept) |       -0.14 | 0.51 | [-1.15, 0.90] | -0.27 | 29 | 0.787
wt          |        0.57 | 1.16 | [-1.88, 2.87] |  0.49 | 29 | 0.623
mpg         |        3.16 | 1.57 | [ 0.56, 7.02] |  2.01 | 29 | 0.044
> 
> (glm_c_ci <- ci_wald(glm_c, dof = Inf)) # ci_waldで信頼区間を算出
    Parameter CI     CI_low   CI_high
1 (Intercept) 95 -1.1316755 0.8571893
2   scale(wt) 95 -1.7011940 2.8418010
3  scale(mpg) 95  0.0822833 6.2347296
> 
> (glm_c_param <- merge(glm_c_param, glm_c_ci, by = "Parameter")) # model_parametersの結果とマージ
    Parameter Coefficient        SE   CI_low.x CI_high.x          z df_error         p CI   CI_low.y CI_high.y
1 (Intercept)  -0.1372431 0.5073728 -1.1518352 0.9009334 -0.2704976       29 0.7867775 95 -1.1316755 0.8571893
2  scale(mpg)   3.1585064 1.5695304  0.5621041 7.0240230  2.0123894       29 0.0441789 95  0.0822833 6.2347296
3   scale(wt)   0.5703035 1.1589486 -1.8768307 2.8735462  0.4920870       29 0.6226579 95 -1.7011940 2.8418010
> 
> ggplot(data = glm_c_param, aes(x = Parameter)) + 
+   geom_pointrange(aes(y = Coefficient, ymin = CI_low.x, ymax = CI_high.x), col = "red", position = position_nudge(x = -0.05)) + 
+   geom_pointrange(aes(y = Coefficient, ymin = CI_low.y, ymax = CI_high.y), col = "blue", position = position_nudge(x = 0.05)) + 
+   geom_hline(yintercept = 0, lty = 2) + ylab("Estimate")
```

![glm_c.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/535720/57b4ffae-7420-8323-f161-3cd40e039861.png)

描画の参考：
[Nudge points a fixed distance - ggplot2](https://ggplot2.tidyverse.org/reference/position_nudge.html)


## glmer
`glmer`の結果を描画すると以下のようになります。

```R
> form <- TICKS ~ YEAR + HEIGHT + (1 | BROOD) + (1 | LOCATION)
> glmer1  <- glmer(form, family = "poisson", data = grouseticks)
 警告メッセージ: 
1:  checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  で: 
  Model failed to converge with max|grad| = 0.0397798 (tol = 0.002, component 1)
2:  checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  で: 
  Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?
> 
> (glmer1_param <- model_parameters(glmer1))
Parameter   | Coefficient |   SE |         95% CI |      z |  df |      p
-------------------------------------------------------------------------
(Intercept) |       11.35 | 0.48 | [10.41, 12.30] |  23.54 | 397 | < .001
YEAR [96]   |        1.17 | 0.23 | [ 0.72,  1.62] |   5.08 | 397 | < .001
YEAR [97]   |       -0.98 | 0.25 | [-1.48, -0.48] |  -3.84 | 397 | < .001
HEIGHT      |       -0.02 | 0.00 | [-0.03, -0.02] | -24.08 | 397 | < .001
> 
> ggplot(data = glmer1_param, aes(x = Parameter)) + 
+   geom_pointrange(aes(y = Coefficient, ymin = CI_low, ymax = CI_high)) + 
+   geom_hline(yintercept = 0, lty = 2) + ylab("Estimate")
```

![glmer1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/535720/1555a719-c5f7-292a-4049-32b6f7e4c4b2.png)



交互作用がある場合にも簡単に描画できます。このように説明変数が複数の場合には`parameters::model_parameters`はとても便利です。

```R
> glmer2 <- glmer(count ~ spp * mined + (1 | site), family = "poisson", data = Salamanders)
 警告メッセージ: 
 checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  で: 
  Model failed to converge with max|grad| = 0.00689362 (tol = 0.002, component 1)
> 
> (glmer2_param <- model_parameters(glmer2))
Parameter                | Coefficient |   SE |         95% CI |     z |  df |      p
-------------------------------------------------------------------------------------
(Intercept)              |       -3.37 | 0.72 | [-4.79, -1.96] | -4.66 | 629 | < .001
spp [PR]                 |        0.91 | 0.83 | [-0.70,  2.53] |  1.11 | 629 | 0.268 
spp [DM]                 |        2.25 | 0.73 | [ 0.81,  3.69] |  3.07 | 629 | 0.002 
spp [EC-A]               |        0.69 | 0.85 | [-0.98,  2.37] |  0.81 | 629 | 0.417 
spp [EC-L]               |        1.70 | 0.76 | [ 0.22,  3.19] |  2.24 | 629 | 0.025 
spp [DES-L]              |        2.52 | 0.72 | [ 1.10,  3.94] |  3.48 | 629 | < .001
spp [DF]                 |        2.52 | 0.72 | [ 1.10,  3.94] |  3.48 | 629 | < .001
mined [no]               |        4.11 | 0.75 | [ 2.64,  5.58] |  5.49 | 629 | < .001
spp [PR] * mined [no]    |       -2.49 | 0.86 | [-4.17, -0.81] | -2.90 | 629 | 0.004 
spp [DM] * mined [no]    |       -2.15 | 0.74 | [-3.61, -0.69] | -2.89 | 629 | 0.004 
spp [EC-A] * mined [no]  |       -1.53 | 0.87 | [-3.24,  0.18] | -1.75 | 629 | 0.080 
spp [EC-L] * mined [no]  |       -1.12 | 0.77 | [-2.62,  0.39] | -1.46 | 629 | 0.145 
spp [DES-L] * mined [no] |       -1.95 | 0.73 | [-3.39, -0.51] | -2.66 | 629 | 0.008 
spp [DF] * mined [no]    |       -2.67 | 0.74 | [-4.11, -1.22] | -3.61 | 629 | < .001
> 
> ggplot(data = glmer2_param, aes(x = Parameter)) + 
+   geom_pointrange(aes(y = Coefficient, ymin = CI_low, ymax = CI_high)) + 
+   geom_hline(yintercept = 0, lty = 2) + ylab("Estimate")
```

![glmer2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/535720/d6520548-0bcf-3fef-97df-64655623ca11.png)


# 参考
[easystats/parameters](https://github.com/easystats/parameters)
[Package ‘parameters’](https://cran.r-project.org/web/packages/parameters/parameters.pdf)
[parameters: a powerful and lightweight alternative to broom to describe your models' coefficients](https://easystats.github.io/blog/posts/parameters_presentation/)
[Confidence Intervals (CI)](https://easystats.github.io/parameters/reference/ci.merMod.html)
[cran/MASS](https://github.com/cran/MASS)
[信頼区間 - Wikipedia](https://ja.wikipedia.org/wiki/%E4%BF%A1%E9%A0%BC%E5%8C%BA%E9%96%93)
[数理情報学９ (数理統計学) 5. 最尤推定量](https://ocw.nagoya-u.jp/files/247/05_mle.pdf)
[geom_pointrangeで推定値と95%信頼区間の図示](https://peripatopsis.jimdofree.com/2020/02/26/geom-pointrange%E3%81%A7%E6%8E%A8%E5%AE%9A%E5%80%A4%E3%81%A895-%E4%BF%A1%E9%A0%BC%E5%8C%BA%E9%96%93%E3%81%AE%E5%9B%B3%E7%A4%BA/)