【R】glhtとemmeansでTukeyの多重比較

R

RでTukeyの多重比較を行うときの、glht()とemmeans()の使い方を紹介します。

glhtとemmeansの概要

glht（General Linear Hypotheses Test）は、Rのmultcompパッケージに含まれる関数です。一般化線型モデル、線型混合効果モデル、生存モデルなどのパラメトリックモデルの一般線型仮説と多重比較を行います。

emmeans(Estimated marginal means)は、Rのemmeansパッケージに含まれる関数です。線形モデル内の指定された因子または因子の組み合わせの推定限界平均を計算します。また、オプションでそれらの間の比較または対比も計算します。（推定限界平均 は最小二乗平均とも呼ばれます。）

両方ともオプションでTukeyの多重比較を実行できます。

https://cran.r-project.org/web/packages/multcomp/multcomp.pdf

https://cran.r-project.org/web/packages/emmeans/emmeans.pdf

glhtとemmeansの違い

どちらもグループ間の比較を行うツールですが、計算の仕組みやインターフェースに違いがあります。

glht()はz分布を使用するのに対し、emmeans()はt分布を使用します。

z分布とt分布は、サンプルサイズや母集団標準偏差によって使い分けます。

z分布はサンプルサイズが大きい場合（一般的に30≦n）や母集団標準偏差が既知の場合に使われます。t分布はサンプルサイズが小さい時（n<30）や母集団標準偏差が未知の時に使用されます。

https://www.researchgate.net/post/What-are-the-differences-between-these-ways-to-carry-out-posthoc-tests-on-mixed-effects-models
www.researchgate.net
The t-Distribution
The t-distribution describes the standardized distances of sample means to the population mean when the population stand...
www.jmp.com

コード例

TukeyHSD()、glht()、emmeans()を使ってTukeyの多重比較を実行してみます。

次のデータセットと分散分析モデルを使用します。

# サンプルデータ
set.seed(123)
data <- data.frame(
  value = c(value1 <- rnorm(20, mean = 5, sd = 1),
            value2 <- rnorm(20, mean = 5.5, sd = 1),
            value3 <- rnorm(20, mean = 6, sd = 1)),
  group = factor(rep(c("group1", "group2", "group3"), each = 20))
)

# 分散分析モデル
lm <- aov(data$value ~ data$group, data = data)

TukeyHSDの計算結果

TukeyHSD()のコードと計算結果は次のようになります。

TukeyHSD(lm)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = data$value ~ data$group, data = data)

$`data$group`
                   diff         lwr      upr     p adj
group2-group1 0.3071190 -0.39465579 1.008894 0.5468419
group3-group1 0.9648614  0.26308661 1.666636 0.0045735
group3-group2 0.6577424 -0.04403243 1.359517 0.0705975

glhtの計算結果

glht()でp値の出力を得るには、summary()を使います。

コードと計算結果は次のようになります。

library(multcomp)
summary(glht(lm, linfct = mcp("data$group" = "Tukey")))
	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = data$value ~ data$group, data = data)

Linear Hypotheses:
                     Estimate Std. Error t value Pr(>|t|)   
group2 - group1 == 0   0.3071     0.2916   1.053  0.54683   
group3 - group1 == 0   0.9649     0.2916   3.309  0.00462 **
group3 - group2 == 0   0.6577     0.2916   2.255  0.07062 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

glht()では実行毎に小数点の末尾の出力が変わりますが結果の解釈に影響するほどではなさそうです。

Tukeyの他に、multcomp::contrMat()で設定されている比較方法に対応しています。

"Dunnett", "Sequen", "AVE", "Changepoint", "Williams", "Marcus", "McDermott", "UmbrellaWilliams", "GrandMean"

アルファベットの添え字は次のように出力します。

# アルファベットの添え字
cld(glht(lm, linfct = mcp("data$group" = "Tukey")))
group1 group2 group3 
   "a"   "ab"    "b"

emmeansの計算結果

emmeans()でp値の出力を得るには、pairs()を使います。

コードと計算結果は次のようになります。

library(emmeans)
pairs(emmeans(lm, pairwise ~ "group", adjust = "tukey"))
 contrast        estimate    SE df t.ratio p.value
 group1 - group2   -0.307 0.292 57  -1.053  0.5468
 group1 - group3   -0.965 0.292 57  -3.309  0.0046
 group2 - group3   -0.658 0.292 57  -2.255  0.0706

P value adjustment: tukey method for comparing a family of 3 estimates

adjust = "tukey"の他に、emmeans::summary.emmGrid()で設定されているp値の調整方法に対応しています。

"scheffe", "sidak", "bonferroni", "dunnettx", "mvt", "none"

アルファベットの添え字は次のように出力します。

# アルファベットの添え字
cld(emmeans(lm, pairwise ~ "group", adjust = "tukey"), Letters = letters)
 group  emmean    SE df lower.CL upper.CL .group
 group1   5.14 0.206 57     4.73     5.55  a    
 group2   5.45 0.206 57     5.04     5.86  ab   
 group3   6.11 0.206 57     5.69     6.52   b   

Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

デフォルトではアルファベットの代わりに数字が出力されるので、cld()の引数"Letters = letters"を追加することでアルファベット出力しています。

結果の比較

TukeyHSD()、glht()、emmeans()それぞれ同様の計算結果が得られました。

TukeyHSDglhtemmeans
group1 - group20.54684190.546830.5468
group1 - group30.00457350.004620.0046
group2 - group30.07059750.070620.0706

まとめ

TukeyHSD()はaov()による分散分析モデルで多重比較を簡単に実行できる関数ですが、cld()を使ったアルファベットの添え字出力に対応していません。

glht()はTukeyの多重比較だけでなく、multcomp::contrMat()で設定されている比較方法に対応できる柔軟性があります。また、cld()を使ったアルファベットの添え字出力に対応していますが、交互作用を含む複雑なモデルでアルファベットの添え字がおかしくなることがあります。

glht and emmeans returning crazy compact letters for unbalanced dataset in R
HiI am having a problem when trying to get glht or emmeans to define compact letters for a dataset with unequal sample s...
stackoverflow.com

emmeans()は、非常に多くのモデルに対応している上、emmeans::summary.emmGrid()に設定されているp値の調整方法に対応できるため、より柔軟性の高い分析ができます。ただし、出力される平均値と標準誤差は推定平均に基づいているため注意が必要です。

特徴TukeyHSDglhtemmeans
対応モデルaovlm, glm, aovlm, glm, aovなど多数
柔軟性低い高い高い
用途単純なANOVAカスタマイズ可能な比較モデル条件を考慮した比較
cld()対応なしありあり

それぞれの関数を使う際、分析の目的やモデルの複雑さに応じて選択するとよいでしょう。

今回はglhtとemmeansを使ってTukeyの多重比較を実行しました。

少しでもあなたの解析が楽になりますように！

タイトルとURLをコピーしました