More than 3 years have passed since last update.

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

Last updated at 2020-05-25Posted at 2020-05-25

解析結果の視覚的な確認

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

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

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

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

インストール

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

CRAN

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

install.packages("parameters")

GitHub

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

remotes::install_github("easystats/easystats")

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

`parameters::model_parameters`について

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

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

> 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

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

> 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を使わない方法で信頼区間を算出しているようです（参考）。よって、「Estimate ± 1.959964... x SE」というよく見る計算式による計算結果とは必ずしも一致しません。対してparameters::ci_wald(dof = Inf)による計算結果は「Estimate ± 1.959964... x SE」による計算結果と一致します（参考。厳密にはt分布が使用されていますが自由度無限大なので正規分布のquantileと一致しています）。

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

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

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

推定値の視覚化

ここでは以下のパッケージを使用します。

library(parameters)
library(tidyverse)
library(lme4)
library(glmmTMB)

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

基本的には、パラメータ推定→parameters::model_parametersの実行→ggplot2で描画、という流れで推定値と信頼区間をプロットします。
推定値と信頼区間はggplot2::geom_pointrangeというとても便利な関数を使って描画します（こちらのページで紹介されています）。またエラーバーを描画するgeom_*関数は過去の記事で列挙してあります。

glm

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

> 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_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")

描画の参考：
ggplot2で縦軸と横軸をひっくり返したい
 Reverse the order of a categorical axis in ggplot2

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

> 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")

描画の参考：
Nudge points a fixed distance - ggplot2

glmer

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

> 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")

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

> 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")