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の符号順位決定についての詳細な説明は省くが以下のページで詳しく述べられている。
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が大きいにもかかわらず差があると出ている。
[Pが大きいにもかかわらず差があると出ている。]
間違えていますよ。
[alternative hypothesis: true location shift is not equal to 0]は，対立仮説はこれこれですよというのを書いているだけで，「対立仮説が正しい」と言っているわけではない。
また，階級ごとに平均値を求め？それを男女比較しているわけですが，それはおかしいでしょう。
t検定に対応するためには，元のデータをつかって wilcox 検定をやらないといけないでしょう。つまり，wilcox.test(spnbmd ~ gender, bone) とする。