統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

反復測定分散分析のGreenhouse-GeisserやHuynh-Feldtとは何か?

Toukei Consul Banner

KH Coder Consul Banner

繰り返し測定分散分析に登場するGreenhouse-Geisserの方法とかHuynh-Feldtの方法とかは何をしているのか?

一言で言えば、同じ対象者を複数回測定するために起きる相関を考慮した検定方法で、自由度を調整する方法だ。

統計ソフトで計算してしまうと何をしているのか不明のまま使うことになり、ちょっと心もとない。

One by oneで解説しているPDFを見つけたので、one by oneでやってみる。ANOVA君の結果と照らし合わせながら確認する。

今回の教科書的PDF

Herv´ e Abdi The Greenhouse-Geisser Correction. In Neil Salkind (Ed.), Encyclopedia of Research Design. Thousand Oaks, CA: Sage. 2010.

http://www.utd.edu/~herve/abdi-GreenhouseGeisser2010-pretty.pdf

データの準備

まずデータを準備する。要因はAの一つだけで、Aのグループは4グループある。4回繰り返して測定しているという意味。対象者は5名。

a1 <- c(76,60,58,46,30)
a2 <- c(64,48,34,46,18)
a3 <- c(34,46,32,32,36)
a4 <- c(26,30,28,28,28)

dat <- data.frame(a1,a2,a3,a4)
matrix(summary(dat)[4,])
dat
(A.bar <- matrix(c(mean(a1),mean(a2),mean(a3),mean(a4))))
(s.bar <- matrix(c(
 mean(unlist(c(dat[1,]))),
 mean(unlist(c(dat[2,]))),
 mean(unlist(c(dat[3,]))),
 mean(unlist(c(dat[4,]))),
 mean(unlist(c(dat[5,])))
)))
matrix(c(sd(a1),sd(a2),sd(a3),sd(a4)))
matrix(c(var(a1),var(a2),var(a3),var(a4)))

データ一覧、要因Aのグループごとの平均、対象者ごとの平均、要因Aグループごとの標準偏差、分散を計算してみると、以下の通り。

> dat
  a1 a2 a3 a4
1 76 64 34 26
2 60 48 46 30
3 58 34 32 28
4 46 46 32 28
5 30 18 36 28

> (A.bar <- matrix(c(mean(a1),mean(a2),mean(a3),mean(a4))))
     [,1]
[1,]   54
[2,]   42
[3,]   36
[4,]   28

> (s.bar <- matrix(c(
+  mean(unlist(c(dat[1,]))),
+  mean(unlist(c(dat[2,]))),
+  mean(unlist(c(dat[3,]))),
+  mean(unlist(c(dat[4,]))),
+  mean(unlist(c(dat[5,])))
+ )))
     [,1]
[1,]   50
[2,]   46
[3,]   38
[4,]   38
[5,]   28

> matrix(c(sd(a1),sd(a2),sd(a3),sd(a4)))
          [,1]
[1,] 17.146428
[2,] 17.146428
[3,]  5.830952
[4,]  1.414214

> matrix(c(var(a1),var(a2),var(a3),var(a4)))
     [,1]
[1,]  294
[2,]  294
[3,]   34
[4,]    2

ANOVA君で実行してみる

こちらから最新のANOVA君をダウンロードしてお借りする。

riseki.php.xdomain.jp

source("C:\\Users\\happy\\Downloads\\anovakun_482.txt")

anovakun(dat, "sA", 4, nopost=TRUE)

単純な反復測定分散分析は sA と表現する。ちなみに、二要因混合計画 Split Plot Design は AsBと書く。詳しくは以下を参照。

toukeier.hatenablog.com

結果は以下の通り。DESCRIPTIVE STATISTICSは、上記の値と一致した計算結果が表示されている。今回 one by one でトライするのは ANOVA TABLE の結果を得ようというものだ。

> anovakun(dat, "sA", 4, nopost=TRUE)

[ sA-Type Design ]

This output was generated by anovakun 4.8.2 under R version 3.5.1.
It was executed on Sun Feb 10 17:36:33 2019.

 
<< DESCRIPTIVE STATISTICS >>

----------------------------
  A   n     Mean     S.D. 
