今回はglmでの多重比較を実行する方法をメモしておきます。
すでにいろいろ記事はありますが、自分の学習の記録として残しておきます。
今回はR組み込みデータの中からまず殺虫スプレーの効果のデータを使って試します。
IS<-InsectSprays
head(IS,10)
summary(IS)
> head(IS,10)
count spray
1 10 A
2 7 A
3 20 A
4 14 A
5 14 A
6 12 A
7 10 A
8 23 A
9 17 A
10 20 A
> summary(IS)
count spray
Min. : 0.00 A:12
1st Qu.: 3.00 B:12
Median : 7.00 C:12
Mean : 9.50 D:12
3rd Qu.:14.25 E:12
Max. :26.00 F:12
虫のカウントデータでスプレーは6種類の効果を見ているようです。
早速多重比較をやっていきます。
multcompパッケージをダウンロード。
install.packages("multcomp")
library(multcomp)
カウントデータということで、まずは誤差構造にポアソン分布を仮定してモデル作成し多重比較を実行します。
fit<-glm(count~spray,family=poisson,IS)
multicomp<-glht(fit,mcp(spray="Tukey"))
summary(multicomp)
> summary(multicomp)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glm(formula = count ~ spray, family = poisson, data = IS)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
B - A == 0 0.05588 0.10574 0.528 0.99446
C - A == 0 -1.94018 0.21389 -9.071 < 0.001 ***
D - A == 0 -1.08152 0.15065 -7.179 < 0.001 ***
E - A == 0 -1.42139 0.17192 -8.268 < 0.001 ***
F - A == 0 0.13926 0.10367 1.343 0.74319
C - B == 0 -1.99606 0.21315 -9.364 < 0.001 ***
D - B == 0 -1.13740 0.14961 -7.602 < 0.001 ***
E - B == 0 -1.47727 0.17101 -8.638 < 0.001 ***
F - B == 0 0.08338 0.10215 0.816 0.96090
D - C == 0 0.85866 0.23864 3.598 0.00376 **
E - C == 0 0.51879 0.25261 2.054 0.29105
F - C == 0 2.07944 0.21213 9.803 < 0.001 ***
E - D == 0 -0.33987 0.20189 -1.683 0.51948
F - D == 0 1.22078 0.14815 8.240 < 0.001 ***
F - E == 0 1.56065 0.16973 9.195 < 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
とりあえず結果を出すことはできました。
が、ここで注意点としては、ポアソン分布を仮定したものの、カウントデータのモデルでは過分散になっている可能性が高いことです。
過分散になっていると信頼区間が狭くなり、本当は有意でないパラメータが有意であると判断されてしまう（第1種の過誤）場合が多くなります。
モデルの中身をチェックします。
> summary(fit)
Call:
glm(formula = count ~ spray, family = poisson, data = IS)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.67415 0.07581 35.274 < 2e-16 ***
sprayB 0.05588 0.10574 0.528 0.597
sprayC -1.94018 0.21389 -9.071 < 2e-16 ***
sprayD -1.08152 0.15065 -7.179 7.03e-13 ***
sprayE -1.42139 0.17192 -8.268 < 2e-16 ***
sprayF 0.13926 0.10367 1.343 0.179
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 409.041 on 71 degrees of freedom
Residual deviance: 98.329 on 66 degrees of freedom
AIC: 376.59
Number of Fisher Scoring iterations: 5
仮定した誤差構造が、データの母集団が持つ誤差を正確に表現できているかを検定するためにはDispersion Parameterが規定の値と一致しているかを調べる必要があります。
データが十分に大きい場合は、残差逸脱度（Residual deviance）を残差自由度（degrees of freedom）で割った値でDispersion Parameterを推定できます。
早速計算してみます。
> 98.329/66
[1] 1.489833
大東（2010）によればDispersion Parameterが実際にいくつであれば過分散なのか明確な基準はないそうですが「1.5程度であれば過分散ではないと判断してもかまわないかもしれない」とのことです。微妙ですね。
Dispersion Parameterですが、Pearsonのカイ二乗を残差自由度で割った値の方が正確という説が多いそうです。Pearsonのカイ二乗とは、観測値と推定値の差の二乗を推定値で割ったものの総和です。
計算してみます。
> sum(residuals(fit,type ="pearson")^2)/fit$df.residual
[1] 1.507713
若干値が大きくなりました。
調べてみるとこのPearsonのカイ二乗を使って過分散をチェックできるperformanceパッケージというのがありました。
ダウンロードして実行します。
install.packages("performance")
library(performance)
check_overdispersion(fit)
> check_overdispersion(fit)
# Overdispersion test
dispersion ratio = 1.508
Pearson's Chi-Squared = 99.509
p-value = 0.005
やはり過分散のようです。
負の二項分布でもう一度多重比較をやってみます。
library(MASS)
fit_nb<-glm.nb(count~spray,data=IS)
multicomp2<-glht(fit_nb,mcp(spray="Tukey"))
summary(multicomp2)
> summary(multicomp2)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glm.nb(formula = count ~ spray, data = IS, init.theta = 28.09950714,
link = log)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
B - A == 0 0.05588 0.13082 0.427 0.99805
C - A == 0 -1.94018 0.22733 -8.535 < 0.001 ***
D - A == 0 -1.08152 0.16920 -6.392 < 0.001 ***
E - A == 0 -1.42139 0.18838 -7.545 < 0.001 ***
F - A == 0 0.13926 0.12914 1.078 0.88325
C - B == 0 -1.99606 0.22664 -8.807 < 0.001 ***
D - B == 0 -1.13740 0.16827 -6.759 < 0.001 ***
E - B == 0 -1.47727 0.18755 -7.877 < 0.001 ***
F - B == 0 0.08338 0.12793 0.652 0.98592
D - C == 0 0.85866 0.25076 3.424 0.00744 **
E - C == 0 0.51879 0.26408 1.964 0.34747
F - C == 0 2.07944 0.22568 9.214 < 0.001 ***
E - D == 0 -0.33987 0.21608 -1.573 0.60098
F - D == 0 1.22078 0.16697 7.311 < 0.001 ***
F - E == 0 1.56065 0.18639 8.373 < 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
ほとんど変わりませんが、Ｐ値は若干大きくなっており、より適切な検定になったのではないでしょうか。
念のためもう一度過分散をチェックしてみます。
> check_overdispersion(fit_nb)
# Overdispersion test
dispersion ratio = 1.138
Pearson's Chi-Squared = 75.098
p-value = 0.207
No overdispersion detected.
まあでも元データをみれば最初からだいたいのことはわかりそうなんですよね。
ヒストグラムを描いてみます。
>hist(IS$count)
箱ひげ図も描いてみます。
> boxplot(count~spray,IS)
平均と分散をチェックします。
> mean(IS$count)
[1] 9.5
> var(IS$count)
[1] 51.88732
やっぱりだいぶ分散が大きいですね。ポアソン分布を選択する基準が平均＝分散なのでこれだけ分散が大きいと最初から負の二項分布を仮定してもよさそうです。
説明変数は1個なので一元配置の分散分析（one-way ANOVA）じゃだめなの？という意見もありそうなので、とりあえず誤差の正規性と等分散性もチェックしました。
performanceパッケージでこれらをチェックできます。
fit<-lm(count~spray,data=IS)
check_normality(fit)
check_homogeneity(fit)
> check_normality(fit)
Warning: Non-normality of residuals detected (p = 0.022).
> check_homogeneity(fit)
Warning: Variances differ between groups (Bartlett Test, p = 0.000).
はやりglmで正規分布以外の誤差構造を選択した方がよさそうです。
データが正規分布に従わない場合の、等分散性の検証法としてLevene testというのもあるようです。
> leveneTest(count ~ spray,data=IS)
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 3.8214 0.004223 **
66
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
なんとperformanceパッケージにどの確率分布にすればよいか、ランダムフォレストのアルゴリズムで計算してくれる関数がありました。
実行しようとするとrandomForestパッケージのダウンロードが求められるのでイエス（y）を入力します。
> check_distribution(fit)
# Distribution of Model Family
Predicted Distribution of Residuals
Distribution Probability
normal 44%
tweedie 28%
cauchy 25%
Predicted Distribution of Response
Distribution Probability
neg. binomial (zero-infl.) 69%
beta-binomial 25%
poisson (zero-infl.) 3%
やはり負の二項分布がよいだろうという結果になってますね。
せっかくなのでグラフに描画してみます。
以下の記事が参考になりました。
#結果をデータフレームに格納
alphabet = cld(multicomp2,decreasing=F)
library(tidyverse)
abc = alphabet[["mcletters"]][["Letters"]]
plot =
ggplot(data = IS, aes(x = spray, y = count))+
geom_boxplot()+
stat_summary(geom = 'text', label =abc, fun = max, vjust = -1)+
theme_classic()
print(plot)
～参考文献～
大東健太郎. (2010). 線形モデルから一般化線形モデル (GLM) へ. 雑草研究, 55(4), 268-274.
他にも以下の記事が参考になりました。