統計講座

多重比較検定

B大学の学生の中から500人を無作為抽出を行い、スマートフォン利用時間について調べたデータからスマートフォン利用時間に関する学部差を知りたい。

下記データより、B大学の学生のスマートフォン利用時間について、学部間の違いについてデータから得られる知見を述べよ。

元データ
faculty	time
science	193
law	166
science	293
social	215
…

science	law	social	literature	agriculture	economy
205.32	169.44	18085	171.36	194.70	178.36

分散分析は複数グループの集合の中で、カテゴリにより（グループが「原因」となって）、平均が異なる、ということを示す(ex.血液型によって判定結果が変わる!)。マクロな視点。

その際(あるいはそれとは別に)複数カテゴリの集合の中で具体的にどの2グループの平均差が有意なのかを調べたくなる(A型とB型には差があるのだろうか?それ以外は?)。ミクロな視点。それには分散分析は役に立たず、さりとて通常の平均の差の検定は使えない(多重比較)。

2変数の比較を組み合せるのは「検定の繰り返し」となり、有意性が不当に出やすくなる。

多重検定

※分かりやすくするために有意水準α = 0.1とする。

集合A,B,Cの中から（対立）仮説「B ≠ C」を検定したい。

このとき「A ≠ B」「B ≠ C」「C ≠ A」、少なくともどれか一つを誤って有意判定する確率は1 - (1 - 0.1)^3→0.271である。つまり危険度0.271となり、有意水準0.1を上回ってしまう。

「いや、待て。おかしいやん。オレは「B ≠ C」だけしか言ってないがな。他の組合せの話なんてしてへんやん。なんでB,Cの話しかしてないのに他の組合せの確率関係あるねん」

「では君はA,B,C三つある中でなんでそのB,Cの組合せだけに着目したのかね？君はこっそりデータを見て、B ≠ Cに有意差が出ることを確認した上で、「B ≠ C」という仮説を立てたのであって、もしB ≠ Cには有意差が出ず、A ≠ Cに有意差が出たら今度はそれを仮説に立てて、君は有意性を主張するんじゃないのか？君は三つの組み合わせの中でどれか一つでも有意であればそれを持ち出して、有意性を主張することが出来る立場にいるのだから、逆にその三つの仮説の持つリスクをすべて背負う義務もあるのだよ」

「ひぃ！」

本当に仮説を立てる前にこっそりデータを確認したかどうかはこの際関係がない。当該組合せだけに着目する根拠をデータの外部で説明できれば2変数間の平均の差の検定で構わない。でも例えば「経済学部・法学部・理工学部・文学部・社会学部」があるなかで文学部と理工学部だけに着目する理由を説明するのは結構大変。というので、セットの中から二組を比較する場合、多重比較となるケースはかなりある。

というわけで多重比較なら結果を補正しなければならない。そのために様々な手法が開発されている。代表的な方法をかいつまんで紹介する。

p値を求めた後で修正する方法

検定を重ねれば、それだけ危険度が上がってしまうので、検定回数に応じてp値を修正する。

Bonferroniの調整

検定統計量に何をおくかは定めていない
「p値 * 検定回数」とする（=有意水準/検定回数）
有意水準0.1の検定を3回繰り返す → (1-0.1/3)^3 = 0.9032963 ほぼ0.9となる。
検定の回数が増えると有意性が出にくい(検出力が低い)

Holm法、Benjamini & Hochberg(BH)法

Bonferroni調整があまりに検出力が低いのを修正したもの。

p値を求める際の分布を変更する方法

p値を求める際、群数によって有意性が出にくくなるような分布を用いる。

Tukey-Kramer法

Studentのt検定を拡張して検定統計量を求める(分散分析の考え方を援用)
t分布の代わりにグループ数に影響される分布(スチューデント化した範囲(studentized range)の分布)を用いる
等分散性を前提

Games-Howell法

Welchのt検定を用いて検定統計量を求める
t分布の代わりにスチューデント化した範囲の分布を用いる
不等分散に対して頑健

なお上記検定はいずれも本質的には分散分析とは別の独立した検定である(検出力を上げるために、分散分析の結果を利用することもある)。分散分析の結果にかかわらず、この検定を行うことは可能である。ただ多重比較検定で有意であったとしても、検証群中のある二群の平均の差が有意であったということを論証するに過ぎず、その群を分割したカテゴリーによって（カテゴリーが「原因」で）平均が異なるという主張は分散分析から得られるものである。その意味で多重検定と分散分析はミクロとマクロで相互補完的な関係にあると言える。