----------------------------
 a1   5  54.0000  17.1464 
 a2   5  42.0000  17.1464 
 a3   5  36.0000   5.8310 
 a4   5  28.0000   1.4142 
----------------------------


<< SPHERICITY INDICES >>

== Mendoza's Multisample Sphericity Test and Epsilons ==

-------------------------------------------------------------------------
 Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
-------------------------------------------------------------------------
      A  0.0078      6.5990   5 0.2786 ns  0.3333 0.4460 0.5872 0.3333 
-------------------------------------------------------------------------
                              LB = lower.bound, GG = Greenhouse-Geisser
                             HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller


<< ANOVA TABLE >>

---------------------------------------------------------
 Source         SS  df        MS  F-ratio  p-value     
---------------------------------------------------------
      s  1152.0000   4  288.0000                       
---------------------------------------------------------
      A  1800.0000   3  600.0000   5.3571   0.0142 *   
  s x A  1344.0000  12  112.0000                       
---------------------------------------------------------
  Total  4296.0000  19  226.1053                       
             +p < .10, *p < .05, **p < .01, ***p < .001


output is over --------------------///

まず、自由度調整をしない素の計算を one by oneでやってみる

全平方和(SS)、要因Aの平方和(SSA)、対象者内平方和(SSs)を計算して、要因A×対象者内平方和(SSsA)を計算する。SSsAのバラツキを超えて、SSAが大きければ、Aは意味がある(異なる)と言えるわけだ。

X <- c(a1,a2,a3,a4)
(X.bar <- mean(X))
mean(unlist(dat))

(SS <- sum((X-X.bar)^2))

(SSA <- length(a1)*sum((A.bar - X.bar)^2))
a <- length(dat)
a-1

(SSs <- length(dat[1,])*sum((s.bar - X.bar)^2))
s <- length(a1)
s-1

(SSsA <- SS-SSA-SSs)
(a-1)*(s-1)

(F <- SSA/(a-1)/(SSsA/((a-1)*(s-1))))
(pF <- pf(F, a-1, (a-1)*(s-1), lower.tail=FALSE))

結果は、以下の通り。ANOVA君の出力における、ANOVA TABLE内の SS (Totalは4296, sは1152, Aは1800, sxAは1344), df (sは4, Aは3, s x Aは12), F-ratio (5.3571), p-value (0.0142) に一致した値が見て取れる。

> (SS <- sum((X-X.bar)^2))
[1] 4296
> 
> (SSA <- length(a1)*sum((A.bar - X.bar)^2))
[1] 1800
> a <- length(dat)
> a-1
[1] 3
> 
> (SSs <- length(dat[1,])*sum((s.bar - X.bar)^2))
[1] 1152
> s <- length(a1)
> s-1
[1] 4
> 
> (SSsA <- SS-SSA-SSs)
[1] 1344
> (a-1)*(s-1)
[1] 12
> 
> (F <- SSA/(a-1)/(SSsA/((a-1)*(s-1))))
[1] 5.357143
> (pF <- pf(F, a-1, (a-1)*(s-1), lower.tail=FALSE))
[1] 0.01423304

自由度調整をやってみる

上記の計算だと、分母のバラツキは、同一対象が繰り返し測定していることが考慮されていない。同一対象者では同種のテストや検査を繰り返し行えば、似たような値になるのはよく知られている。それを考慮して自由度調整を行うのがGreenhouse-Geisserの方法とHuynh-Feldtの方法だ。

Greenhouse-Geisserの方法

まずANOVA君の結果

gg=TRUEと付けると、Greenhouse-Geisserの方法でANOVA TABLEを出力してくれる。

anovakun(dat, "sA", 4, nopost=TRUE, gg=TRUE)

ANOVA TABLEの結果はこちら。p valueが0.0142から0.0600に上昇しているのがわかる。dfも少数点以下があって、小さい値だ。Greenhouse-Geisserの方法はとてもConservative(控えめ)な方法で、結果は統計学的有意ではなくなった。

<< ANOVA TABLE >>

== Adjusted by Greenhouse-Geisser's Epsilon ==

-----------------------------------------------------------
 Source         SS   df         MS  F-ratio  p-value     
