CRAN: https://cran.r-project.org/web/packages/performance/index.html rdrr.io: https://rdrr.io/cran/performance/
YouTubeでの解説(英語)
コードの使用法
https://rdrr.io/cran/performance/f/README.md
回帰モデルを構築するパッケージモデルの品質と適合度の指標を計算するためにR2（二乗）、RMSE（二乗平均平方根誤差）、クラス内相関係数（ICC）などの指標や（混合）モデルの過分散、ゼロインフレ、収束、特異性などをチェックする関数が含まれている。
library(performance)
線形回帰 R二乗
model <- lm(mpg ~ wt + cyl, data = mtcars) r2(model)
Show in New Window
# R2 for Linear Regression
R2: 0.830
adj. R2: 0.819
二項ロジット 疑似決定係数
model <- glm(am ~ wt + cyl, data = mtcars, family = binomial) r2(model)
Tjurの読み方はおそらくテュアー。
# R2 for Logistic Regression Tjur's R2: 0.705
文献はこちら。
https://www.tandfonline.com/doi/abs/10.1198/tast.2009.08210
- Tjur, T. (2009). Coefficients of determination in logistic regression models - A new proposal: The coefficient of discrimination. The American Statistician, 63(4), 366-372.
R二乗がどのくらい必要なのかは疑問だが、NagelkerkeのR2やCox. & shellのR2よりはよいかもしれない。McFadden'sやEfron's R2もある。
順序ロジット 疑似決定係数
library(MASS) data(housing) model <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) r2(model)
結果。
Nagelkerke's R2: 0.108
また、r2_bayes()、r2_coxsnell()、r2_nagelkerke()などの関数により、指定したいR2乗を直接算出できる。
混合モデルの場合、条件付きR2乗とマージナルR2乗が返される。マージナルR2乗は、固定効果の分散のみを考慮し、モデルの分散が固定効果部分のみによってどれだけ説明されるかを示すものである。条件付きR2乗は、固定効果とランダム効果の両方を考慮し、モデルの分散のどれくらいが「完全な」モデルによって説明されるかを示す。
frequentistの混合モデルでは、r2() (resp. r2_nakagawa()) は平均ランダム効果分散を計算するので、r2()はランダムスロープや入れ子ランダム効果など、より複雑なランダム効果構造を持つ混合モデルにも適している (Johnson 2014; Nakagawa, Johnson, and Schielzeth 2017)。
ベイズ
set.seed(123) library(rstanarm) model <- stan_glmer(Petal.Length ~ Petal.Width + (1 | Species), data = iris, cores = 4) r2(model)
結果。
# Bayesian R2 with Compatibility Interval
Conditional R2: 0.953 (95% CI [0.940, 0.963])
Marginal R2: 0.824 (95% CI [0.725, 0.898])
線形混合モデル
library(lme4) model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy) r2(model)
結果。
# R2 for Mixed Models
Conditional R2: 0.799
Marginal R2: 0.279
級内相関係数(ICC)
R二乗と同様に、ICCは説明される分散に関する情報を提供し、「母集団におけるグループ化構造によって説明される分散の割合」として解釈することができる（Hox 2010）。
icc() は、'stanreg'モデルを含む様々な混合モデルオブジェクトの ICC を計算する。
library(lme4) model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy) icc(model)
結果。
# Intraclass Correlation Coefficient
Adjusted ICC: 0.722
Conditional ICC: 0.521
brmsfitクラスのモデルは、、、
library(brms) set.seed(123) model <- brm(mpg ~ wt + (1 | cyl) + (1 + wt | gear), data = mtcars) icc(model)
# Intraclass Correlation Coefficient
Adjusted ICC: 0.930
Conditional ICC: 0.771
モデル診断
過大分散のチェック。 過大分散は、観測されたデータの分散が、モデルの仮定から期待される分散（ポアソンの場合、分散は結果の平均にほぼ等しい）よりも大きい場合に発生する。 check_overdispersion() は、カウントモデル（混合モデルを含む）が過大分散であるかどうかをチェックする。
library(glmmTMB) data(Salamanders) model <- glm(count ~ spp + mined, family = poisson, data = Salamanders) check_overdispersion(model)
結果。
# Overdispersion test
dispersion ratio = 2.946
Pearson's Chi-Squared = 1873.710
p-value = < 0.001
Overdispersion detected.
過剰分散は，分散パラメータをモデル化するか（すべてのパッケージで可能ではない），別の分布族を選択することで修正できる（準ポアソンや負の二項分布など， (Gelman and Hill 2007)を参照）．
ゼロインフレのチェック
ゼロインフレは（準）ポアソンモデルにおいて、観測されたゼロの量が予測されたゼロの量より多い場合に示され、そのためモデルはゼロをアンダーフィットしていることになる。このような場合、負の二項モデルまたはゼロインフレートモデルを使用することが推奨される。
適合したモデルにゼロインフレートがあるかどうかを調べるには check_zeroinflation() を使用する。
model <- glm(count ~ spp + mined, family = poisson, data = Salamanders) check_zeroinflation(model)
結果。
# Check for zero-inflation
Observed zeros: 387
Predicted zeros: 298
Ratio: 0.77
特異なモデル適合をチェックする
特異なモデル適合とは、分散共分散行列のいくつかの次元が正確にゼロとして推定されていることを意味する。これは過度に複雑なランダム効果構造を持つ混合モデルでよく発生する。
check_singularity() は、混合モデル（lme, merMod, glmmTMB または MixMod クラス）の特異性をチェックし、モデルの適合が特異である場合に TRUE を返す。
# library
library(lme4)
data(sleepstudy)
# prepare data
set.seed(123)
sleepstudy$mygrp <- sample(1:5, size = 180, replace = TRUE)
sleepstudy$mysubgrp <- NA
for (i in 1:5) {
filter_group <- sleepstudy$mygrp == i
sleepstudy$mysubgrp[filter_group] <- sample(1:30, size = sum(filter_group), replace = TRUE)
}
# fit strange model model <- lmer(Reaction ~ Days + (1 | mygrp/mysubgrp) + (1 | Subject), data = sleepstudy) check_singularity(model)
結果。
boundary (singular) fit: see ?isSingular [1] TRUE
シンギュラーフィットの問題を解決するための解決法は、こちら。 https://easystats.github.io/performance/reference/check_singularity.html
異種分散性のチェック
線形モデルは、一定の誤差分散（homoskedasticity）を仮定している。check_heteroscedasticity()関数は、この仮定が破られているかどうかを評価する。
data(cars) model <- lm(dist ~ speed, data = cars) check_heteroscedasticity(model)
結果。
Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.031).
モデルチェックの包括的な可視化
performanceには、check_collinearity()、check_normality()、check_heteroscedasticity()など、モデルの仮定をチェックする関数が多数用意されている。包括的なチェックを行うには、 check_model() を使用する。
# defining a model model <- lm(mpg ~ wt + am + gear + vs * cyl, data = mtcars) # checking model assumptions check_model(model)
モデル性能の要約
model_performance() は、回帰モデルのモデル性能の指標を計算する。モデルオブジェクトに依存するが、典型的な指標は r-2乗、AIC、BIC、RMSE、ICC または LOOIC であろう。
線形回帰
m1 <- lm(mpg ~ wt + cyl, data = mtcars) model_performance(m1)
結果。
# Indices of model performance AIC | BIC | R2 | R2 (adj.) | RMSE | Sigma ----------------------------------------------------- 156.010 | 161.873 | 0.830 | 0.819 | 2.444 | 2.568
ロジスティック回帰
m2 <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial") model_performance(m2)
結果。
# Indices of model performance AIC | BIC | Tjur's R2 | RMSE | Sigma | Log_loss | Score_log | Score_spherical | PCP -------------------------------------------------------------------------------------------- 31.298 | 35.695 | 0.478 | 0.359 | 0.934 | 0.395 | -14.903 | 0.095 | 0.743
線形混合モデル
library(lme4) m3 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy) model_performance(m3)
結果。
# Indices of model performance AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma ---------------------------------------------------------------------------------- 1755.628 | 1756.114 | 1774.786 | 0.799 | 0.279 | 0.722 | 23.438 | 25.592
モデル比較
compare_performance()関数は、複数のモデル（異なるタイプのモデルを含む）の性能と品質を比較するために使用することができる。
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12) outcome <- gl(3, 1, 9) treatment <- gl(3, 3) m4 <- glm(counts ~ outcome + treatment, family = poisson()) compare_performance(m1, m2, m3, m4)
結果。
Name | Model | AIC | AIC weights | BIC | BIC weights | RMSE | Sigma | Score_log | Score_spherical | R2 | R2 (adj.) | Tjur's R2 | Log_loss | PCP | AICc | AICc weights | R2 (cond.) | R2 (marg.) | ICC | Nagelkerke's R2 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- m1 | lm | 156.010 | < 0.001 | 161.873 | < 0.001 | 2.444 | 2.568 | | | 0.830 | 0.819 | | | | | | | | | m2 | glm | 31.298 | 1.000 | 35.695 | 1.000 | 0.359 | 0.934 | -14.903 | 0.095 | | | 0.478 | 0.395 | 0.743 | | | | | | m3 | lmerMod | 1763.986 | < 0.001 | 1783.144 | < 0.001 | 23.438 | 25.592 | | | | | | | | 1764.471 | | 0.799 | 0.279 | 0.722 | m4 | glm | 56.761 | < 0.001 | 57.747 | < 0.001 | 3.043 | 1.132 | -2.598 | 0.324 | | | | | | | | | | | 0.657
モデル性能の一般的な指標
また、モデルの性能を表す合成指標を簡単に計算し、最適なモデルから悪いモデルへと並べ替えることもできる。
compare_performance(m1, m2, m3, m4, rank = TRUE)
結果。
# Comparison of Model Performance Indices Name | Model | RMSE | Sigma | AIC weights | BIC weights | Performance-Score -------------------------------------------------------------------------------- m2 | glm | 0.359 | 0.934 | 1.000 | 1.000 | 100.00% m4 | glm | 3.043 | 1.132 | < 0.001 | < 0.001 | 46.89% m1 | lm | 2.444 | 2.568 | < 0.001 | < 0.001 | 46.09% m3 | lmerMod | 23.438 | 25.592 | < 0.001 | < 0.001 | 0.00%
モデル性能の指標の視覚化
最後に可視化機能を提供する。
plot(compare_performance(m1, m2, m4, rank = TRUE))
テストモデル
test_performance() (およびそのベイズ版である test_bf) は、入力に基づき、最も関連性が高く適切なテストを実行する。(たとえば、モデルがネストされているかどうかなど)。
データ。
set.seed(123) data(iris) lm1 <- lm(Sepal.Length ~ Species, data = iris) lm2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris) lm3 <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris) lm4 <- lm(Sepal.Length ~ Species * Sepal.Width + Petal.Length + Petal.Width, data = iris)
test_performance(lm1, lm2, lm3, lm4)
結果。
Name | Model | BF | Omega2 | p (Omega2) | LR | p (LR) ------------------------------------------------------------ lm1 | lm | | | | | lm2 | lm | > 1000 | 0.69 | < .001 | -6.25 | < .001 lm3 | lm | > 1000 | 0.36 | < .001 | -3.44 | < .001 lm4 | lm | > 1000 | 0.73 | < .001 | -7.77 | < .001
test_bf(lm1, lm2, lm3, lm4)
結果。
Bayes Factors for Model Comparison
Model BF
[lm2] Species + Petal.Length 3.45e+26
[lm3] Species * Sepal.Width 4.69e+07
[lm4] Species * Sepal.Width + Petal.Length + Petal.Width 7.58e+29
* Against Denominator: [lm1] Species
* Bayes Factor Type: BIC approximation
