モデルの品質（model quality）のチェック

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

これらの関数による解析結果の妥当性を評価するために、タイトルに上げたようないくつかの事項を確認したい場合があります。
AICや対数尤度を確認するのとは異なり、これらを確認するにはオブジェクトからデータを抽出したりパッケージを使い分けたりする必要があります。よって結構面倒です。

しかしeasystatsパッケージ群に含まれるperformanceパッケージを使うとかなり簡単にこれらの事項を確認できます。ここではその一例を紹介したいと思います。

※　指標や検定について詳しい説明はしません（統計は勉強中なので踏み込んだ説明が出来ません。ごめんなさい）。あくまでパッケージの紹介と手順の確認のためのメモとして見てもらえればと思います。
※　指標の根拠や算出方法については、performanceパッケージのマニュアル（あるいはレファレンス）やGitHubにあるソースコードで確認できます。レポート等に使用されたい方はそちらを参照してください。

インストール

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

CRAN

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

install.packages("performance")

GitHub

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

remotes::install_github("easystats/easystats")

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

各事項の確認

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

library(performance)
library(MASS)
library(lme4)
library(glmmTMB)

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

R2の確認

performanceパッケージでは総称的関数であるr2で各モデルに対する適切な決定係数（あるいは疑似決定係数）を計算できます。

各モデルに対応するR2は以下のようになっています：

線形モデル：R2および調整済みR2
一般化線形モデル：Nagelkerke’s R2
ロジスティックモデル：Tjur’s R2
混合モデル：Nakagawa’s R2

参照：Compute the model's R2

lm

> rm1 <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
> r2(rm1)
# R2 for Linear Regression

       R2: 0.835
  adj. R2: 0.811

glm

> rm2 <- glm(Days ~ Sex/(Age + Eth*Lrn), data = quine, family = "poisson")
> r2(rm2)
$R2_Nagelkerke
Nagelkerke's R2 
      0.9905114

> # ロジスティックモデル
> rm3 <- glm(cbind(incidence, size - incidence) ~ period, family = binomial, data = cbpp)
> r2(rm3)
$R2_Tjur
  Tjur's R2 
0.001039836

glmer

> # 一般化線形混合モデル
> rm4 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), family = binomial, data = cbpp)
> r2(rm4)
# R2 for Mixed Models

  Conditional R2: 0.187
     Marginal R2: 0.085

多重共線性の確認

performanceパッケージでは多重共線性の指標の一つであるVIF統計量をcheck_collinearityで確認できます。
各パラメータに対するVIFの値が示されると同時に、棒グラフで値が可視化できます。VIF<5の場合低い相関、5≤VIF≤10の場合中程度の相関、10<VIFの場合高い相関とみなされます。

参照：Check for multicollinearity of model terms

glm

> cm1 <- glm(mpg ~ wt + cyl + gear + disp, data = mtcars, family = "gaussian")
> co1 <- check_collinearity(cm1)
> co1
# 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(co1)

glmer

> form <- TICKS ~ YEAR + HEIGHT + (1|BROOD) + (1|INDEX) + (1|LOCATION)
> cm2 <- glmer(form, family = "poisson", data = grouseticks)
 警告メッセージ: 
 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?
> co2 <- check_collinearity(cm2)
> co2
# Check for Multicollinearity

Low Correlation

 Parameter  VIF Increased SE
      YEAR 1.00         1.00
    HEIGHT 1.00         1.00
> plot(co2)

正規性の確認

performanceパッケージでは残差の正規性をcheck_normalityで確認できます。
標準化された残差に対してシャピロ検定が実行されます。lmあるいはglmの結果に対してのみ使えます。

参照：Check model for (non-)normality of residuals.

> nm <- glm(mpg ~ wt + cyl + gear + disp, data = mtcars, family = "gaussian")
> no <- check_normality(nm)
OK: Residuals appear as normally distributed (p = 0.230).
> plot(no)

過分散の確認

performanceパッケージではcheck_overdispersionで過分散を確認できます。
過分散の指標としてピアソン残差の2乗和に対するカイ二乗検定の結果が示されます。dispersion rateは期待値との比、Pearson's Chi-Squaredはピアソン残差の2乗和、p-valueはp値です。確率分布がポアソン分布になっているモデルが対象です。
サンプル数（正確には残差自由度）が十分大きい場合に、あくまで過分散のラフな指標として使用可能とのことです（参考）。

ピアソン残差：pr_i = \frac{(y_i - \hat{\mu}_i)}{\sqrt{\hat{V}_i}}\quad  (y_i:観測iにおける観測値;\quad \hat{\mu}_i:観測iにおける予測値;\quad \hat{V}_i:観測iにおける予測値の分散)

（参考）

※　ピアソン残差の2乗和が自由度＝残差自由度のカイ二乗分布に従うことを利用しているのですが、上記の参考を見てもいまいち導出過程がわかりませんでした。

参照：Check overdispersion of GL(M)M's

glm

> po1 <- glm(Days ~ Sex/(Age + Eth*Lrn), data = quine, family = "poisson")
> check_overdispersion(po1)
# Overdispersion test

       dispersion ratio =   10.534
  Pearson's Chi-Squared = 1390.436
                p-value =  < 0.001

Overdispersion detected.

glmer

> form <- TICKS ~ YEAR + HEIGHT + (1|BROOD) + (1|INDEX) + (1|LOCATION)
> po2 <- glmer(form, family = "poisson", data = grouseticks)
 警告メッセージ: 
 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?
> check_overdispersion(po2)
# Overdispersion test

       dispersion ratio =   0.412
  Pearson's Chi-Squared = 163.211
                p-value =       1

No overdispersion detected.

ゼロ過剰の確認

performanceパッケージではcheck_zeloinflationでモデルのゼロ過剰を確認できます。
実際の0の数とモデルから期待される0の数、およびその比が示されます。こちらも確率分布がポアソン分布になっているモデルが対象です。

参照：Check for zero-inflation in count models

glm

> zi1 <- glm(count ~ spp + mined, family = "poisson", data = Salamanders)
> check_zeroinflation(zi1)
# Check for zero-inflation

   Observed zeros: 387
  Predicted zeros: 298
            Ratio: 0.77

Model is underfitting zeros (probable zero-inflation).

glmer

> zi2 <- glmer(count ~ spp + mined + (1|site), family = "poisson", data = Salamanders)
> check_zeroinflation(zi2)
# Check for zero-inflation

   Observed zeros: 387
  Predicted zeros: 311
            Ratio: 0.80

Model is underfitting zeros (probable zero-inflation).

備考

自己相関を確認するcheck_autocorrelationもあるのですが、これはデータに順序があると仮定してダービーワトソン検定を行うもののようです。残差がAR(1) かどうかを確認するためのものなので注意が必要です。

参照：Check model for independence of residuals.

performanceパッケージはglmmTMB等の他のパッケージの関数にも対応しています。
また収束や外れ値などの確認もできるので気になる方は確認してみてください。

glmやglmerのR2、多重共線性、正規性、過分散、ゼロ過剰を確認する (easystats)

モデルの品質（model quality）のチェック

インストール

各事項の確認

R2の確認

lm

glm

glmer

多重共線性の確認

glm

glmer

正規性の確認

過分散の確認

glm

glmer

ゼロ過剰の確認

glm

glmer

備考

参考