-----------------------------------------------------------
      s  1152.0000    4   288.0000                       
-----------------------------------------------------------
      A  1800.0000 1.34  1345.4082   5.3571   0.0600 +   
  s x A  1344.0000 5.35   251.1429                       
-----------------------------------------------------------
  Total  4296.0000   19   226.1053                       
               +p < .10, *p < .05, **p < .01, ***p < .001
One by oneでやってみる

やり方は、分散共分散行列を作って \epsilon を推定する。母集団分散共分散行列(Boxの方法で使用。今回は割愛)からサンプル分散共分散行列(Greenhouse-Geisserの方法で使用)を作るという方法を使う。

A1 <- cbind(a1-mean(a1),a2-mean(a2),a3-mean(a3),a4-mean(a4))
(dat.AA <- t(A1)%*%A1/(length(a1)-1))

(t.bar <- mean(unlist(dat.AA)))
(ta.bar <- colMeans(dat.AA))
(t.bar.diff <- ta.bar-t.bar)

(mat.ta.bar <- matrix(rep(ta.bar,4),4,4))
(t.mat.ta.bar <- t(matrix(rep(ta.bar,4),4,4)))

(s.AA <- (dat.AA-t.bar)-(mat.ta.bar-t.bar)-(t.mat.ta.bar-t.bar))

(epsilon.hat <- sum(diag(s.AA))^2/((a-1)*sum(s.AA^2)))

epsilon.hat*(a-1)
epsilon.hat*(a-1)*(s-1)
(pF.GG <- pf(F, epsilon.hat*(a-1), epsilon.hat*(a-1)*(s-1), lower.tail=FALSE))

以下のdat.AAが母集団の分散共分散行列で、s.AAがサンプルの分散共分散行列だ。epsilon.hatがGreenhouse-Geisserの方法の \epsilonの推定値。F値の分子の自由度(epsilon.hat x (a-1)=1.337884)と分母の自由度(epsilon.hat x (a-1) x (s-1)=5.351536)及びp値(pF.GG=0.05999477)は、上記のANOVA君の結果と同一だ(ANOVA君の出力ではそれぞれ、1.34, 5,35, 0.0600とある)。

> (dat.AA <- t(A1)%*%A1/(length(a1)-1))
     [,1] [,2] [,3] [,4]
[1,]  294  258    8   -8
[2,]  258  294    8   -8
[3,]    8    8   34    6
[4,]   -8   -8    6    2

> (s.AA <- (dat.AA-t.bar)-(mat.ta.bar-t.bar)-(t.mat.ta.bar-t.bar))
     [,1] [,2] [,3] [,4]
[1,]   90   54  -72  -72
[2,]   54   90  -72  -72
[3,]  -72  -72   78   66
[4,]  -72  -72   66   78

> (epsilon.hat <- sum(diag(s.AA))^2/((a-1)*sum(s.AA^2)))
[1] 0.4459613

> epsilon.hat*(a-1)
[1] 1.337884

> epsilon.hat*(a-1)*(s-1)
[1] 5.351536

> (pF.GG <- pf(F, epsilon.hat*(a-1), epsilon.hat*(a-1)*(s-1), lower.tail=FALSE))
[1] 0.05999477

ちなみに、 \epsilon の推定値 0.4459613 は、ANOVA君のSPHERICITY INDICIESの項にもGGの欄に同じ数値 0.4460 が出ている。

<< SPHERICITY INDICES >>

== Mendoza's Multisample Sphericity Test and Epsilons ==

-------------------------------------------------------------------------
 Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
-------------------------------------------------------------------------
      A  0.0078      6.5990   5 0.2786 ns  0.3333 0.4460 0.5872 0.3333 
-------------------------------------------------------------------------
                              LB = lower.bound, GG = Greenhouse-Geisser
                             HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller

Huynh-Feldtの方法

まずANOVA君の結果

今度はHuynh-Feldtの方法で自由度調整を行ってみる。

anovakun(dat, "sA", 4, nopost=TRUE, hf=TRUE)

ANOVA君の結果は以下の通り。p値は0.0411で統計学的有意。

<< ANOVA TABLE >>

== Adjusted by Huynh-Feldt-Lecoutre's Epsilon ==

