実装関連事項

各種プログラミング言語の基本的な書き方やソフトウェア等の使用方法について.

R にてチューキー・クレーマー検定 (Tukey-Kramer test) を行う.チューキー・クレーマー検定は,チューキー検定 (Tukey test) とかチューキーの範囲検定 (Tukey range test) とかチューキーのHSD検定 (Tukey honestly significant difference test) 等とも呼ばれる多重比較の検定法のひとつであり,多群のデータ中における各2群間の平均値の差についての検定を行う方法である.アメリカの数学者,John Wilder Tukey によって考案された.日本語表記では,Tukey をテューキーとする場合が多い.R では,コマンド TukeyHSD にて実行することができる.

以下の,サンプルサイズがそれぞれ,12,10,12,11,14,12からなるデータAからFが得られたときの各データ間における平均値の比較をチューキー・クレーマー検定にて行う.帰無仮説 (H0) は対象の各2群間の平均値に差がないことである.

データA301, 311, 325, 291, 388, 402, 325, 361, 287, 261, 238, 362
データB197, 180, 178, 260, 247, 199, 179, 134, 163, 200
データC209, 331, 192, 155, 234, 290, 175, 116, 285, 216, 237, 301
データD343, 247, 316, 395, 324, 138, 245, 228, 214, 374, 235
データE225, 214, 269, 388, 375, 224, 177, 262, 319, 247, 125, 184, 135, 258
データF408, 314, 325, 347, 335, 285, 345, 257, 325, 344, 238, 352

p値の計算

最初に,以下のコマンドにて,上のデータAからFをまとめて変数 vx に格納する.

$ vx=c(301, 311, 325, 291, 388, 402, 325, 361, 287, 261, 238, 362, 197, 180, 178, 260, 247, 199, 179, 134, 163, 200, 209, 331, 192, 155, 234, 290, 175, 116, 285, 216, 237, 301, 343, 247, 316, 395, 324, 138, 245, 228, 214, 374, 235, 225, 214, 269, 388, 375, 224, 177, 262, 319, 247, 125, 184, 135, 258, 408, 314, 325, 347, 335, 285, 345, 257, 325, 344, 238, 352)

次に,以下のコマンドにて,データの各水準にラベルをつける.ラベルは新たな変数 fx に格納する.データAからFのサンプルサイズはそれぞれ,12,10,12,11,14,12であるが,それらが以下のコマンドの c(...) に対応する.データの読み込みは以上であるが,このような読み込み方法ではなく,分散分析を行うときのようにエクセルやテキストファイル形式のデータをデータフレームとしてそのまま読み込み,解析することもできる.

$ fx=factor(rep(c("A", "B", "C", "D", "E", "F"), c(12, 10, 12, 11, 14, 12)))

以上の従属変数 (vx) および独立変数 (fx) を用いて,チューキー・クレーマー検定は,TukeyHSD(aov(従属変数 ~ 独立変数)) のように実行する.このとき,等分散性を仮定できるかどうかを知っておくのは重要.コマンドは,実際には以下のように打つ.

$ TukeyHSD(aov(vx~fx))

これを実行した結果は以下のようになる.

 
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = vx ~ fx)

$fx
           diff         lwr        upr     p adj
B-A -127.300000 -205.208103 -49.391897 0.0001379
C-A  -92.583333 -166.865795 -18.300871 0.0064679
D-A  -42.909091 -118.861030  33.042848 0.5633612
E-A  -78.000000 -149.580371  -6.419629 0.0248489
F-A    1.916667  -72.365795  76.199129 0.9999996
C-B   34.716667  -43.191437 112.624770 0.7793179
D-B   84.390909    4.889428 163.892390 0.0311424
E-B   49.300000  -26.036166 124.636166 0.3986413
F-B  129.216667   51.308563 207.124770 0.0001060
D-C   49.674242  -26.277697 125.626182 0.3992971
E-C   14.583333  -56.997038  86.163705 0.9908014
F-C   94.500000   20.217538 168.782462 0.0051127
E-D  -35.090909 -108.402317  38.220498 0.7235594
F-D   44.825758  -31.126182 120.777697 0.5155242
F-E   79.916667    8.336295 151.497038 0.0199544

結果の diff は比較した群間の平均値の差である.例えば B-A の行の-127.3はデータBの平均がデータAの平均より-127.3大きいということを示している.また,lwr および upr はそれぞれ下方信頼限界および情報信頼限界,すなわち信頼区間の下限値 (lower) と上限値 (upper) を示している.この区間に0を含まない場合 (例えば,B-A は含まないが D-A の場合は含む),比較した2群間の差は少なくとも0ではないといえるので,2群間の平均値に有意差があるということになる.最後の p adj はp値である.もし有意水準を5%と決めていた場合,このp値が0.05より小さいとき,その2群間に有意差があるといえる.そのときは,上述の信頼区間内に0は含んでいないはずである.以上の検定結果より,有意水準5%ではA-B,A-C,A-E,B-D,B-F,C-F,E-F間にそれぞれ有意差があるといえる.

信頼区間の計算

また,表示させる信頼区間を変更するオプションに conf.level がある.0以上で1より小さい数を指定する.デフォルトでは95%信頼区間が表示されるが,これを99%信頼区間を表示させるようにするには以下のように打つ.

$ TukeyHSD(aov(vx~fx), conf.level=0.99)

これを実行した結果は以下のようになる.

 
  Tukey multiple comparisons of means
    99% family-wise confidence level

Fit: aov(formula = vx ~ fx)

$fx
           diff         lwr        upr     p adj
B-A -127.300000 -220.591776 -34.008224 0.0001379
C-A  -92.583333 -181.533552  -3.633115 0.0064679
D-A  -42.909091 -133.858440  48.040258 0.5633612
E-A  -78.000000 -163.714575   7.714575 0.0248489
F-A    1.916667  -87.033552  90.866885 0.9999996
C-B   34.716667  -58.575109 128.008443 0.7793179
D-B   84.390909  -10.808872 179.590690 0.0311424
E-B   49.300000  -40.911987 139.511987 0.3986413
F-B  129.216667   35.924891 222.508443 0.0001060
D-C   49.674242  -41.275107 140.623592 0.3992971
E-C   14.583333  -71.131242 100.297908 0.9908014
F-C   94.500000    5.549782 183.450218 0.0051127
E-D  -35.090909 -122.878329  52.696511 0.7235594
F-D   44.825758  -46.123592 135.775107 0.5155242
F-E   79.916667   -5.797908 165.631242 0.0199544

最初の検定で有意水準5%の場合は有意に差があった E-A において信頼区間に0が含まれており,実際に,p値も0.01 (1%) より大きな値を示している.

また,以下のように TukeyHSD の実行結果を plot() コマンドに供することで,信頼区間を図示することができる.

$ plot(TukeyHSD(aov(vx~fx)))

結果は以下のようになる.

ist_r_tukey_kramer_test_01.svg

オプション一覧

その他のオプションには以下のようなものがある.

オプション詳細
xaov(従属変数ベクトル ~ 独立変数ファクター) を入力する (必須).
orderedordered=T を指定すると,2群間で平均値が大きいほうから小さいほうを引いた場合の diff や lwr および upr を表示.
conf.level表示させる信頼区間の指定.
Hatena Google+