分散分析: あるグループの集合のなかで、グループ間に平均差が生じるのはそのグループの分け方が原因である。
多重比較検定: あるグループの集合のなかで、特定の組み合わせにおいて平均差が異なっている。

t検定 + Bonferroni調整

通常のt検定を用いる。

t検定のp値(両側)に検定回数を掛けたものをp値とする。ただしp値が1を超える時には1とする。

p ← (1 - pt(t,自由度))*2*検定回数

data <- read.csv("https://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/smartphone00-B.csv", fileEncoding = "utf-8")
summary(data)
value <- data$time
group <- factor(data$faculty)

#グループごとの平均
mean.group <- tapply(value,group,mean)
#グループごとの分散
var.group <- tapply(value,group,var)
#グループごとのサンプルサイズ
N.group <- tapply(value,group,length)
#グループごとの自由度
df.group <- tapply(value,group,function(x){length(x)-1})
#グループの数
factor.N <- length(N.group)

#群ごとの比較（グループの全組合せに対して計算）
#平均差
mean.diff <- combn(factor.N,2,function(x){-diff(mean.group[x])})
#標準誤差
se <-combn(factor.N,2,function(x){sqrt(sum((sum(var.group[x]*df.group[x])/sum(df.group[x]))/N.group[x]))})
#自由度
df <- combn(factor.N,2,function(x)sum(df.group[x]))

t.value <- abs(mean.diff)/se #t値
p.value <- (1 - pt(t.value, df))*2*choose(length(N.group),2) #p値
p.value[which(p.value > 1)] <- 1

#結果のまとめ
Bonferroni.result <- cbind(mean.diff,t.value, p.value) 
rownames(Bonferroni.result) <- combn(levels(group),2,paste,collapse=" vs. ")

t検定拡張 + Bonferroni調整

通常のt検定ではなく、Tukey法で使われている検定統計量を用いる。全グループを対象として共通分散と自由度を計算する。

SPSSの「その後の多重比較」で「Bonferroni」を選択した時にはこの方法が用いられているようだ。

共通分散 ← sum(各集団の偏差平方和) / sum(各集団の自由度)
SE ← sqrt(共通分散 / 集団1のサンプルサイズ + 共通分散 / 集団2のサンプルサイズ)
検定統計量t ← abs(集団1の平均 - 集団2の平均)/SE
自由度 ← 全グループの自由度の総和

t検定のp値(両側)に検定回数を掛けたものをp値とする。ただしp値が1を超える時には1とする。

p ← (1 - pt(t,自由度))*2*検定回数

data <- read.csv("https://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/smartphone00-B.csv", fileEncoding = "utf-8")
summary(data)
value <- data$time
group <- factor(data$faculty)

#グループごとの平均
mean.group <- tapply(value,group,mean)
#グループごとのサンプルサイズ
N.group <- tapply(value,group,length)
#グループごとの自由度
df.group <- tapply(value,group,function(x){length(x)-1})
#グループの数
factor.N <- length(N.group)

#残差自由度
residual.df <- sum(df.group)
#残差変動
residual.ss <- sum(tapply(value,group,function(x){sum((x -mean(x))^2)}))
#残差分散(共通分散)
residual.var <-  residual.ss / residual.df

#群ごとの比較（グループの全組合せに対して計算）
#平均差
mean.diff <- combn(factor.N,2 ,function(x){-diff(mean.group[x])})
#標準誤差
se <- combn(factor.N, 2, function(x){sqrt(sum(residual.var/N.group[x]))})  

t.value <- abs(mean.diff)/se #t値
p.value <- (1 - pt(t.value, residual.df))*2*choose(length(N.group),2) #p値
p.value[which(p.value > 1)] <- 1

#結果のまとめ
Bonferroni.result <- cbind(mean.diff,t.value, p.value) 
rownames(Bonferroni.result) <- combn(levels(group),2,paste,collapse=" vs. ")

Tukey-Kramer

二グループの各組合せに対して、Studentのt検定とほぼ同じ計算を行う。ただし共通分散と自由度は全グループを対象として計算する。こうして分散分析の結果を援用することで検出力が上がる。