-----------------------------------------------------------
 Source         SS   df         MS  F-ratio  p-value     
-----------------------------------------------------------
      s  1152.0000    4   288.0000                       
-----------------------------------------------------------
      A  1800.0000 1.76  1021.8341   5.3571   0.0411 *   
  s x A  1344.0000 7.05   190.7424                       
-----------------------------------------------------------
  Total  4296.0000   19   226.1053                       
               +p < .10, *p < .05, **p < .01, ***p < .001
One by oneでやってみる

Huynh-Feldtの方法は、Greenhouse-Geisserの方法の \epsilonの推定値( \hat{\epsilon})を使ってより適切に修正した方法なのだ。Huynh-Feldtの \tilde{\epsilon}は、 \hat{\epsilon}よりも大きくなり、統計学的有意に出やすい。事実、今回はGreenhouse-Geisserの方法では有意でなかったが、Huynh-Feldtの方法だと有意になった。

(epsilon.tilde <- (s*(a-1)*epsilon.hat-2)/((a-1)*(s-1-(a-1)*epsilon.hat)))
epsilon.tilde*(a-1)
epsilon.tilde*(a-1)*(s-1)
(pF.HF <- pf(F, epsilon.tilde*(a-1), epsilon.tilde*(a-1)*(s-1), lower.tail=FALSE))

以下の結果は、ANOVA君の結果と一致している。F値の分子の自由度(epsilon.tilde x (a-1)=1.761538)、分母の自由度(epsilon.tilde x (a-1) x (s-1)=7.046154)、p値(pF.HF=0.04114125)はそれぞれ、ANOVA君の結果では、1.76, 7.05, 0.0411 と出力されている。

> (epsilon.tilde <- (s*(a-1)*epsilon.hat-2)/((a-1)*(s-1-(a-1)*epsilon.hat)))
[1] 0.5871795

> epsilon.tilde*(a-1)
[1] 1.761538

> epsilon.tilde*(a-1)*(s-1)
[1] 7.046154

> (pF.HF <- pf(F, epsilon.tilde*(a-1), epsilon.tilde*(a-1)*(s-1), lower.tail=FALSE))
[1] 0.04114125

Huynh-Feldt-Lecoutreの方法(要因数が2以上の場合Split Plot Designなど)

Huynh-Feldtの方法は、要因が2つ以上の場合、例えばSplit Plot Design(二要因混合計画「AsB」)などの場合は適切ではないことがLecoutreによって指摘された。修正の計算スクリプトは以下の通り。ANOVA君においては、この点は修正済み(出力にHF = Huynh-Feldt-Lecoutreとの記載あるので確認できる)。

g <- 1 #要因の数今回は要因Aのみで1
(epsilon.tilde.lecoutre <- ((s-g+1)*(a-1)*epsilon.hat-2)/((a-1)*(s-g-(a-1)*epsilon.hat)))

今回のように要因が1つの場合は結果は一致する。上記のepsilon.tildeと同じことがわかる。

> (epsilon.tilde.lecoutre <- ((s-g+1)*(a-1)*epsilon.hat-2)/((a-1)*(s-g-(a-1)*epsilon.hat)))
[1] 0.5871795

まとめ

反復測定分散分析で登場する対象者内相関を考慮した自由度調整方法、Greenhouse-Geisserの方法とHuynh-Feldtの方法(Huynh-Feldt-Lecoutreの方法)をone by oneで確認してみて、ANOVA君の結果と照合した。

Greenhouse-Geisserの方法はConservativeすぎるとのことで、多重比較のBonferroniの方法を思い出す。ここはどちらかと言われれば、Huynh-Feldt-Lecoutreの方法を使うのがよいと言える。

おまけ:読み方をカタカナにすると・・・

今回登場する先生方の名前はなんて発音するのか想像がつかない、もしくは想像はなんとなくつくが自信がない。以下の4名の先生方の発音を調べてみた。カタカナ表記をすると、Geisserはガイサー、Huynhはフイン、Feldtはフェルト、Lecoutreはルクトル。結構意外!?


How to pronounce Geisser


🔴 How to Pronounce Huynh


How to Pronounce Feldt - PronounceNames.com

eow.alc.co.jp