Wilcoxonの符号付順位決定

boneのspnbmdを使って年齢別に男女差があるかを確認してみる。検定ではWilcoxonの符号付順位決定を利用する。

boneは特定個人について骨密度を計測している。データにはIDと骨密度以外に性別と計測時の年齢が登録されている。

> str(bone)
'data.frame':   485 obs. of  5 variables:
 $ idnum  : int  1 1 1 2 2 2 3 3 3 4 ...
 $ age    : num  11.7 12.7 13.8 13.2 14.3 ...
 $ gender : Factor w/ 2 levels "female","male": 2 2 2 2 2 2 2 2 2 1 ...
 $ spnbmd : num  0.01808 0.06011 0.00586 0.01026 0.21053 ...

年齢で層別化して男女別に平均値をプロットしてみると以下のようになる。

この図についてこちらで詳細な手順について記載している。

パラメトリックな検定は正規分布と等分散を前提とするが、boneデータを男女別に年齢で分けてみてみると、正規分布とはいえないので平均値の検定であるt検定は利用できない。しかし男女別に違いはありそうなので、統計的にその違いをだしてみたい。

ノンパラメトリックな検定方法はいくつかあるがここではWilcoxonの符号付順位決定を利用してみる。Wilcoxonの符号順位決定についての詳細な説明は省くが以下のページで詳しく述べられている。

• 符号検定とWilcoxonの符号付き順位検定
• ウィルコクソンの符号付順位和検定

RにおけるWilcoxonの検定はwilcox.test()で実施する。これについては下記のページで説明されている。
ウィルコクソンの符号付順位和検定

boneについてはこのままでは検定に利用できないのでフォーマットを変える。変換の目的は以下の２つである。

• 男女別年齢別に平均値を出すこと
• 年齢別をkeyとして男女のspnbmdを横展開すること

まずcutを利用して年齢別に分ける。

> bone$age_cat <- cut(bone$age, breaks=seq(9, 26, 1), right=FALSE, ordered_result=TRUE)
> head(bone)
  idnum   age gender      spnbmd age_cat
1     1 11.70   male 0.018080670 [11,12)
2     1 12.70   male 0.060109290 [12,13)
3     1 13.75   male 0.005857545 [13,14)
4     2 13.25   male 0.010263930 [13,14)
5     2 14.30   male 0.210526300 [14,15)
6     2 15.30   male 0.040843210 [15,16)

次にaggragate()を利用して年齢層別(age_cat), 性別(gender)で平均値を計算する。同時にフィールド名を変更する。このフィールド名の変更は次のdcast()の前処理となる。

bone01 <- aggregate(spnbmd~age_cat + gender, bone, mean)
colnames(bone01) <- c("age_cat", "variable", "value")
> head(bone01)
  age_cat variable      value
1  [9,10)   female 0.05042776
2 [10,11)   female 0.06930900
3 [11,12)   female 0.09346837
4 [12,13)   female 0.09673106
5 [13,14)   female 0.07934586
6 [14,15)   female 0.03964870
> tail(bone01)
   age_cat variable        value
29 [20,21)     male  0.008247127
30 [21,22)     male -0.010319843
31 [22,23)     male  0.006744858
32 [23,24)     male -0.002232583
33 [24,25)     male  0.002787990
34 [25,26)     male  0.009482647

上記では年齢別性別をkeyとしてspnbmdの平均値が計算されている。これでも検定に利用できるが扱いやすくするために年齢をkeyとしてspnbmdの男女を横展開する。横展開しるためにはreshape2パッケージのdcast()を利用する。

dcastではkeyとなるフィールドと横展開するフィールドをパラメータに指定できるが、紐づく値のフィールドについては”value”という名前であることを前提としている。このため前のステップで以下のようにspnbmdについては名前を変えている。

colnames(bone01) <- c("age_cat", "variable", "value")

dcastを利用して展開する。

bone02 <- dcast(bone01, age_cat ~ variable)
head(bone02)
> head(bone02)
  age_cat     female       male
1  [9,10) 0.05042776 0.03006022
2 [10,11) 0.06930900 0.03806144
3 [11,12) 0.09346837 0.03981337
4 [12,13) 0.09673106 0.06545330
5 [13,14) 0.07934586 0.07054018
6 [14,15) 0.03964870 0.08225971

Wilcox.test()を実行する。

> wilcox.test(bone02[,2], bone02[,3])

        Wilcoxon rank sum test

data:  bone02[, 2] and bone02[, 3] 
W = 138, p-value = 0.8384
alternative hypothesis: true location shift is not equal to 0

Pが大きいにもかかわらず差があると出ている。