• 検索結果がありません。

多重比較補正付きの検定関数

ドキュメント内 ためになった他の人のサイト script of (ページ 78-83)

一つのデータ

X

から複数の検定を行うとする.各検定の帰無仮説を

H

1

, H

2

, . . . , H

k

とし,集合

C

i

(α)

を帰無仮説

H

i の下での有意水準

α

の棄却域とする.全ての帰無仮説

H

i が同時に成り立つという帰無仮説

H

を考える.ある帰無仮説

H

i が棄却されれば

H

も棄却される

(

おなじことだが,全ての

H

i が棄却されなければ

H

も棄却されない

)

よう にするのが一見自然と思われる

(

これは

i

C

i

(α)

H

の棄却域にすることと同値

)

.し かし,こうした検定の有意水準は普通

α

にはならず,一般的にいえることは

P { X ∈ ∪

i

C

i

(α) } ≤ X

i

P { X ∈ C

i

(α) } = kα

だけである

(

検定の多重比較問題

)

.ここで

P { X ∈ ∪

i

C

i

(α/k) } ≤ k(α/k) = α

であるから

i

C

i

(α/k)

は帰無仮説

H

の有意水準

α

の棄却域を与える.つまり,各検定

H

i を有意水準

α/k

で行えば,その結果は同時帰無仮説

H

の有意水準

α

の検定として使 える.これは

Bonferroni

の検定の多重比較補正と呼ばれる方法である.これは検定の多 重比較補正の基本となる考え方であるが,当然

k

が少し大きくなると

α/k

が小さくなり すぎ,実際上無意味となる.

検定の多重比較補正法としては,より精緻な方法が幾つも提案されているが,いずれも 個別の検定をある有意水準

α

i で行えば,それが結果として同時帰無仮説

H

の有意水準

α

の検定となるような補正有意水準

α

1

, α

2

, . . . , α

k を与える

(

一般に

α

i

≤ α )

.同じこ とを

p

値の言葉で表現すると,本来の

p

p

1

, p

2

, . . . , p

k の多重比較補正値と呼ばれる 数値

p

1

, p

2

, . . . , p

k

(

一般に

p

i

≥ p

i

)

を定め,ある

i

p

i

≤ α

であれば,同時帰無仮説

H

が有意水準

α

で棄却されたことになるようにする.例えば

Bonferroni

の多重比較補 正法に対しては,

p

値の多重比較補正を

p

i

= kp

iと取れば良い

:

p

i

≤ α ⇐⇒ p

i

≤ α/k ⇐⇒ H

i を有意水準

α/k

で棄却.

3.3.1 p 値のベクトルの多重比較補正 p.adjust()

p.adjust()

p

値のベクトルを与えて,幾つかの方法を用いて多重比較補正された

p

値のベクトルを返す.

書式:

p.adjust(p, method=p.adjust.methods, n=length(p)) p.adjust.methods #

可能な手法名の文字ベクトル 引数:

p p

値のベクトル

method

補正方法.以下の文字列のいずれか

"holm","hochberg","hommel","bonferroni","BH","BY","fdr","none"

n

比較の個数

返り値: 補正された

p

値のベクトル

補正方法の一つは

Bonferroni

補正

("bonferroni")

で,

p

値は比較回数倍される.

それほど保守的でない補正法がそれぞれ

Holm ("holm")

Hochberg ("hochberg")

Hommel ("hommel")

そして

Benjamini & Hochberg ("fdr")

で与えられる.無補正オ プション

("none")

もある.補正法をオプションとして持ち,それを引き数

p.adjust

に指定する必要がある手法の便宜用に,手法名の文字列ベクトル

p.adjust.methods

が ある.

最初の

4

つの方法はファミリ単位の誤差率を強力にコントロールするように作られてい

る.

Bonferroni

補正はそのままでは,如何なる場合も使用可能な

Holm

法より劣るので

使用する理由は無い.

Hochberg

Hommel

の方法は,検定が独立か,非負の連関を持 つ

(Sarkar, Sarkar & Chang)

という仮定下で有効である.

Hommel

法は

Hochberg

法よ りもより強力であるが,差異は普通小さく,

Hochberg

法は

p

値をより高速に計算する.

Benjamini & Hochberg

"fdr"

法は棄却された仮説中の誤った発見の平均の割合であ る「誤発見率

(false discovery rate)

」をコントロールする.誤発見率はファミリ単位の誤 判別率よりも厳しい条件であり,従って

