MENU
| 統計ER
  1. ホーム
  2. R
  3. R で Greenhouse-Geisser 補正と Huynh-Feldt 補正の違いを説明する

R で Greenhouse-Geisser 補正と Huynh-Feldt 補正の違いを説明する

R 分散分析 反復測定分散分析

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

＞＞もう統計で悩むのは終わりにしませんか？ 

↑1万人以上の医療従事者が購読中

ANOVA君で実行

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

http://riseki.php.xdomain.jp/index.php?ANOVA%E5%90%9Briseki.php.xdomain.jp

source("C:\\Users\\happy\\Downloads\\anovakun_482.txt")
anovakun(dat, "sA", 4, nopost=TRUE)

単純な反復測定分散分析は sA と表現する。

ちなみに、分割プロットデザイン は AsBと書く。

詳しくは以下を参照。

あわせて読みたい
R で分割プロット分散分析を行う方法 分割プロット分散分析 split-plot design は、群分けのある反復測定分散分析のこと 分割プロット分散分析を R で行う方法。 分割プロット分散分析のサンプルデータ 今回…

結果は以下の通り。

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 --------------------///

自由度補正をしない計算

全平方和（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の結果はこちら。

<< 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 の結果

やり方は、分散共分散行列を作ってイプシロン $\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

Huynh-Feldt 補正 Step by step でやってみる

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君の結果と一致している。

> (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 Geisser

🔴 How to Pronounce Huynh

How to Pronounce Feldt – PronounceNames.com

https://eow.alc.co.jp/search?q=Le+Coutre

参考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

R 分散分析 反復測定分散分析
よかったらシェアしてね！
  • URLをコピーしました！
  • URLをコピーしました！

リサーチクエスチョン探し？データ分析？論文投稿？、、、で、もう悩まない！

第1章：臨床研究ではなぜ統計が必要なのか？計画することの重要性
  • 推定ってどんなことをしているの？
  • 臨床研究を計画するってどういうこと？
  • どうにかして標本平均を母平均に近づけられないか？
第2章：研究目的をどれだけ明確にできるのかが重要
  • データさえあれば解析でどうにかなる、という考え方は間違い
  • 何を明らかにしたいのか？ という研究目的が重要
  • 研究目的は4種類に分けられる
  • 統計専門家に相談する上でも研究目的とPICOを明確化しておく
第3章：p値で結果が左右される時代は終わりました
  • アメリカ統計協会（ASA）のp値に関する声明で指摘されていること
  • そうは言っても、本当に有意差がなくてもいいの…？
  • なぜ統計専門家はp値を重要視していないのか
  • 有意差がない時に「有意な傾向があった」といってもいい？
  • 統計を放置してしまうと非常にまずい
第4章：多くの人が統計を苦手にする理由
  • 残念ながら、セミナー受講だけで統計は使えません。
  • インプットだけで統計が使えない理由
  • どうやったら統計の判断力が鍛えられるか？
  • 統計は手段なので正解がないため、最適解を判断する力が必要
第5章：統計を使えるようになるために今日から何をすれば良いか？
  • 論文を読んで統計が使えるようになるための5ステップ
第6章：統計を学ぶために重要な環境
  • 統計の3つの力をバランスよく構築する環境

『統計を身につける5つのステップ』を
通常1,650円 → 無料でプレゼント

以下のボタンをクリックして、画面に出てくる指示に従って、必要事項を記入してください。

通常価格：1,650円　→　無料でお届け

【無料】統計を使いこなすための秘策を学ぶ

この記事を書いた人

toukei-erのアバター toukei-er

統計 ER ブログ執筆者

元疫学研究者

統計解析が趣味

関連記事

コメント

コメント一覧 （2件）

コメントする

目次