統計検定に関する事柄を忘れないようにまとめます.

R にてダネット検定 (Dunnett test) を行う.ダネット検定は,ひとつの対照群 (コントロール群) と2つ以上の処理群からなる多群のデータ中における対照群と処理群の各2群間の平均値の差についての検定を行う方法である.カナダの統計学者,Charles Dunnett によって開発された.R では,パッケージ multcomp の関数 glht() にて実行することができる.

最初に,Rを起動させ,以下のコマンドにてパッケージをインストールし,そのパッケージを読み込む.既にインストール済みの場合,1行目のコマンドは不要.

install.packages("multcomp", repos="http://cran.ism.ac.jp/")
library(multcomp)

以下のサンプルサイズが5からなる対照群と,それぞれ,4,6,6からなるデータX,YおよびZが得られたとき,CとX,CとY,CとZの各データ間における平均値の比較をダネット検定にて行う.帰無仮説 (H0) は対象の各2群間の平均値に差がないことである.

データC56, 48, 72, 60, 55
データX60, 62, 76, 84
データY78, 53, 62, 44, 90, 57
データZ77, 72, 83, 81, 91, 83

p値の計算

データを読み込む.チューキー・クレーマー検定で行っているように,データをベクトルにて,データラベル (水準名) を因子にて読み込んでも良いが,以下のように読み込むのも簡単である.最初に,データを以下のようにテキストファイル等で成型する.

fx	vx
C	56
C	48
C	72
C	60
C	55
X	60
X	62
X	76
X	84
Y	78
Y	53
Y	62
Y	44
Y	90
Y	57
Z	77
Z	72
Z	83
Z	81
Z	91
Z	83

次に,以上のデータをコピーし,以下のコマンドで変数 dx に読み込む.これは,clipboard の部分をファイル名に置き換えてファイルから読み込んでも良い.

$ dx=read.table("clipboard",header=T)

この変数 dx から fx 行を新たな変数 fx に,vx 行を新たな変数 vx に,以下のコマンドにて格納する.

$ fx=factor(dx$fx)
$ vx=dx$vx

以上で読み込んだデータをダネット検定する.ダネット検定は,従属変数 vx および独立変数 fx を用いて,summary(glht(aov(従属変数 ~ 独立変数), linfct=mcp(独立変数="Dunnett"))) のようにコマンドを打ち,実行する.複雑なコマンドであるが,分散分析を行う関数である aov(),結果の要約を表示する関数である summary(),実際のダネット検定の計算を行う関数である glht() に分解できる.実際には以下のように打つ.

$ summary(glht(aov(vx~fx),linfct=mcp(fx="Dunnett")))

結果は以下のようになる.ダネット検定では,両側検定に加え,右側検定と左側検定もできる.右側検定および左側検定はオプションにそれぞれ,alternative="g" または alternative="l" を加える.以上では何も指定していないのでデフォルトの両側検定が実行される.

 
         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = vx ~ fx)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)  
X - C == 0   12.300      7.902   1.557   0.3054  
Y - C == 0    5.800      7.132   0.813   0.7544  
Z - C == 0   22.967      7.132   3.220   0.0132 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

結果の X - C と表示されている行に処理群のデータXと対照群のデータCの比較の結果が示されている.YおよびZについても同様である.重要なのは Pr(>|t|) の列であり,ここに示されている値がp値である.もし,有意水準を5%と決めていた場合,ZとCの比較におけるp値は0.05以下であるので,これら2群の平均値には有意な差があるといえる.以上においては,対照群と処理群の比較が正確に行われているが,この glht はエラーなのか仕様なのか,水準名 (ここでの C,X,Y,Z) をアルファベット順 (水準名が数字の場合,降順) で並べたときの最も若いものを対照群にするので,そこに注意しなければ正確な対比較が行えない.

信頼区間の計算

また,以下の confint() コマンドを用いると信頼区間を算出することができる.オプション level にて表示させる信頼区間の幅を変更できる.以下のように level=0.99 を指定すると99%信頼区間が表示される.

$ confint(glht(aov(vx~fx),linfct=mcp(fx="Dunnett")),level=0.99)

これを実行した結果は以下のようになる.結果の lwr および upr はそれぞれ下方信頼限界および情報信頼限界,すなわち信頼区間の下限値 (lower) と上限値 (upper) を示している.この場合は,全ての対比較において,信頼区間に0が含まれているので有意差は出ていないことが確認できる.

 
         Simultaneous Confidence Intervals

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = vx ~ fx)

Quantile = 3.3492
99% family-wise confidence level
 

Linear Hypotheses:
           Estimate lwr      upr     
X - C == 0  12.3000 -14.1636  38.7636
Y - C == 0   5.8000 -18.0879  29.6879
Z - C == 0  22.9667  -0.9213  46.8546

Hatena Google+