R言語でノンパラメトリック検定であるウィルコクソンの順位和検定とマン・ホイットニーのU検定を実行する関数やその実行例を紹介します。
実行例では、正規性をもたないデータセットに対し、ウィルコクソンの順位和検定やマン・ホイットニーのU検定で位置母数に関する検定を行う例について見ていきます。
この記事で扱うRスクリプトは以下からダウンロードできます。
また、ウィルコクソンの順位和検定について詳しく見たい方は次の記事を参照してください。
-
【統計学】ウィルコクソンの順位和検定
ノンパラメトリックの検定の一種であるウィルコクソンの順位和検定について解説する。 ウィルコクソンの順位和検定の検定統計量と棄却域を紹介し、それらを導出する方法を紹介する。 順位和統計量の期待値と分散を ...
続きを見る
ポチップ
ポチップウィルコクソンの順位和検定
ノンパラメトリック検定であるウィルコクソンの順位和検定の概要と検定を実行する関数を紹介します。
検定の概要
ウィルコクソンの順位和検定の概要は以下の通りです。
次に示すように、分布（モーメント）を仮定していないのが特徴です。
パラメトリック検定で代表的なt検定では母集団分布に正規性を仮定しましたが、ウィルコクソンの順位和検定ではしないため、標本数が十分であればどのようなデータに対しても位置母数に関する検定を行うことができます。
ウィルコクソンの順位和検定
\(x_1, \ldots, x_{n_1}\)と\(y_1, \ldots, y_{n_2}\)をそれぞれ確率密度関数\(f(x)\)と\(g(y)\)分布関数\(F(x)\)と\(G(y)\)をもつ母集団からの標本数\(n_1\)、\(n_2\)の標本とする。また、これらの標本\(x_1, \ldots, x_{n_1}\)と\(y_1, \ldots, y_{n_2}\)を小さい順に並べ、ぞれぞれの母集団に対応する順位を\(r_1, \ldots, r_{n_1}\)、\(s_1, \ldots, s_{n_2}\)とする。今、母集団分布の尺度母数は等しいと仮定する。すなわち\(F(x) = G(x + \Delta)\)。ここで、「2つの母集団分布の位置母数が等しい」という次の仮説を考える。
この仮説の検定統計量として次を用いる。
ここに、\(W\)は次で定義される順位和統計量である。
また、\(W\)の期待値と分散は次で与えられる。
有意水準\(\alpha\)の検定の棄却域は次のとおりである。
関数wilcox.test
続いてウィルコクソンの順位和検定の関数を紹介します。
次のwilcox.testによって、ウィルコクソンの順位和検定を実行することができます。
2標本検定の場合はマン・ホイットニーのU検定と同じ結果が得られます。
wilcox.test(x, …)
# S3 method for default
wilcox.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, exact = NULL, correct = TRUE, conf.int = FALSE, conf.level = 0.95, …)
# S3 method for formula wilcox.test(formula, data, subset, na.action, …)
また、wilcox.testの引数は以下の通りです。
| x | numeric型の観測値のベクトル。有限で無い値は除外される。 |
| y | numeric型の観測値のベクトル。二標本検定の際に用いる。有限で無い値は除外される。 |
| alternative | character型。"two.sided"で両側検定、"greater"で右片側検定、"less"で左片側検定を指定できる。 |
| mu | numeric型。帰無仮説の下での位置パラメータ。 |
| paired | logical型。TRUEのとき対応のある検定をする。 |
| exact | logical型。TRUEのとき正確なp値を算出する。 |
| correct | logical型。TRUEのときp値の計算の際に正規近似に連続性補正を適用する。 |
| conf.int | logical型。TRUEのとき信頼区間を算出する。 |
| conf.level | numeric型。信頼区間の信頼水準。 |
| formula | lhs~rhsの形から成るformula型。lhsは観測値を表すnumeric型の変数名であり、rhsはグループをlevelsにもつfactor型の変数名。 |
| data | formulaで記述した変数をもつデータフレームまたは行列。 |
| subset | どの行を検定に用いるか指定するためのベクトル。 |
| na.action | 欠測値NAが含まれているときどうするか。 |
実行例
wilcox.testの実行例を紹介します。
1標本検定と2標本検定の2つに分けて実行例をみていきます。
ヒストグラムのプロットなどで次のパッケージggplot2を使うので事前にインストールしといてください。
1 2 | install.packages("ggplot2") library(ggplot2) |
1標本検定
1標本のウィルコクソンの順位和検定を実行例をみていきます。
データセットとして次のtreesを用います。
1 2 | #1標本検定 dataset <- trees |
データセットtreesは木のGirth（周囲の長さ）、Height（高さ）、Volume（重さ）の3つの変数から構成されます。
1 2 3 4 5 6 7 8 | > head(dataset) Girth Height Volume 1 8.3 70 10.3 2 8.6 65 10.3 3 8.8 63 10.2 4 10.5 72 16.4 5 10.7 81 18.8 6 10.8 83 19.7 |
今、Volumeの位置母数（中央値）に関心があるとします。
Volumeがどのようなデータであるのか把握するために、まずヒストグラムを描いてみます。つぎを実行するとVolumeのヒストグラムがプロットされます。
1 2 3 | ggplot(dataset, aes(x = Volume, y = ..density..)) + #ヒストグラム geom_histogram(position = "identity", color = "gray", bins = 10, alpha = 0.8) + geom_density(stat = "density", position = "identity", color = "black", fill = "black", alpha = 0.4) |
上の図のヒストグラムからVolumeは左右対称的に分布しておらず、正規性を持たないことが確認できます。
このようなデータに対しt検定などの母集団分布に正規性を仮定する検定を用いるのは不適切であるため、ウィルコクソンの順位和検定により位置母数の検定を行いたいと思います。
「Volumeの中央値は20であるか」という仮説検定は、wilcox.testを用いることで次のように行えます。
1 2 | mu0 <- 20 testResult <- wilcox.test(dataset$Volume, mu = mu0) |
testResultを参照すると次のように検定結果がコンソールに表示されます。
1 2 3 4 5 6 7 | > testResult Wilcoxon signed rank test with continuity correction data: dataset$Volume V = 399, p-value = 0.003184 alternative hypothesis: true location is not equal to 20 |
ポイント
「p-value = 0.003184」と書いてあることから、Volumeの中央値は20であるとは言えないという結論が得られました（Volumeの中央値は20ではないことが分かりました）。
p値や検定統計量などの検定結果を取得したいときは、次のようにtestResultの後ろに$を付ければ大丈夫です。
1 2 3 | stat <- testResult$statistic #検定統計量 pValue <- testResult$p.value #p値 nullValue <- testResult$null.value #帰無仮説の下での中央値 |
補足として、これらの検定結果をcsvファイルなどに保存する場合は、data.frameを使って検定結果をデータフレームに格納し、write.tableやwrite.csvを用いて外部ファイルへ保存しましょう。以下、csvファイルへの保存例です。
1 2 | resultTable <- data.frame(stat = stat, p.value = pValue, null.value = nullValue, row.names = NULL) write.csv(resultTable , "ウィルコクソンの順位和検定.csv", row.names = FALSE) |
ポチップ
2標本検定
次に、2標本のウィルコクソンの順位和検定（マン・ホイットニーのU検定）の実行例について見ていきます。
ggplot2で作成したプロットを結合させるために、次のgirdExtraを用います。インストールしといてください。
1 2 | install.packages("gridExtra") library(gridExtra) |
データセットとして、次のCO2を使います。
1 2 | #2標本検定 dataset <- CO2 |
CO2は uptake（二酸化炭素の吸収量）やTreatment（凍らしたか、凍らしてないか）などの要素から成ります。
1 2 3 4 5 6 7 8 | > head(dataset) Plant Type Treatment conc uptake 1 Qn1 Quebec nonchilled 95 16.0 2 Qn1 Quebec nonchilled 175 30.4 3 Qn1 Quebec nonchilled 250 34.8 4 Qn1 Quebec nonchilled 350 37.2 5 Qn1 Quebec nonchilled 500 35.3 6 Qn1 Quebec nonchilled 675 39.2 |
uptakeをTreatmentのlevels別にヒストグラムと箱ひげ図をプロットすると次のようになります。
1 2 3 4 5 6 | plot_hist <- ggplot(dataset, aes(x = uptake, color = Treatment, fill = Treatment)) + geom_histogram(position = "identity", alpha = 0.75, bins = 10) plot_box <- ggplot(dataset, aes(x = uptake, color = Treatment, fill = Treatment)) + geom_boxplot(alpha = 0.75) grid.arrange(plot_hist, plot_box) |
また、次を実行することで、下の図のようなuptakeの順位に関するグラフを作成することができます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | dataset_plot <- cbind(dataset, rank = rank(dataset$uptake, ties.method = "first"))[order(dataset$uptake), ][seq_len(nrow(dataset)) %% 5 == 0, ] ggplot(dataset) + geom_point(aes(x = uptake, y = 1, color = Treatment)) + geom_hline(aes(yintercept = 1)) + geom_segment(data = dataset_plot, aes(x = uptake, y = 1, color = Treatment), xend = dataset_plot$uptake, yend = 1.15) + geom_text(data = dataset_plot, aes(x = uptake, y = 1.2, color = Treatment, label = rank), size = 3.75) + geom_text(data = data.frame(x = rep(10, 2), y = c(2.7, 2.5), Treatment = as.factor(c("nonchilled", "chilled"))), aes(x = x, y = y, color = Treatment, label = paste0("rank sum: ", c(sum(rank(dataset$uptake, ties.method = "first" )[dataset$Treatment == "nonchilled"]), sum(rank(dataset$uptake, ties.method = "first" )[dataset$Treatment == "chilled"]))))) + scale_y_continuous("", limits = c(0.5, 3), breaks = NULL) |
ヒストグラムと上の順位のプロットから、Treatmentlがnonchilledであるときはデータの順位が高い傾向があり、chilledであるとき順位が低い傾向があるのが見てとれます。
順位に関するプロットを描いたのは、ウィルコクソンの順位和検定（マン・ホイットニーのU検定）の検定統計量が順位和に基づいているからです。
この後、「Treatmentに関してuptakeは同じである」仮説検定を実行していきますが、この帰無仮説の下では、Treatmentがnonchilledまたはchilledのグループのuptakeの順位和がその期待値より大きく下回ったり、上回ったりすることはないはずです。
「Treatmentによってuptakeは同じである」という仮説検定を実行したいと思います。この2標本の検定は、次のようにwilcox.testの引数xとyを指定するか、formulaとdataを指定します。2通りの書き方があります。
1 2 3 4 5 | testResult <- wilcox.test(uptake ~ Treatment, data = dataset) uptake_nonchilled <- dataset$uptake[dataset$Treatment == "nonchilled"] uptake_chilled <- dataset$uptake[dataset$Treatment == "chilled"] wilcox.test(uptake_nonchilled, uptake_chilled) #別の書き方 |
ウィルコクソンの順位和検定の結果は次のようになります。
1 2 3 4 5 6 7 | > testResult Wilcoxon rank sum test with continuity correction data: uptake by Treatment W = 1187.5, p-value = 0.006358 alternative hypothesis: true location shift is not equal to 0 |
「p-value = 0.006358 < 0.05」より、有意水準0.05の下で帰無仮説は棄却され、Treatmentによってuptakeは同じであるとはいえないことが分かりました。
検定結果を取得するには、1標本のきと同様にtestResultのうしろに$を付けます。
1 2 3 | stat <- testResult$statistic #検定統計量 pValue <- testResult$p.value #p値 nullValue <- testResult$null.value #帰無仮説の下での母平均 |
これらの検定結果をまとめるときは次のように関数data.frameを使ってデータフレームに格納すると便利です。
1 2 3 4 | > resultTable <- data.frame(stat = stat, p.value = pValue, null.value = nullValue, row.names = NULL) > resultTable stat p.value null.value 1 1187.5 0.006358173 0 |
まとめ
R言語でウィルコクソンの順位和検定やマン・ホイットニーのU検定を実行する関数や実行例を見てきました。
関数wilcox.testを用いることでウィルコクソンの順位和検定を行うことが可能です。
関数wilcox.testは1標本検定と2標本検定のどちらにも対応しています。
2標本検定の場合は、マン・ホイットニーのU検定の結果と一致します。