Benjamini & Hochberg

法は他の手法よりもよ り強力な方法である.

関連:

stats

パッケージ中の

pairwise.t.test()

等の

pairwise.*()

関数.

> p.adjust.methods # 可能な手法名の文字ベクトル

[1] "holm" "hochberg" "hommel" "bonferroni" "fdr" "none"

> x <- rnorm(10, m=c(rep(0,5),rep(3,5))) # 平均0,3の正規乱数5個ずつ

> p <- 2*pnorm(-abs(x)) # テスト用のp値ベクトル

> round(p, 3)

[1] 0.165 0.282 0.621 0.435 0.696 0.040 0.001 0.102 0.001 0.010

> round(p.adjust(p), 3) # 既定のHolm法による補正

[1] 0.825 1.000 1.000 1.000 1.000 0.280 0.006 0.615 0.006 0.084

> round(p.adjust(p,"bonferroni"), 3) # Bonferoni法による補正 [1] 1.000 1.000 1.000 1.000 1.000 0.399 0.006 1.000 0.007 0.105

> round(p.adjust(p,"fdr"), 3) # false discovery rate法による補正 [1] 0.275 0.402 0.690 0.544 0.696 0.100 0.003 0.205 0.003 0.035

3.3.2 多重比較補正を伴う対毎の比率の比較 pairwise.prop.test()

pairwise.prop.test()

は,多重比較補正を伴う対毎の比率の比較を行う.

書式:

pairwise.prop.test(x, n, p.adjust.method=p.adjust.methods, ...)

引数:

x

成功回数のベクトル.または,それぞれ成功・失敗回数を

2

つの列に持つ行列

n

試行回数のベクトル.もし

x

が行列なら無視される

p.adjust.method p

値を補正する方法

(p.adjust()

を参照

) ... prop.test()

に引き渡される追加引き数

返り値: クラス

"pairwise.htest"

を持つオブジェクト 関連:

prop.test(), p.adjust()

# 患者4グループ中の喫煙者数のデータ.各グループ対毎に喫煙者の割合が同一という帰無仮説の検定

# prpo.test()の参考例を参照(4グループ込みの比率の同一性の検定)

> smokers <- c(83, 90, 129, 70)

> patients <- c(86, 93, 136, 82)

# 各グループ対毎の喫煙者数比率の同一性を検定し,p値をHolm法で多重比較補正

> pairwise.prop.test(smokers, patients)

$method

[1] "Pairwise comparison of proportions"

$data.name

[1] "smokers out of patients"

# 対毎の比較の補正済みp値.第2,4グループ間の補正p値が0.1未満なので帰無仮説は10%有意

$p.value

1 2 3

2 1.0000000 NA NA

3 1.0000000 1.00000000 NA 4 0.1185648 0.09321728 0.1237680

$p.adjust.method # Holm法で多重比較補正(既定)

[1] "holm"

attr(,"class") [1] "pairwise.htest"

(以下省略)

> pairwise.prop.test(smokers, patients, p.adjust.method="fdr") # "fdr"法で多重比較補正 Pairwise comparisons using Pairwise comparison of proportions

data: smokers out of patients

1 2 3

2 1.000 - 3 0.965 0.965 -4 0.062 0.062 0.062

P value adjustment method: fdr (以下省略)

3.3.3 多重比較補正を伴うグループ水準間の対毎の t 検定による比較

pairwise.t.test()

pairwise.t.test()

は,多重比較補正を伴うグループ水準間の,対毎の

t

検定による

比較を行う.

書式:

pairwise.t.test(x, g, p.adjust.method = p.adjust.methods, pool.sd = TRUE, ...)

引数:

x

応答ベクトル

g

グルーピング用のベクトルか因子

p.adjust.method p

値を補正する方法

(p.adjust()

を参照

)

pool.sd

プールした標準偏差の使用を許すか許さないか

... t.test()

に引き渡される追加引き数

返り値: クラス

"pairwise.htest"

を持つオブジェクト

関連:

t.test(), p.adjust()

# NY市の毎日の大気観測値データairquality.各月対毎にオゾン量平均が等しいという帰無仮説のt検定

> attach(airquality) # 成分を名前で参照できるようにする

> Month <- factor(Month, labels = month.abb[5:9]) # 月数を省略文字名で因子化

> pairwise.t.test(Ozone, Month)

$method # プールした標準偏差を用いたt検定

