【R言語】得られたデータから効果量を計算する（ノンパラメトリック二群比較編）

こんにちは。プログラミング初心者のえいこです。

最近はもっぱら効果量とか検定力とかサンプルサイズについて勉強しています。

前回の記事では、パラメトリックな二群比較の場合の効果量の計算のコードを書きました。

今回はノンパラメトリックだった場合の効果量の計算をやっていきたいと思います。

パラメトリックとノンパラメトリックの違いについては、こちらの記事にまとめています。

t検定やTukey法など耳なじみのある検定は、ほとんどがパラメトリックな検定です。

パラメトリック検定がよく使われる理由としては、
・実験を始める前にデータの正規性を前提にした方がサンプルサイズなどを決めやすい。
・直感的にわかりやすい
・統計的有意差を出しやすい

などがあるのではないでしょうか？

ノンパラメトリック検定って意外と使えるんじゃないかと思っている今日この頃。

その理由は、
・データの正規性を確認しなくても使える
・もし正規分布だったとしてもノンパラメトリック検定だからダメという理由はない

汎用範囲が広いのがノンパラメトリック検定の特長です。

とりあえずはノンパラメトリック検定をしておけばOKというのが良いところ。

耳なじみがない人が多いので、それってどんな検定なの？と言われてた時に答えられるようにしておいた方が良いでしょう。（私も含めて）

ノンパラメトリックな二群検定って？

検定の初歩は２つの群を比べること。

パラメトリックな検定はt検定ですが、ノンパラメトリックの場合に一般的に使われるのはMann-WhitneyのU検定（ウィルコクソンの順位和検定も実質的には同じ方法というWikipediaの記事もあります）。

Mann-WhitneyのU検定はこちらの記事でもやっていて、デフォルトのRの関数である"wilcox.test()"で検定しました。

少し概念的なことをサクッとまとめておくと（参考にしたサイトはこちら）、

U検定のUというのはどういうものかと言うと、帰無仮説の下ではその分布がわかっている統計量で、

AとBの二つの群を比べたときに

（Aの方が値が大きい割合）×（AとBのサンプル数の積）
（Bの方が値が大きい割合）×（AとBのサンプル数の積）

を比較したうちの小さい方のこと。

愚直にUを計算する方法がこちら、（こちらのサイトを参考にしました）

n1=3 #groupAのサンプル数
n2=3　#groupBのサンプル数
(r=rank(data$value))　#データの上から順位をつける
(W1=sum(r[1:n1])-n1*(n1+1)/2) #groupAの方が大きい割合
(W2=sum(r[n1+1:n2])-n2*(n2+1)/2) #groupBの方が大きい割合
(W=max(W1,W2)) #W1とW2で大きい方をWとする
(U=n1*n2-W) #Uの値を計算する

を実行すると、

> n1=3 #groupAのサンプル数
> n2=3　#groupBのサンプル数
> (r=rank(data$value))　#データの上から順位をつける
[1] 3 2 1 4 6 5
> (W1=sum(r[1:n1])-n1*(n1+1)/2) #groupAの方が大きい割合
[1] 0
> (W2=sum(r[n1+1:n2])-n2*(n2+1)/2) #groupBの方が大きい割合
[1] 9
> (W=max(W1,W2))　#W1とW2で大きい方をWとする
[1] 9
> (U=n1*n2-W) #Uの値を計算する
[1] 0

と計算していくと、U=0でした。

以前、Rで計算した時の結果を引っ張ってきましょう。

Wilcoxon rank sum test

data:  data$value by data$treatment
W = 0, p-value = 0.1
alternative hypothesis: true location shift is not equal to 0

Wで示されていますが、W=0となっています。

（そういえば「Wilcoxの順位和検定」と一番上に書いてありますね）

ノンパラメトリック検定の時の効果量の計算

ノンパラメトリックの時の効果量の計算方法はどうやら2パターンあるようです。

Rでの効果量の計算方法

・Cliffのdelta
・ウィルコクソン順位和検定の場合はZ/√N

Cliffのdeltaについて、先に計算してみましょう。

Cliffのdeltaを計算する？

Cliffのdeltaは効果量の指標の一つで（Choenのdみたいなもの）ノンパラメトリックだったらどんな検定でも、Cliffのdeltaで効果量が示せるようです。（参考にしたサイトはこちら）

効果量の考え方は、AとBの２つの群のデータを総当たりで大小比較をして、

（Aの方が大きい割合）-（Bの方が大きい割合）

を計算したもの。（実質的には-1から1の範囲に入る）

Cliffのdeltaの計算はデフォルトのRには関数としてないので、"effsize"というパッケージをインストールして使います。

install.packages("effsize")

として、１回インストールしておきます。

で、こんな感じにコードを打って実行すると...

library(effsize)
cliff.delta(data$value ~ data$treatment)

Cliff's Delta

delta estimate: -1 (large)
95 percent confidence interval:
    lower      upper 
-1.0000000 -0.1438497 

警告メッセージ: 
cliff.delta.default(x, y, conf.level, use.unbiased, use.normal,  で: 
 The samples are fully disjoint, using approximate Confidence Interval estimation

Cliffのdeltaの値は-1で(Large)という結果でした。

警告メッセージが出ていて、サンプルは不連続でバラバラですと書かれていますが...（あまり気にしなくても良いんですかね？）

Cliffのdeltaの効果量の大きさの指標が、こちら。（Romano (2006)の論文を参考にしたようです）

Z/√Nを計算する

「ノンパラメトリックの効果量を計算する」として日本理学療法士学会のサイトでは、Z/√Nで計算する方法を紹介しています。

私の第一印象は、「効果量はサンプルサイズ（N）に依存しないのでは？」ということ。

とりあえず計算はしてみましょう。

まず検定総計量であるZの値が必要なのですが...

Rの"wilcox.test()"関数では、Zの値が出力されません。

どうやら、"coin"というパッケージを使うとZ値が出力されるようなので、"coin"パッケージをインストールします。（参考にしたサイトはこちら）

install.packages("coin")
library(coin)

で、検定をしてます。

coinパッケージでのMann-WhitneyのU検定は"wilcox_test()"と.ではなく_なので注意が必要です。

> library(coin)
要求されたパッケージ survival をロード中です 
> wilcox_test(data$value ~ data$treatment)

	Asymptotic Wilcoxon-Mann-Whitney Test

data:  data$value by data$treatment (placebo, X)
Z = -1.964, p-value = 0.04953
alternative hypothesis: true mu is not equal to 0

と、Zの値が出てきました。

このZの値を取り出して、計算してみようと思います。

検定結果から数値を取り出す方法は、基本的にはこちらの記事と同じようなことをしています。

書いたコードはこちら↓

wt<-wilcox_test(data$value ~ data$treatment, distribution="exact")
wt
z1<-statistic(wt) #Z（検定統計量）
z1
#効果量を計算する
effect.size<- abs(z1/sqrt(length(data$treatment)))
effect.size

> z1<-statistic(wt) #Z（検定統計量）
> z1
[1] -1.963961
> #効果量を計算する
> effect.size<- abs(z1/sqrt(length(data$treatment)))
> effect.size
[1] 0.8017837

で、効果量は0.8と出たので結果は（large）。描いたグラフの横にr=0.80と書けますね。

最終的にはCliffのdeltaと同じlargeという結果に落ち着きました。

次は、多群比較の場合の効果量について勉強してみようと思います。

それでは、また！