共通分散 ← sum(各集団の偏差平方和) / sum(各集団の自由度)
SE ← sqrt(共通分散 / 集団1のサンプルサイズ + 共通分散 / 集団2のサンプルサイズ)
検定統計量t ← abs(集団1の平均 - 集団2の平均)/SE

このtに2の平方根を掛けたもの(なんで？)はスチューデント化した範囲の分布（Q 分布）という確率分布に従う。

Q分布はグループの数と自由度(残差自由度）によって決まる(グループ数が増えれば同じ検定統計量のp値が大きくなる)。

自由度 ← 全グループの自由度の総和
p ← 1 - ptukey(t*sqrt(2), グループの数, 自由度)
グループ数が2の時、1 - ptukey(t*sqrt(2), 2, df) == (1 - pt(t, df))*2

data <- read.csv("https://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/smartphone00-B.csv", fileEncoding = "utf-8")
summary(data)
value <- data$time
group <- factor(data$faculty)

#グループごとの平均
mean.group <- tapply(value,group,mean)
#グループごとのサンプルサイズ
N.group <- tapply(value,group,length)
#グループごとの自由度
df.group <- tapply(value,group,function(x){length(x)-1})
#グループの数
factor.N <- length(N.group)

#残差自由度
residual.df <- sum(df.group)
#残差変動
residual.ss <- sum(tapply(value,group,function(x){sum((x -mean(x))^2)}))
#残差分散(共通分散)
residual.var <-  residual.ss / residual.df

#群ごとの比較（グループの全組合せに対して計算）
#平均差
mean.diff <- combn(factor.N,2 ,function(x){-diff(mean.group[x])})
#標準誤差
se <- combn(factor.N, 2, function(x){sqrt(sum(residual.var/N.group[x]))})  

t.value <- abs(mean.diff)/se #t値
p.value <- 1 - ptukey(t.value*sqrt(2), factor.N, residual.df) #p値

#結果のまとめ
Tukey.result <- cbind(mean.diff,t.value, p.value) 
rownames(Tukey.result) <- combn(levels(group),2,paste,collapse=" vs. ")

Games-Howell

Welchのt検定の検定統計量とTukey法の「スチューデント化した範囲の分布」を用いる。

SE ← sqrt(集団1の分散 / 集団1のサンプルサイズ + 集団2の分散 / 集団2のサンプルサイズ)
検定統計量t ← abs(集団1の平均 - 集団2の平均)/SE
自由度 ← 集団1，2によるWelchの等価自由度
p ← 1 - ptukey(t*sqrt(2), グループの数, 自由度)

data <- read.csv("https://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/smartphone00-B.csv", fileEncoding = "utf-8")
summary(data)
value <- data$time
group <- factor(data$faculty)

#グループごとの平均
mean.group <- tapply(value,group,mean)
#グループごとの分散
var.group <- tapply(value,group,var)
#グループごとのサンプルサイズ
N.group <- tapply(value,group,length)
#グループの数
factor.N <- length(N.group)

#群ごとの比較（グループの全組合せに対して計算）
#平均差
mean.diff <- combn(factor.N,2 ,function(x){-diff(mean.group[x])})
#標準誤差
se <- combn(factor.N, 2, function(x){sqrt(sum(var.group[x]/N.group[x]))})
#自由度
df <- combn(factor.N, 2, function(x){ sum(var.group[x]/N.group[x])^2/sum((var.group[x]/N.group[x])^2/(N.group[x]-1)) })

t.value <- abs(mean.diff)/se #t値
p.value <- 1 - ptukey(t.value*sqrt(2), factor.N, df) #p値

#結果のまとめ
GH.result <- cbind(mean.diff,t.value, df, p.value) 
rownames(GH.result) <- combn(levels(group),2,paste,collapse=" vs. ")

Rでの多重比較

Bonferroni法、Holm法、Benjamini & Hochberg法

p値を求めた後で修正する方法

p.values <- c(0.049,0.051,0.065)

#Bonferroni
p.adjust(p.values, "bonferroni")

#Holm
p.adjust(p.values, "holm")

#Benjamini & Hochberg
p.adjust(p.values, "BH")

Tukey-Kramer法

スチューデント化した範囲の分布を用いる方法

data <- read.csv("https://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/smartphone00-B.csv", fileEncoding = "utf-8")
summary(data)
x <- data$time
group <- factor(data$faculty)

#Tukey-Kramer
TukeyHSD(aov(x~group))