【統計】クラスカル・ウォリス検定(Kruskal-Wallis test)

クラスカル・ウォリス検定(Kruskal-Wallis test)について。

【目次】

 

クラスカル・ウォリス検定(KW検定)は3標本以上を評価するノンパラメトリックな手法。
Wilcoxonの順位和検定・マンホイットニーのU検定(以下、Wilcoxonの順位和検定等)は2標本。
 

帰無仮説、対立仮設


  • 帰無仮説  H_{0}: 全ての標本の確率分布に差がない(ズレていない)。
  • 対立仮設  H_{1}: 何れかの標本の確率分布に差がある(ズレている)。

Wilcoxonの順位和検定等と同じく、確率分布の位置のズレを見る。正規分布に依存しないのが特徴。

 

計算式等


順位和

統計量の算出前に、先ず、比較したい全標本の観測値に全体順位(Rank)をつけて、標本ごとの合計(順位和)を算出する。観測値にタイ(同一順位)がある場合は、順位の平均値をつける。この方法はWilcoxonの順位和検定等と同じ。

 

H統計量

\qquad H=\dfrac{12}{N\left( N+1\right) }\sum \left(\dfrac{Ri^{2}}{n_{i}}\right) - 3 \left( N+1 \right)

H統計量は近似的に  \chi ^{2} 分布( df = \left( k-1 \right))に従うため、 \chi^{2} 分布を使って検定を行う。

また、タイ(同一順位)の補正を行う場合、上で求めたH統計量にタイの補正値(D)を除算する。

\qquad D=1- \dfrac{ \sum \left( t_{i}^{3} - t_{i} \right) }{ \left( n-1 \right) n \left( n + 1\right) }

\qquad H_{adj} = \dfrac{H}{D}

Wilcoxonの順位和検定等と同じくタイ(同一順位)の補正式第2項の t が同一値の個数を示す。
例えば、{10,10,11,12,12,12} のデータがあるとき、10が2つ、11が1つ、12が3つになるため、それぞれ  t=2、t=1、t=3 となり、  t^{3}-t の計算で  2^{3}-2=6、1^{1}-1=0、3^{3}-3=24 になる。タイがない数値は解が0になり、タイにだけ働く。
なお、タイがまったくない場合は第2項の解全体が0になり、D=1-0=1になる。 H_{adj} = \dfrac{H}{D} = \dfrac{H}{1} = H で等価。
  

計算例


今回は、架空の4標本に違いがあるか評価する。

 A {10,10,11,12,12,12}
 B {12,13,14,17,19}
 C {12,14,15,19,20}
 D {16,17,18,20}

Excelでの計算:

f:id:cochineal19:20201218215746p:plain

Kruskal-Wallis test
Excelの簡単な説明。
・A列:標本のグループ名。
・B列:観測値。小さい順にならべている。
・C列:同値の数を上から数ている(セルC4:=IF(B3<>B4,1,D3+1))。
・D列:同値の最終行のみ表示している(セルD4:=IF(B4=B5,"",C4))。
・E列:分散の補正項用(セルE4:=IF(D4<>"",D4^3-D4,""))
・G列:標本Aの順位(セルG4:=IF($A4=G$3, RANK.AVG($B4,$B:$B,1), ""))
・H列:標本Bの順位(セルH4:=IF($A4=H$3, RANK.AVG($B4,$B:$B,1), ""))
・I列 :標本Cの順位(セルI4 :=IF($A4=I$3, RANK.AVG($B4,$B:$B,1), ""))
・J列:標本Dの順位(セルJ4:=IF($A4=J$3, RANK.AVG($B4,$B:$B,1), ""))

 

H統計量:

順位和は  R_{A}=24,\ R_{B}=57.5,\ R_{C}=65.5,\ R_{D}=63 だったため、 \dfrac{R^{2}}{n} はそれぞれ  A=\dfrac{24^{2}}{6}=96,\ B=\dfrac{57.5^{2}}{5}=661.25,\ C=\dfrac{65.5^{2}}{5}=858.05,\ D=\dfrac{63^{2}}{4}=992.25 になる(Excel L3:P7)。
全体数はN=20のため、公式通り、H統計量は  H=\dfrac{12}{20\times21}\times \left(96+661.25+858.05+992.25 \right) - (3 \times 21)=0.0286\times 2607.55 - 63 = 11.5014 となる(Excel L9:M15)。

 
タイの補正:

同一値 t の計算( \sum t_{i}^{3}-t_{i})は 150 になる(Excel D:E, P12)。
(今回は 12 が5つあり、 5^{3}-5=125-5=120 と大きな値がある。)
公式通り、 D=1- \dfrac{150}{19 \times 20 \times 21} = 0.9812 となり、このDで除算して  H_{adj}= \dfrac{11.5014}{0.9812}=11.7218 が得られる(Excel P13:P15)。

 

p値:

計算式等で述べたとおり、H統計量は  \chi^{2} 分布( df=k-1)に従うため、次のようになる。
\quad \verb|タイ補正なし:|\ p=0.0093 \ \left( H=11.5014,\ df=3\right)
\quad \verb|タイ補正あり:|\ p=0.0084 \ \left( H_{adj}=11.7218,\ df=3\right)
今回の架空データでは、標本間に差がある(p<0.01)。(タイ補正有無問わず)

 

Rでの実行: 

> ads <- data.frame(
+   AVAL    = c(10,10,11,12,12,12,12,12,13,14,14,15,16,17,17,18,19,19,20,20)
+   ,TRT01P = c("A","A","A","A","B","A","A","C","B","B","C","C","D","B","D","D","C","B","C","D"))
> ads$TRT01PF <- factor(ads$TRT01P)
> kruskal_test(AVAL ~ TRT01PF, data=ads, method="exact")

	Asymptotic Kruskal-Wallis Test

data:  AVAL by TRT01PF (A, B, C, D)
chi-squared = 11.722, df = 3, p-value = 0.0084

> kruskal.test(x=ads$AVAL, g=ads$TRT01PF)

	Kruskal-Wallis rank sum test

data:  ads$AVAL and ads$TRT01PF
Kruskal-Wallis chi-squared = 11.722, df = 3, p-value = 0.0084 

デフォルトで使える kruskal.test、coinパッケージのkruskal_testの2つ関数がありますが、どちらもタイ補正ありの結果を得る。
 

SASでの実行:

タイ補正ありの結果が得られる。

data ads;
  do AVAL=10,10,11,12,12,12; TRT01PN=1; output; end;
  do AVAL=12,13,14,17,19; TRT01PN=2; output; end;
  do AVAL=12,14,15,19,20; TRT01PN=3; output; end;
  do AVAL=16,17,18,20; TRT01PN=4; output; end;
run;
proc sort data=ads: by AVAL TRT01PN; run;
proc npar1way data=ads wilcoxon;
  class TRT01PN;
  var   AVAL;
  output out=OutputDS wilcoxon;
run;

f:id:cochineal19:20201219144836p:plain

 

プログラムコード


■ Rのコード

#-- coinパッケージ
library(coin)
kruskal_test(AVAL ~ Group, data=ADS, method="exact")

#-- デフォルトのstatsパッケージ kruskal.test(x=ADS$AVAL, g=ADS$Group)

 

SASのコード

proc npar1way data=[InputDS] wilcoxon;
  class Group;
  var   AVAL;
  output out=[OutputDS] wilcoxon; /*_KW_,P_KWを取得。*/
run;

 

Pythonのコード

整備中

 

2標本でもKW検定


2標本のデータについて、Wilcoxonの順位和検定とKW検定を実行すると、どちらも同じ p値が得られる。

> ads <- data.frame(
+   TRT01P <- c("A","A","A","A","A","A","A","A","B","B","B","A","A","B","B","A","B","B","B","B")
+   ,AVAL  <- c(4,9,10,11,12,13,13,15,16,16,17,18,18,19,19,20,22,22,23,25)
+ )
> ads$TRT01PF <- factor(ads$TRT01P)
> wilcox_test(AVAL ~ TRT01PF, data=ads, method="exact")

	Asymptotic Wilcoxon-Mann-Whitney Test

data:  AVAL by TRT01PF (A, B)
Z = -2.9305, p-value = 0.003384
alternative hypothesis: true mu is not equal to 0

> kruskal_test(AVAL ~ TRT01PF, data=ads)

	Asymptotic Kruskal-Wallis Test

data:  AVAL by TRT01PF (A, B)
chi-squared = 8.5878, df = 1, p-value = 0.003384  

参考


https://www.dataanalytics.org.uk/adjustment-for-tied-ranks-in-the-kruskal-wallis-test/

https://www.lexjansen.com/wuss/2018/26_Final_Paper_PDF.pdf

 


サイトマップ

cochineal19.hatenablog.com

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。