Sponsored Links
データの読み込みと準備
サンプルデータは、統計教育推進委員会が提供してくださっているいろいろな樹木の葉の長さと幅のデータ[1]を使います。このデータにはイチョウ、ヤブツバキ、クスノキ、ハナミヅキの葉の長さと幅のデータが収められています。今回は、イチョウ、ヤブツバキ、クスノキ、ハナミヅキそれぞれの葉の長さを使います。
> x <- read.table(pipe("pbpaste"), header=TRUE) # クリップボードからデータを読み込み、xに代入。ちなみにこちらはMac用のコードでWindowsは、read.table("clipboard", header=TRUE)。
> x # 読み込んだデータを表示。
icho tsubaki kusunoki hanamizuki
1 4.0 6.8 7.0 8.0
2 3.2 8.2 8.2 7.0
3 2.5 6.5 8.7 7.5
4 2.7 6.7 8.0 7.5
5 3.2 7.0 6.9 7.0
〜〜〜〜〜〜〜〜〜〜（以下略）〜〜〜〜〜〜〜〜〜〜分散分析を行うために、データの表現方式をアンスタック形式からスタック形式に変えます。上のがアンスタック形式です。
> y <- stack(x) # アンスタックからスタックに変換し、yに代入。
> y # yを表示。
values ind
1 4.0 icho
2 3.2 icho
〜〜〜〜〜〜〜〜〜〜（以下略）〜〜〜〜〜〜〜〜〜〜
51 6.8 tsubaki
52 8.2 tsubaki
〜〜〜〜〜〜〜〜〜〜（以下略）〜〜〜〜〜〜〜〜〜〜
101 7.0 kusunoki
102 8.2 kusunoki
〜〜〜〜〜〜〜〜〜〜（以下略）〜〜〜〜〜〜〜〜〜〜
151 8.0 hanamizuki
152 7.0 hanamizuki
〜〜〜〜〜〜〜〜〜〜（以下略）〜〜〜〜〜〜〜〜〜〜
> str(y) # yの構造を確認。
'data.frame': 200 obs. of 2 variables:
$ values: num 4 3.2 2.5 2.7 3.2 4 3.6 3.3 2.9 2.8 ...
$ ind : Factor w/ 4 levels "icho","tsubaki",..: 1 1 1 1 1 1 1 1 1 1 ...
> colnames(y) <- c("leaf_l", "breed") # 変数名がピンとこないので、葉の長さをleaf_l、品種をbreedに変更。
> str(y) # yの変数名が変わっていることを確認。
'data.frame': 200 obs. of 2 variables:
$ leaf_l: num 4 3.2 2.5 2.7 3.2 4 3.6 3.3 2.9 2.8 ...
$ breed : Factor w/ 4 levels "icho","tsubaki",..: 1 1 1 1 1 1 1 1 1 1 ...
等分散性の検定
一元配置分散分析はそれぞれの水準が等分散を仮定して計算します。等分散でない場合はノンパラメトリック検定としてKruskal-Wallis検定を行います。等分散性の検定にはF検定の一種であるLevene検定を使うことにします。こちらはcarパッケージを使うと便利です。Levene検定は二つ以上の母集団が等分散であることを帰無仮説としますので、p>0.05なら帰無仮説は棄却されず、母集団は等分散であると考えます。p<0.05なら等分散でないことになります。
> library(car)
> leveneTest(y$leaf_l, y$breed, center = mean)
Levene's Test for Homogeneity of Variance (center = mean)
Df F value Pr(>F)
group 3 11.996 3.035e-07 ***
196 p<0.05未満で等分散ではありませんでしたので、等分散ではなく、ノンパラメトリック検定としてKruskal-Wallis検定を行います。
クラスカル=ウォリス検定の実行
kruskal.test(変数 ~ 因子)により計算できる。
> kruskal.test(y$leaf_l ~ y$breed)
Kruskal-Wallis rank sum test
data: y$leaf_l by y$breed
Kruskal-Wallis chi-squared = 107.93, df = 3, p-value < 2.2e-16p値がとても小さいので、有意な群間差が認められました。
多重比較検定の実行
クラスカル=ウォリス検定の後の多重比較検定には、Mann–Whitney U検定で得られるp値をBonferroniやHolmによる補正をして用いる方法があります。計算方法は以下の通りです。
# Bonferroniの方法
pairwise.wilcox.test(y$leaf_l, y$breed, p.adj="bonferroni", exact=F)
# Holmの方法
pairwise.wilcox.test(y$leaf_l, y$breed, p.adj="holm", exact=F)
Post hocとしてノンパラメトリックの2群の比較検定を繰り返す方法以外にも、Steel-Dwass法やGames-Howelという方法があるのですね[1]。
ここではSteel-Dwass法を紹介します。NSM3というパッケージを使ってみます。
> library(NSM3) # パッケージNSM3を読み込む。
> pSDCFlig(y$leaf_l, y$breed) # pSDCFligという関数を使う。
Ties are present, so p-values are based on conditional null distribution.
Group sizes: 50 50 50 50
Using the Monte Carlo (with 10000 Iterations) method:
For treatments icho - tsubaki, the Dwass, Steel, Critchlow-Fligner W Statistic is 12.0758.
The smallest experimentwise error rate leading to rejection is 0 .
For treatments icho - kusunoki, the Dwass, Steel, Critchlow-Fligner W Statistic is 10.9889.
The smallest experimentwise error rate leading to rejection is 0 .
For treatments icho - hanamizuki, the Dwass, Steel, Critchlow-Fligner W Statistic is 11.5301.
The smallest experimentwise error rate leading to rejection is 0 .
For treatments tsubaki - kusunoki, the Dwass, Steel, Critchlow-Fligner W Statistic is -1.5078.
The smallest experimentwise error rate leading to rejection is 0.7176 .
For treatments tsubaki - hanamizuki, the Dwass, Steel, Critchlow-Fligner W Statistic is -5.5426.
The smallest experimentwise error rate leading to rejection is 6e-04 .
For treatments kusunoki - hanamizuki, the Dwass, Steel, Critchlow-Fligner W Statistic is -2.5306.
The smallest experimentwise error rate leading to rejection is 0.2811 .
警告メッセージ:
1: factorial(sum(outp$n)) で: 'gammafn' 中の値が範囲を超えています
2: factorial(sum(outp$n)) で: 'gammafn' 中の値が範囲を超えています Steel-Dwass法は、統計学的手法の話題 – 生物科学研究所のこちらのページを参考にさせていただきました。pSDCFligではmethod=”Exact”で正確検定を指定できるようですが、今回使ったサンプルサイズが大きかったからか上手く動かなかったので、何も設定しませんでしたが、何も設定しないとmethod=”Monte Carlo”として計算されるようです。それでも計算にかなりの負荷がかかるらしく、数分固まっていました。
Sponsored Links