[1] "t tests with pooled SD"

$data.name

[1] "Ozone and Month"

# 多重比較補正済みp値.(Jul,Sep),(Aug,Sep)間の補正p値が0.005より小さく,0.5%有意

$p.value

May Jun Jul Aug

Jun 1.0000000000 NA NA NA

Jul 0.0002638036 0.05112741 NA NA

Aug 0.0001949061 0.04987333 1.000000000 NA Sep 1.0000000000 1.00000000 0.004878798 0.003878108

$p.adjust.method # Holm法で多重比較補正法(既定)

[1] "holm"

attr(,"class") [1] "pairwise.htest"

# Bonferroni法で多重比較補正すると(Aug,Sep)間のみでp値が0.005以下

> pairwise.t.test(Ozone, Month, p.adj = "bonf")

$method # プールした標準偏差を用いたt検定

[1] "t tests with pooled SD"

$data.name

[1] "Ozone and Month"

$p.value

May Jun Jul Aug

Jun 1.0000000000 NA NA NA

Jul 0.0002931151 0.10225483 NA NA

Aug 0.0001949061 0.08312222 1.000000000 NA Sep 1.0000000000 1.00000000 0.006969712 0.004847635

$p.adjust.method # Bonferroni法で多重比較補正

[1] "bonferroni"

attr(,"class") [1] "pairwise.htest"

> pairwise.t.test(Ozone, Month, pool.sd = FALSE)

$method # プールした標準偏差を用いないt検定

[1] "t tests with non-pooled SD"

$data.name

[1] "Ozone and Month"

$p.value # 今度は1%有意という結果

May Jun Jul Aug

Jun 1.0000000000 NA NA NA

Jul 0.0002646679 0.01527179 NA NA

Aug 0.0019517090 0.02134659 1.000000000 NA Sep 0.8632062149 1.00000000 0.005888826 0.01720974

$p.adjust.method # Holm法で多重比較補正

[1] "holm"

attr(,"class") [1] "pairwise.htest"

> detach() # 成分を消去

3.3.4 Wilcoxon のランク和検定に対する多重比較補正

pairwise.wilcox.test()

pairwise.wilcox.test()

はグループ水準の各組合せ毎に

Wilcoxon

のランク和検定 を実施し,

p

値に対する多重比較補正を行う.

書式:

pairwise.wilcox.test(x,g,p.adjust.method=p.adjust.methods,...)

引数:

x

応答ベクトル

g

グルーピング用ベクトル,もしくは因子

p.adjust.method p

値を補正する方法

(p.adjust()

を参照

) ... wilcox.test()

に引き渡される追加引き数

返り値: クラス属性

"pairwise.htest"

を持つオブジェクト 関連:

wilcox.test()

p.adjust()

# NY市の大気観測値データairquality使用.月の組合せ毎にwilcox.test()を適用し多重比較補正

> Month <- factor(Month, labels = month.abb[5:9]) # 5,6,7,8,9月名を因子に変換

> pairwise.wilcox.test(Ozone, Month)

$method

[1] "Wilcoxon rank sum test"

$data.name

[1] "Ozone and Month"

# 多重比較補正されたp値.(May,Jul)間で0.0003以下なので0.05%有意

$p.value

May Jun Jul Aug

Jun 0.5775049432 NA NA NA

Jul 0.0002996390 0.08481748 NA NA

Aug 0.0010872705 0.12953880 1.000000000 NA Sep 0.4743723336 1.00000000 0.005954083 0.02273569

$p.adjust.method # Holm法で補正(既定)

[1] "holm"

attr(,"class") [1] "pairwise.htest"

Warning messages: # データにタイがあるので警告がでる(省略)

# Bonferroni法で多重比較.再び0.05%有意という結果

> pairwise.wilcox.test(Ozone, Month, p.adj = "bonf")

$method

[1] "Wilcoxon rank sum test"

$data.name

[1] "Ozone and Month"

$p.value

May Jun Jul Aug

Jun 1.0000000000 NA NA NA

Jul 0.0002996390 0.1413625 NA NA

Aug 0.0012080783 0.2590776 1.000000000 NA Sep 1.0000000000 1.0000000 0.007442604 0.03247955

$p.adjust.method [1] "bonferroni"

attr(,"class") [1] "pairwise.htest"

Warning messages: # データにタイがあるので警告がでる(省略)

> detach()

ドキュメント内 ためになった他の人のサイト script of (ページ 78-83)