Rによるカテゴリカルデータ分析事例(1) {coin}パッケージによるウィルコクソンの順位和検定とクラスカル・ウォリス検定

先日（2015年12月26日土曜日）に開催された Kobe.R 第23回で発表を行いました。当日はちょっと風邪気味だったため、あまり準備ができず、テキトーな内容になってしまいました…。そのため、ここで補足をしたいと思います。

kobexr.doorkeeper.jp

　
発表タイトルは「Rによるカテゴリカルデータ分析事例」。現在、私が分析を進めているデータを使った内容です。扱ったデータの構成は

• 性別（0が男性、1が女性の名義尺度）
• 年齢層（1が20代、2が30代、3が40代、4が50代、5が60代の名義尺度）
• あるサービスを受けた満足度（1～7の順序尺度）
• 変数x（1～5の順序尺度）

というもの。実際のデータセットはもっと多くて50個ほど変数があります。サンプルサイズは500です。今回使用するデータは↓においてあるので、興味がある方はやってみてください。 https://app.box.com/s/0jikm4cfeu4iopdhthfkyqlrw4z7ibed

　
分析によって調べたいことは次の3つです。

1. 満足度および変数xは性別によって異なるか。
2. 満足度および変数xは年齢層によって異なるか。
3. 変数xは満足度に正の影響を与えているか（変数xが高いほど、満足度も高いか）。

　
今回は１と２を扱います。３は次回のエントリーに回します。では、やってみます。有意水準は5%とします。

０．使用するデータの準備と確認

データの構成は次のようになっています。

> d <- blogyou.data.20151229
> str(d)
'data.frame':   500 obs. of  4 variables:
 $ gender : int  0 0 0 0 0 1 1 0 0 1 ...
 $ age_seg: int  3 4 4 4 3 1 2 4 3 3 ...
 $ cs     : int  5 7 4 6 6 6 5 6 7 5 ...
 $ x      : int  3 2 1 1 3 4 2 5 2 3 ...

　
性別と年齢層を因子型に変えて、もう1度構成を確認します。

d$gender <- as.factor(d$gender)
d$age_seg <- as.factor(d$age_seg)
> str(d)
'data.frame':   500 obs. of  4 variables:
 $ gender : Factor w/ 2 levels "0","1": 1 1 1 1 1 2 2 1 1 2 ...
 $ age_seg: Factor w/ 5 levels "1","2","3","4",..: 3 4 4 4 3 1 2 4 3 3 ...
 $ cs     : int  5 7 4 6 6 6 5 6 7 5 ...
 $ x      : int  3 2 1 1 3 4 2 5 2 3 ...

１．満足度および変数xは性別によって異なるか？

性別は2値の名義尺度、満足度と変数xは順序尺度です。まずはクロス表を作ってデータを眺めてみます。

　
●満足度と性別の関係

> xtabs(~gender + cs, data=d)

	1	2	3	4	5	6	7
男性(0)	1	3	6	38	66	59	35
女性(1)	4	1	13	41	90	80	63

　
ん～、どうでしょう？　女性のほうが高い満足度にある人が多いような…。図示してみます。

> library(lattice)
> histogram(~factor(cs)|gender, data = d)
## cs は離散値のため、factor() で因子型に直しています。以下、x も同様です。

では、実際に性別によって満足度の差があるか検定します。上図より正規分布に従っているようには見えないので、ノンパラメトリック検定のウィルコクソンの順位和検定（マン・ホイットニーのU検定）を行います。帰無仮説は「性別によって差がない（差が０）」、対立仮説が「性別によって差がある（差が０でない）」です。

　
さて、ここで{coin}パッケージを使います。ウィルコクソンの順位和検定を行うには、wilcox.test や {exactRankTests} の wilcox.exact を使うのですが、{exactRankTests}をインストールして使おうとしたら「{coin}を使え！」とコンソール画面で怒られたので、{coin} を使うことにします。

## {coin}パッケージがない場合は、インストールしてください。
library(coin)
> kekka <- wilcox_test(cs ~ gender, distribution = "exact", data = d)
## distribution = "exact"とすると、パソコンによっては結果が出るまで時間がかかります。
> kekka
    Exact Wilcoxon-Mann-Whitney Test
data:  cs by gender (0, 1)
Z = -1.0953, p-value = 0.2737
alternative hypothesis: true mu is not equal to 0

　
帰無仮説は棄却できませんでした。したがって、性別によって満足度に違いがあるとは言えませんでした。次に、変数ｘと性別の関係と調べます。

　
●変数ｘと性別の関係
先ほどと同じ検定を行いますが、その前にクロス表とヒストグラムで確認します。

> xtabs(~gender + x, data =d)

	1	2	3	4	5
男性(0)	52	30	51	53	22
女性(1)	47	56	63	74	52

histogram(~factor(x)|gender, data =d)

女性は特に変数ｘが5の人が多いように見えます。検定します。

> kekka2 <- wilcox_test(x ~ gender, distribution = "exact", data = d)
> kekka2

    Exact Wilcoxon-Mann-Whitney Test

data:  x by gender (0, 1)
Z = -2.2005, p-value = 0.02769
alternative hypothesis: true mu is not equal to 0

変数x は性別によって違いがあることがわかりました。

２．満足度および変数xは年齢層によって異なるか。

年齢層は5つの水準がある名義変数です。年齢層と満足度および変数xの関係をクロス表とヒストグラムで確認します。

> xtabs(~age_seg + cs, data =d)

	1	2	3	4	5	6	7
20代(1)	0	1	4	18	32	25	20
30代(2)	4	2	3	18	34	22	17
40代(3)	0	1	6	19	30	25	19
50代(4)	0	0	0	14	31	37	18
60代(5)	1	0	6	10	29	30	24

> xtabs(~age_seg + x, data =d) 

	1	2	3	4	5
20代(1)	24	23	17	19	17
30代(2)	20	22	24	19	15
40代(3)	25	10	24	29	12
50代(4)	18	16	21	29	16
60代(5)	12	15	28	31	14

histogram(~factor(cs)|age_seg, data =d)
histogram(~factor(x)|age_seg, data =d)

どうも年齢層による違いは満足度、変数ｘともに見られそうにありません。検定します。水準が3つ以上なので、先ほどのウィルコクソンの順位和検定（マン・ホイットニーのU検定）ではなく、クラスカル・ウォリス検定を行います。

## kruskal_test() は {coin} パッケージにあります。
kekka3 <- kruskal_test(cs ~ age_seg, data = d)
kekka4 <- kruskal_test(x ~ age_seg, data = d)

　
満足度および変数xの年齢層による違いは次の結果から差があるとはいえませんでした。

> kekka3

    Approximative Kruskal-Wallis Test

data:  cs by age_seg (1, 2, 3, 4, 5)
chi-squared = 8.0373, p-value = 0.0895

> kekka4

    Asymptotic Kruskal-Wallis Test

data:  x by age_seg (1, 2, 3, 4, 5)
chi-squared = 5.4405, df = 4, p-value = 0.245

次回に続きます。今回、参考にしたものはこちらの2冊です。

カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)

人文・社会科学のためのカテゴリカル・データ解析入門