Greenhouse-Geisser 補正と Huynh-Feldt 補正の違いの解説。
- Greenhouse-Geisser 補正と Huynh-Feldt 補正とは何か？
- Greenhouse-Geisser 補正と Huynh-Feldt 補正 サンプルデータの準備
- Greenhouse-Geisser 補正と Huynh-Feldt 補正 ANOVA君で実行
- Greenhouse-Geisser 補正と Huynh-Feldt 補正 自由度補正をしない計算
- Greenhouse-Geisser 補正と Huynh-Feldt 補正 自由度調整をした場合
- まとめ
- おまけ：各研究者の名前の読み方
- 参考PDF
Greenhouse-Geisser 補正と Huynh-Feldt 補正とは何か？
繰り返し測定分散分析に登場する Greenhouse-Geisser 補正とかHuynh-Feldt 補正とかは何をしているのか？
一言で言えば、同じ対象者を複数回測定するために起きる相関を考慮した検定方法で、自由度を調整する方法だ。
Greenhouse-Geisser 補正と Huynh-Feldt 補正 サンプルデータの準備
まずデータを準備する。
要因は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
Greenhouse-Geisser 補正と Huynh-Feldt 補正 ANOVA君で実行
こちらから最新のANOVA君をダウンロードしてお借りする。
source("C:\\Users\\happy\\Downloads\\anovakun_482.txt") anovakun(dat, "sA", 4, nopost=TRUE)
単純な反復測定分散分析は sA と表現する。
ちなみに、分割プロットデザイン は AsBと書く。
詳しくは以下を参照。
結果は以下の通り。
DESCRIPTIVE STATISTICSは、上記の値と一致した計算結果が表示されている。
今回 step by step でトライするのは 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 --------------------///
Greenhouse-Geisser 補正と Huynh-Feldt 補正 自由度補正をしない計算
全平方和（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 補正とHuynh-Feldt 補正だ。
Greenhouse-Geisser 補正 ANOVA君の結果
gg=TRUEと付けると、Greenhouse-Geisserの方法でANOVA TABLEを出力してくれる。
anovakun(dat, "sA", 4, nopost=TRUE, gg=TRUE)
ANOVA TABLEの結果はこちら。
<< 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
p valueが0.0142から0.0600に上昇しているのがわかる。
dfも少数点以下があって、小さい値だ。
Greenhouse-Geisserの方法はとてもConservative（控えめ）な方法で、結果は統計学的有意ではなくなった。
Greenhouse-Geisser 補正 Step by step の結果
やり方は、分散共分散行列を作って を推定する。
母集団分散共分散行列（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 補正のの推定値。
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
ちなみに、 の推定値 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
Huynh-Feldt 補正 Step by step でやってみる
Huynh-Feldt 補正は、Greenhouse-Geisser 補正のの推定値（
）を使ってより適切に修正した方法なのだ。
Huynh-Feldtのは、
よりも大きくなり、統計学的有意に出やすい。
事実、今回は 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君の結果と一致している。
> (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
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 と出力されている。
Huynh-Feldt-Lecoutre 補正（要因数が 2 以上の場合分割プロットデザインなど）
Huynh-Feldt 補正は、要因が2つ以上の場合、例えば 分割プロットデザイン（ 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 補正）をstep by step で確認してみて、ANOVA君の結果と照合した。
Huynh-Feldt-Lecoutre 補正を使うのがよい。
おまけ：各研究者の名前の読み方
今回登場する先生方の名前はなんて発音するのか想像がつかない、もしくは想像はなんとなくつくが自信がない。
以下の4名の先生方の発音を調べてみた。カタカナ表記をすると、Geisserはガイサー、Huynhはフイン、Feldtはフェルト、Lecoutreはルクトル。
How to Pronounce Feldt - PronounceNames.com
参考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
