確率・統計

(14) 分散分析法(ANOVA)

「二標本の検定(Two Sample Test)」の例として、以前「F-検定(F-Test)」を紹介しました。二つの母集団があり、その分散が等しい場合、母集団の平均は等しいかどうかを F-分布を使って検定することができるのですが、これをさらに多くの母集団に対して一度にまとめて検定する手法として「分散分析法(ANalysis Of VAriance ; ANOVA)」というものがあります。ここでは、ANOVA を中心に紹介をしたいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 一元配置分散分析法(One-way ANalysis Of VAriance ; One-way ANOVA)

k 個の母集団 (例えば装置ごとの製品処理時間やクラス別の学力テスト結果など) があり、各集団を Ai ( i = 1, 2, ..., k ) で表します。Ai から Ni 個の標本を抽出したとして、それを xij ( j = 1, 2, ... Ni ) とします。また、母集団 Ai は正規分布 N( μi, σ2 ) に従う(このとき、各母集団の分散は等しいものと仮定しています)とします。このとき、

μ = Σi{1→k}( niμi ) / N

但し、N = Σi{1→k}( ni )

とすれば、μ は抽出した標本全体に対する平均の期待値を表すことになり、

di = μi - μ

とすると、di は母集団全てにおける平均と各母集団の平均との差分を表すことになります。複数の母集団が表す因子(装置ごとの製品処理時間やクラス別の学力テスト結果など)を A で表すとき、diAi における因子 A の「要因効果」と呼ばれます。

ここで、次の三つの平方和について考えてみます。

ST = Σi{1→k}( Σj{1→ni}( ( xij - m )2 ) )

SC = Σi{1→k}( ni( mi - m )2 )

SE = Σi{1→k}( Σj{1→ni}( ( xij - mi )2 ) )

但し、m = Σi{1→k}( Σj{1→ni}( xij ) ) / N ; mi = Σj{1→ni}( xij ) / ni

ST は、全体の平均 m との差分の平方和であり「全体の変動」を表します。同様に考えて、SC は各集団間での変動を表し、SE は集団内での変動を表しています。ST の式を変形すると

ST=Σi{1→k}( Σj{1→ni}( ( xij - m )2 ) )
=Σi{1→k}( Σj{1→ni}( { ( xij - mi ) + ( mi - m ) }2 ) )
=Σi{1→k}( Σj{1→ni}( ( xij - mi )2 + 2( xij - mi )( mi - m ) + ( mi - m )2 ) )
=Σi{1→k}( Σj{1→ni}( ( xij - mi )2 ) ) + 2Σi{1→k}( ( mi - m )Σj{1→ni}( xij - mi ) ) + Σi{1→k}( ni( mi - m )2 )

と表すことができます。最後の式において、第一項は SE であり、第三項は SC です。残る第二項は、Σj{1→ni}( xij - mi ) がゼロになるので無視することができて、結局

ST = SE + SC

という関係式が得られます。全体の変動が、各集団間の変動と集団内での変動に分解できることをこの式は表しています。SE と SC の期待値を求めると、

E[SE]=E[Σi{1→k}( Σj{1→ni}( ( xij - mi )2 ) )]
=Σi{1→k}( E[Σj{1→ni}( ( xij - mi )2 )] )
=Σi{1→k}( ( ni - 1 )σ2 )
=( N - k )σ2
E[SC]=E[Σi{1→k}( ni( mi - m )2 )]
=Σi{1→k}( niE[mi2 - 2mim + m2] )
=Σi{1→k}( niE[mi2] )- 2E[mΣi{1→k}( nimi )] + NE[m2]
=Σi{1→k}( niE[mi2] )- 2NE[m2] + NE[m2]
=Σi{1→k}( niE[mi2] )- NE[m2]

E[mi2] は、集団内の各要素が互いに独立ならば

E[mi2]=E[{ Σj{1→ni}( xij ) }2] / ni2
={ Σj{1→ni}( E[xij2] ) + 2Σj{1→ni}( Σl{j+1→ni}( E[xijxil] ) ) } / ni2
={ Σj{1→ni}( σ2 + μi2 ) + 2Σj{1→ni}( Σl{j+1→ni}( E[xij]E[xil] ) ) } / ni2
={ ni( σ2 + μi2 ) + 2Σj{1→ni}( Σl{j+1→ni}( μi2 ) ) } / ni2
={ ni( σ2 + μi2 ) + ni( ni - 1 )μi2 } / ni2
=( niσ2 + ni2μi2 ) / ni2
=σ2 / ni + μi2

と求められます。また、E[m2] は、

E[m2]=E[{ Σi{1→k}( nimi ) }2] / N2
=E[Σi{1→k}( ni2mi2 ) + Σi{1→k}( Σl{1→k;l≠i}( ninlmiml ) )] / N2
={ Σi{1→k}( ni2E[mi2] ) + Σi{1→k}( Σl{1→k;l≠i}( ninlE[miml] ) ) } / N2

で求められます。各集団間の要素も互いに独立ならば、

E[miml] = E[mi]E[ml] = μiμl

であり、

E[m2]={ Σi{1→k}( ni2( σ2 / ni + μi2 ) ) + Σi{1→k}( Σl{1→k;l≠i}( ninlE[miml] ) ) } / N2
={ σ2Σi{1→k}( ni ) + Σi{1→k}( ni2μi2 + Σl{1→k;l≠i}( ninlμiμl ) ) } / N2
=[ Nσ2 + Σi{1→k}( niμi{ niμi + Σl{1→k;l≠i}( nlμl ) } ) ] / N2
=[ Nσ2 + Σi{1→k}( niμi{ Σl{1→k}( nlμl ) } ) ] / N2
={ Nσ2 + Σi{1→k}( niμiNμ ) } / N2
=( Nσ2 + N2μ2 ) / N2
=( σ2 + Nμ2 ) / N

となるので、

E[SC]=Σi{1→k}( ni( σ2 / ni + μi2 ) ) - ( σ2 + Nμ2 )
=( k - 1 )σ2 + Σi{1→k}( niμi2 ) - Nμ2

が求める解になります。ところで

Σi{1→k}( ni( μi - μ )2 )=Σi{1→k}( niμi2 - 2niμiμ + niμ2 )
=Σi{1→k}( niμi2 ) - 2μΣi{1→k}( niμi ) + μ2Σi{1→k}( ni )
=Σi{1→k}( niμi2 ) - 2μNμ + μ2N
=Σi{1→k}( niμi2 ) - Nμ2

なので、結局 E[SC] は

E[SC] = ( k - 1 )σ2 + Σi{1→k}( ni( μi - μ )2 )

と表すこともできます。E[SE] と E[SC] の二つの和から、E[ST] は

E[ST] = E[SE] + E[SC] = ( N - 1 )σ2 + Σi{1→k}( ni( μi - μ )2 )

となり、それぞれの平方和の期待値が得られたことになります。未知数 σ2 は SE と SC の期待値を使って

σ2 = E[SE] / ( N - k ) = { E[SC] - Σi{1→k}( ni( μi - μ )2 ) } / ( k - 1 )

になるので、各集団の平均値 μi が全体の平均 μ より外れるほど E[SC] / ( k - 1 ) の値は E[SE] / ( N - k ) よりも大きくなります。そこで、

F = { SC / ( k - 1 ) } / { SE / ( N - k ) }

としたとき、F は各集団間の平均に差がなければ 1 に近い値となり、差が大きくなるほど F は大きくなっていきます。これを利用して平均値に差がないか検定する手法を「(一元配置)分散分析法((One-way) ANalysis Of VAriance ; (One-way) ANOVA)」といいます。

Σj{1→ni}( ( xij - mi )2 ) / σ2 は自由度 ni - 1 の χ2-分布に従います。

SE = Σi{1→k}( Σj{1→ni}( ( xij - mi )2 ) )

より、SE / σ2 は自由度 ni - 1 のχ2-分布に従う確率変数の和の形で表されるので、自由度 Σi{1→k}( ni - 1 ) = N - k の χ2-分布に従うことになります(*1-1)。また、mi は正規分布 N( μi, σ2 / ni ) に従うので、もし全ての集団 Ai の平均値が等しいとすれば、その平均値を μ としたとき ( mi - μ ) / ( σ2 / ni )1/2 = √ni( mi - μ ) / σ は正規分布 N( 0, 1 ) に従います(*1-2)。よって、Σi{1→k}( { √ni( mi - μ ) / σ }2 ) は自由度 k の χ2分布に従いますが、μ の代わりに mi の標本平均 m に置き換えるとこの式は SC / σ2 と等しくなり、このとき自由度を一つ落として k - 1 の χ2-分布に従うことになります(補足1)。これらの結果から、F は自由度 ( k - 1, N - k ) の F-分布に従い(*1-3)、これを利用して検定を行うことができます。この検定は、次に示すような「分散分析表(ANOVA table)」で表現されるのが一般的です。

平方和自由度不偏分散分散比
集団間SC = Σi{1→k}( ni( mi - m )2 )k - 1vC = SC / ( k - 1 )F = { SC / ( k - 1 ) } / { SE / ( N - k ) }
集団内SE = Σi{1→k}( Σj{1→ni}( ( xij - mi )2 ) )N - kvE = SE / ( N - k )
全体ST = Σi{1→k}( Σj{1→ni}( ( xij - m )2 ) )N - 1  

F が自由度 ( k - 1, N - k ) の F-分布に従うのは、各集団が正規分布に従うことが前提となりますが、正規分布でなくとも分散が等しければ、各集団の平均値 μi が全体の平均 μ より外れるほど E[SC] / ( k - 1 ) の値が E[SE] / ( N - k ) より大きくなる傾向にあるので、F の大きさだけを使って簡便的に調べるというような使い方もできます。但し、分散が等しいということが前提条件であることに変わりはないので、あらかじめ各集団の分散を求めて等しいと見てよいかを判断する必要があります。

ところで、ST = SE + SC の関係式から、全体の変動は集団間の変動と集団内の変動の二つに分けられることは前に説明したとおりです。そこで、全体の中で集団間の変動分がどの程度の値であるかを示す値として、下式のような計算を行います。

R2 = SC / ST

R2 の値が小さければ、集団間での変動も小さいことを意味し、変動の大部分は集団内の変動が占めることになるので、集団間における平均のばらつきは小さいということになり、逆に R2 の値が大きければ、集団間の変動が大きいことになるので、平均に差がある傾向は強くなります。この値は「寄与率(Contribution Ratio)」と呼ばれ、その平方根 R は「相関比(Correlation Ratio)」ともいいます。式から明らかなように、0 ≤ R2 ≤ 1 が成り立ちます。

*1-1)(6) 標本分布」の「カイ二乗分布に対する性質」を参照

*1-2)(5) 正規分布」の「3) 大数の法則(Law of Large Numbers)と中心極限定理(Central Limit Theorem)」を参照

*1-3)(6) 標本分布」の「F-分布に対する性質」を参照


ANOVA 検定を行うためのサンプル・プログラムを以下に示します。

/*
  oneWayANOVA : 一元配置分散分析(one-way ANOVA)

  const vector< vector<double> >& data : データ列
  double a : 危険率
  double threshold : binSearchで棄却域を求める時のしきい値

  戻り値 : True ... 仮説を棄却 , False ... 仮説を保留 または データがない
*/
bool oneWayANOVA( const vector< vector<double> >& data, double a, double threshold )
{
  cout << "***** One-way ANOVA *****" << endl << endl;

  unsigned int k = data.size(); // 集団の数
  if ( k < 2 ) {
    cerr << "At least two classes are necessary." << endl;
    return( false );
  }

  vector<double> ni( k ); // 各集団のデータ数
  vector<double> mi( k ); // 各集団の平均
  double m = 0;           // 全体の平均
  double n = 0;           // 全体のデータ数

  // 平均値の算出
  for ( unsigned int i = 0 ; i < k ; ++i ) {
    ni[i] = data[i].size();
    n += ni[i];
    if ( ni[i] == 0 ) {
      cerr << "No samples in data[" << i << "]." << endl;
      return( false );
    }
    mi[i] = sum( data[i] );
    m += mi[i];
    mi[i] /= (double)( ni[i] );
  }
  m /= (double)n;

  cout << "#Class(k) = " << k << " ; #Data(N) = " << n << endl;
  cout << "Mean. : " << m << endl;
  cout << "N     of class : ";
  for ( unsigned int i = 0 ; i < k - 1 ; ++i )
    cout << ni[i] << ", ";
  cout << ni[k-1] << endl;
  cout << "Mean. of class : ";
  for ( unsigned int i = 0 ; i < k - 1 ; ++i )
    cout << mi[i] << ", ";
  cout << mi[k-1] << endl << endl;

  // 変動の算出
  double sc = 0; // 各集団間の変動
  double se = 0; // 集団内の変動
  for ( unsigned int i = 0 ; i < k ; ++i ) {
    Deviation<double> dev( mi[i] );
    sc += (double)( ni[i] ) * pow( mi[i] - m, 2 );
    se += sum( data[i], dev );
  }

  // F値の算出
  double f = ( sc * (double)( n - k ) ) / ( se * (double)( k - 1 ) );

  // 結果出力
  cout << "<ANOVA Table>" << endl;
  cout << "\t\tSS\tDOF\tUnbiased Var." << endl;
  cout << "Between-class\t" << sc << "\t" << k - 1 << "\t" << sc / (double)( k - 1 ) << "\t\tF = " << f << endl;
  cout << "Intra-class\t" << se << "\t" << n - k << "\t" << se / (double)( n - k ) << endl;
  cout << "Total.\t\t" << sc + se << "\t" << n - 1 << endl << endl;

  cout << "R^2 = " << 100.0 * sc / ( sc + se ) << " %" << endl << endl; // 寄与率の出力

  // 検定
  FDistribution fDist( k - 1, n - k );
  double f0 = binSearch( fDist, 1.0 - a, threshold );
  cout << "F0 = " << f0 << endl << endl;

  double p = 1.0 - fDist.lower_p( f );
  cout << "p-value = " << p << endl << endl;

  if ( f0 < f ) {
    cout << "Null hypothesis ( average of all samples are same ) is rejected." << endl;
    return( true );
  }

  cout << "Null hypothesis ( average of all samples are same ) is accepted." << endl;
  return( false );
}

変数 sc は「各集団間の変動」、se は「集団内の変動」をそれぞれ表します。出力結果として、集団間の変動(Between-class)、集団内の変動(Intra-class)、全変動(Total) それぞれに対し、平方和(SS)、自由度(DOF)、不偏分散(Unbiased Var.) (但し、不偏分散については全変動に対する値は除いています)を求め、その右側に F 値を出力しています。前述の通り、このような出力内容は「分散分析表(ANOVA table)」と呼ばれ、前章の「重相関係数」でも利用されています。実は、重相関係数を利用した回帰の有意性検定も ANOVA 検定の一種とみなすことができます。

平均値や変動を求めるときに利用している関数 sum は、配列内の要素からなる和を求めるための関数で、第二引数として Deviation<double>型のオブジェクト dev を利用したものは、dev 構築時に渡した変数 m (これは平均値を表します)との差分の平方和を計算することができるようになっています。これらは、「(7) 標本の抽出と要約」にある「要約統計量算出用サンプル・プログラム」で紹介しています。

検定に利用している binSearch は、F-分布において P(0≤f≤f0) = 1.0 - a を満たす信頼区間の上限を二部探索を利用して求めるためのサブルーチンであり、「(8) 推定」の中の「平均値の区間推定用サンプル・プログラム」にて内容を確認することができます。また、F-分布を構築するクラスは「(6) 標本分布」の中の「サンプル・プログラム」にあります。


分散分析の例を以下に示します。キャベツを栽培するときの条件として、二種類の肥料を使い、肥料を用いない場合と併せて三通りの方法を用いたとします。その他の条件は全く同等として、収穫したキャベツの質量(g)を調べたところ、以下の結果が得られました。

肥料なし 肥料A  肥料B 
1149115111449
2136514961512
3151115691502
4143616321389
514291497 
6 1592 
平均144615501463
平方和13200160909594
不偏分散330032193198

不偏分散の値を見る限り、ほぼ等しいと見てよさそうなので、厳密な検定はせずにこのまま ANOVA 検定を行ってみましょう。まずは ST, SC, SE をそれぞれ求めると、次のようになります。

SC = 3.360 x 104 ; SE = 3.889 x 104 ; ST = 7.248 x 104

この結果を元に、分散分析表を作成すると、以下のようになります。

平方和自由度不偏分散分散比
集団間SC = 3.360 x 1042vC = 1.680 x 104F = 5.184
集団内SE = 3.889 x 10412vE = 3.241 x 103
全体ST = 7.248 x 10414  

F 値が 5.184 であるという結果は、平均には差があると考えてもよさそうに思えます。キャベツの質量が正規分布に従うと考えることができれば、自由度 ( 2, 12 ) の F-分布において F ≥ 5.184 になる確率は 0.02384 で、危険率 1 % の場合は保留となりますが、危険率 5 % とすれば(全ての集団において平均が等しいという)仮説は棄却され、いずれかの集団間において差があると判断されることになります。データを見る限りにおいては、「肥料 B」は「肥料なし」の場合と比べて平均値に差はないように見える反面、「肥料 A」は他の二つの集団よりも質量が大きく、肥料により収穫量は多くなると判断できそうに見えます。しかし、標本数が少ないため、本当に正規分布に従うと見て問題ないかどうか、疑問が残ります。以下は、それぞれの集団別に色分けしたヒストグラムです。各集団については、正規分布に従うと判断することはグラフを見てもできないことがわかります。実際にはもっと多くの標本を使って検定を行う必要があるでしょう。

集団別色分けヒストグラム

寄与率は、3.360 x 104 / 7.248 x 104 = 46.35 % になり、全体の半分程度が集団間の変動で占めるという結果になります。

以下は、サンプル・プログラムを使い、危険率 5% で検定を行った結果です。

***** One-way ANOVA *****

#Class(k) = 3 ; #Data(N) = 15
Mean. : 1492.07
N     of class : 5, 6, 4
Mean. of class : 1446.4, 1549.5, 1463

<ANOVA Table>
               SS      DOF     Unbiased Var.
Between-class  33598.2 2       16799.1          F = 5.18402
Intra-class    38886.7 12      3240.56
Total.         72484.9 14

R^2 = 46.352 %

F0 = 3.88529

p-value = 0.0238408

Null hypothesis ( average of all samples are same ) is rejected.

2) 二元配置分散分析法(Two-way ANOVA)

標本が二つの因子 A, B の影響を受けていると仮定します。例えば、キャベツの収量に影響する要因として肥料だけでなく品種もいっしょに問題にするような場合、前に示した手法は利用できません。そこで、次のようなモデルを考えます。

二つの要因 Ar ( r = 1, 2, ... l ) と Bc ( c = 1, 2, ... k ) に対し、その組み合わせ ( Ar, Bc ) ごとに一つの標本 xrc を割り当てるとします。各要因 Ar, Bc ごとの平均は

mr. = Σc{1→k}( xrc ) / k

m.c = Σr{1→l}( xrc ) / l

で表し、全体の平均を

m = Σr{1→l}( Σc{1→k}( xrc ) ) / kl

とすると、右に示した表のように表すことができます。

要因B1B2...Bk平均
A1x11x12...x1km1.
A2x21x22...x2km2.
:::...::
Alxl1xl2...xlkml.
平均m.1m.2...m.km

標本 xrc は全て、等しい分散を持つ正規分布 N( μrc, σ2 ) に従い、互いに独立であるとします。ここで μrc は、全体の平均 μ に対して、各因子によって変化する量を αr, βc として

μrc = μ + αr + βc

αr = μr. - μ = Σc{1→k}( μrc ) / k - μ

βc = μ.c - μ = Σr{1→l}( μrc ) / l - μ

で表せると仮定します。このとき、

Σr{1→l}( αr ) = Σr{1→l}( Σc{1→k}( μrc ) / k - μ ) = 0

Σc{1→k}( βc ) = Σc{1→k}( Σr{1→l}( μrc ) / l - μ ) = 0

となるので、μrc の総和は、μrc = μ + αr + βc の左辺と右辺それぞれから

Σr{1→l}( Σc{1→k}( μrc ) ) = klμ

Σr{1→l}( Σc{1→k}( μ + αr + βc ) ) = klμ

と求められ、等号は成り立つことが示されます。そこで、標本 xrc

xrc = μ + αr + βc + εrc

として、εrc が正規分布 N( 0, σ2 ) に従う確率変数であるとすれば、xrc は正規分布 N( μ + αr + βc, σ2 ) = N( μrc, σ2 ) に従うことになります。αr は、因子 Ar それぞれに、βc は因子 Bc それぞれに固有な量であり、各標本はこれらの量で加法的に表すことができると仮定しているわけです。

このモデルにおける変動量として、以下の値を定義します。

ST=Σr{1→l}( Σc{1→k}( ( xrc - m )2 ) )
SR=r{1→l}( ( mr. - m )2 )
SC=c{1→k}( ( m.c - m )2 )
SE=Σr{1→l}( Σc{1→k}( { ( xrc - m ) - ( mr. - m ) - ( m.c - m ) }2 ) )
=Σr{1→l}( Σc{1→k}( ( xrc - mr. - m.c + m )2 ) )

ST は「全体の変動」を表します。SR は各行の平均 mr. と全体の平均 m との差の平方和なので「行間変動」と呼ばれ、列間の変動とは無関係な変動量になります。同様な考え方から、SC は「列間変動」になります。最後の SE は、各標本 xrc から、全体の平均 m と、それぞれの行・列の平均と全体平均の差 ( mr. - m ) , ( m.c - m ) を除外した値の平方和であり、mr., m.c のいずれとも無関係な量になるので「誤差変動」と呼ばれます。モデル式に当てはめると、SR は αr の、SC は βc の平方和であり、SE は標本から μ と αr と βc を除外した値、すなわち εrc の平方和を表す量になると考えることができます。

SE を展開すると、

SE=Σr{1→l}( Σc{1→k}( ( xrc - mr. - m.c + m )2 ) )
=Σr{1→l}( Σc{1→k}( { ( xrc - m ) - ( mr. - m ) - ( m.c - m ) }2 ) )
=Σr{1→l}( Σc{1→k}( ( xrc - m )2 + ( mr. - m )2 + ( m.c - m )2
 - 2( xrc - m )( mr. - m ) - 2( xrc - m )( m.c - m ) + 2( mr. - m )( m.c - m ) ) )
=ST + SR + SC
 - 2Σr{1→l}( Σc{1→k}( ( xrc - m )( mr. - m ) + ( xrc - m )( m.c - m ) - ( mr. - m )( m.c - m ) ) )

になります。ST, SR, SC 以外の各項を計算すると、

Σr{1→l}( Σc{1→k}( ( xrc - m )( mr. - m ) ) )=Σr{1→l}( ( mr. - m )Σc{1→k}( xrc - m ) )
=Σr{1→l}( ( mr. - m )・k( mr. - m ) )
=r{1→l}( ( mr. - m )2 ) = SR
Σr{1→l}( Σc{1→k}( ( xrc - m )( m.c - m ) ) )=c{1→k}( ( m.c - m )2 ) = SC
Σr{1→l}( Σc{1→k}( ( mr. - m )( m.c - m ) ) )=Σr{1→l}( ( mr. - m )Σc{1→k}( ( m.c - m ) ) )
=Σr{1→l}( ( mr. - m )・k( m - m ) ) ) = 0

となるので、

SE=ST + SR + SC - 2( SR + SC - 0 )
=ST - SR - SC
ST=SR + SC + SE

という関係式が得られます。全体の変動 ST は、行間変動 SR、列間変動 SC、誤差変動 SE の和で表されることをこの式は示しています。

ST の期待値 E[ST] は、

E[ST]=E[Σr{1→l}( Σc{1→k}( ( xrc - m )2 ) )]
=E[Σr{1→l}( Σc{1→k}( xrc2 - 2mxrc + m2 ) )]
=Σr{1→l}( Σc{1→k}( E[xrc2] ) ) - 2E[mΣr{1→l}( Σc{1→k}( xrc ) )] + klm2
=Σr{1→l}( Σc{1→k}( E[xrc2] ) ) - 2E[m・klm] + klm2
=Σr{1→l}( Σc{1→k}( E[xrc2] ) ) - klE[m2]

となります。E[m2] は一元配置分散分析にて求めたときと同様な方法で

E[m2]=E[{ Σr{1→l}( kmr. ) }2] / ( kl )2
=E[k2Σr{1→l}( mr.2 ) + k2Σr{1→l}( Σs{1→l;s≠r}( mr.ms. ) )] / ( kl )2
={ Σr{1→l}( E[mr.2] ) + Σr{1→l}( Σs{1→l;s≠r}( E[mr.]E[ms.] ) ) } / l2
={ Σr{1→l}( σ2 / k + μr.2 ) + Σr{1→l}( Σs{1→l;s≠r}( μr.μs. ) ) } / l2
={ lσ2 / k + Σr{1→l}( μr.2 + Σs{1→l;s≠r}( μr.μs. ) ) } / l2
={ lσ2 / k + Σr{1→l}( μr.Σs{1→l}( μs. ) ) } / l2
={ lσ2 / k + Σr{1→l}( μr.・lμ ) } / l2
=( lσ2 / k + l2μ2 ) / l2
=σ2 / kl + μ2

と求められ、第一項目は

Σr{1→l}( Σc{1→k}( E[xrc2] ) )=Σr{1→l}( Σc{1→k}( σ2 + μrc2 ) )
=Σr{1→l}( Σc{1→k}( σ2 + μrc2 ) )
=klσ2 + Σr{1→l}( Σc{1→k}( μrc2 ) )

と求められるので、

E[ST]=klσ2 + Σr{1→l}( Σc{1→k}( μrc2 ) ) - kl( σ2 / kl + μ2 )
=( kl - 1 )σ2 + Σr{1→l}( Σc{1→k}( μrc2 - μ2 ) )
=( kl - 1 )σ2 + Σr{1→l}( Σc{1→k}( ( μrc - μ )2 ) )

という結果が得られます。μrc = μ + αr + βc を代入すると、

E[ST]=( kl - 1 )σ2 + Σr{1→l}( Σc{1→k}( { ( μ + αr + βc ) - μ }2 ) )
=( kl - 1 )σ2 + Σr{1→l}( Σc{1→k}( ( αr + βc )2 ) )
=( kl - 1 )σ2 + Σr{1→l}( Σc{1→k}( αr2 + 2αrβc + βc2 ) )
=( kl - 1 )σ2 + kΣr{1→l}( αr2 ) + 2Σr{1→l}( αrΣc{1→k}( βc ) ) + lΣc{1→k}( βc2 )
=( kl - 1 )σ2 + kΣr{1→l}( αr2 ) + lΣc{1→k}( βc2 )

と表すこともできます。SR の期待値 E[SR] は、

E[SR]=E[kΣr{1→l}( ( mr. - m )2 )]
=E[Σr{1→l}( kmr.2 - 2kmmr. + km2 )]
=r{1→l}( E[mr.2] ) - 2k・E[mΣr{1→l}( mr. )] + kl・E[m2]
=r{1→l}( E[mr.2] ) - 2k・E[lm2] + kl・E[m2]
=r{1→l}( E[mr.2] ) - kl・E[m2]

となり、E[mr.2] についても一元配置分散分析でのやり方と同様にして

E[mr.2]=E[{ Σc{1→k}( xrc ) }2] / k2
={ Σc{1→k}( E[xrc2] ) + 2Σc{1→k}( Σd{c+1→k}( E[xrcxrd] ) ) } / k2
={ Σc{1→k}( σ2 + μr.2 ) + 2Σc{1→k}( Σd{c+1→k}( E[xrc]E[xrd] ) ) } / k2
={ k( σ2 + μr.2 ) + 2Σc{1→k}( Σd{c+1→k}( μr.2 ) ) } / k2
={ k( σ2 + μr.2 ) + k( k - 1 )μr.2 } / k2
=( kσ2 + k2μr.2 ) / k2
=σ2 / k + μr.2

という結果が得られるので、

E[SR]=r{1→l}( σ2 / k + μr.2 ) - kl( σ2 / kl + μ2 )
=( l - 1 )σ2 + kΣr{1→l}( μr.2 ) - klμ2
=( l - 1 )σ2 + kΣr{1→l}( μr.2 - μ2 )
=( l - 1 )σ2 + kΣr{1→l}( ( μr. - μ )2 )

となります。これも μr. = μ + αr を代入すると、

E[SR]=( l - 1 )σ2 + kΣr{1→l}( { ( μ + αr ) - μ }2 )
=( l - 1 )σ2 + kΣr{1→l}( αr2 )

と表すことができます。同様にして SC の期待値 E[SC] は

E[SC]=( k - 1 )σ2 + lΣc{1→k}( ( μ.c - μ )2 )
=( k - 1 )σ2 + lΣc{1→k}( βc2 )

となるので、SE の期待値 E[SE] は

E[SE]=E[ST - SR - SC]
={ ( kl - 1 )σ2 + kΣr{1→l}( αr2 ) + lΣc{1→k}( βc2 ) }
 - { ( l - 1 )σ2 + kΣr{1→l}( αr2 ) } - { ( k - 1 )σ2 + lΣc{1→k}( βc2 ) }
=( kl - l - k + 1 )σ2
=( k - 1 )( l - 1 )σ2

になり、SE は αr, βc の値、すなわち行間・列間の変動に影響されない一定の値であることがこの結果からも示されたことになります。

ST の期待値は、行間・列間の両方の変動に影響し、変動が大きいほど値は大きくなる傾向にあります。また、SR の期待値は行間の、SC の期待値は列間の変動が大きくなるほど大きくなります。しかし、SE の期待値はこれらの変動に影響しないので、

FR = { SR / ( l - 1 ) } / { SE / ( k - 1 )( l - 1 ) } = ( k - 1 )SR / SE

FC = { SC / ( k - 1 ) } / { SE / ( k - 1 )( l - 1 ) } = ( l - 1 )SC / SE

とすれば、FR は行間の変動がない場合 ( つまり、μ1. = μ2. = ... = μr. = μ の場合 )、1 と等しくなり、変動が大きくなるほど値が大きくなります。また、FC は列間の変動がない場合 ( つまり、μ.1 = μ.2 = ... = μ.c = μ の場合 )、1 と等しくなり、変動が大きくなるほど値が大きくなります。

一元配置分散分析にて SC / σ2 = Σi{1→k}( ni( mi - m )2 ) / σ2 が自由度 k - 1 の χ2-分布に従う(補足1)ことを利用すれば、xrc が全て互いに独立であるという仮定のもとで SR / σ2 は ni = k となった特別な場合と考えられるので、もし行間変動がなければ自由度 l - 1 の χ2-分布に従い、同様に SC / σ2 は列間変動がないとき自由度 k - 1 の χ2-分布に従うことになります。

最後に、SE / σ2 が自由度 ( k - 1 )( l - 1 ) の χ2-分布に従うことを利用すると(補足2)、FR は自由度 ( l - 1, ( k - 1 )( l - 1 ) ) の、また FC は自由度 ( k - 1, ( k - 1 )( l - 1 ) ) の F-分布にそれぞれ従うので(*2-1)、以下のような分散分析表を利用した検定を行うことができます。これを「二元配置分散分析法(Two-way ANOVA)」といいます。

平方和自由度不偏分散F値
行間変動SR = kΣr{1→l}( ( mr. - m )2 )l - 1vR = SR / ( l - 1 )FR = { SR / ( l - 1 ) } / { SE / ( k - 1 )( l - 1 ) }
列間変動SC = lΣc{1→k}( ( m.c - m )2 )k - 1vC = SC / ( k - 1 )FC = { SC / ( k - 1 ) } / { SE / ( k - 1 )( l - 1 ) }
誤差変動SE = Σr{1→l}( Σc{1→k}( ( xrc - mr. - m.c + m )2 ) )( k - 1 )( l - 1 )vE = SE / ( k - 1 )( l - 1 ) 
全体ST = Σr{1→l}( Σc{1→k}( ( xrc - m )2 ) )kl - 1  

一元配置分散分析の場合と同様に、集団間の変動の度合いを示す寄与率を計算することができます。但し、今回は行間と列間の二つの変動が存在するので、次のような計算式を用います。

行間変動における寄与率 RR2 = SR / ( SR + SE )

列間変動における寄与率 RC2 = SC / ( SC + SE )

どちらの場合も、行間・列間の変動が小さいほど値は小さく、逆に大きいほど寄与率は大きくなります。極端な場合、行間・列間の変動がゼロであれば、値はゼロになり、誤差変動が無視できるほど行間・列間の変動が大きければほぼ 1 に等しくなります。

*2-1)(6) 標本分布」の「F-分布に対する性質」を参照


二元配置分散分析を行うためのサンプル・プログラムを以下に示します。

/*
  twoWayANOVA : 二元配置分散分析(two-way ANOVA)

  const vector< vector<double> >& data : データ列
  double a : 危険率
  double threshold : binSearchで棄却域を求める時のしきい値

  戻り値 : True ... 行間・列間変動のいずれかの仮説が棄却 , False ... どちらの仮説も保留 または データが条件を満たさない
*/
bool twoWayANOVA( const vector< vector<double> >& data, double a, double threshold )
{
  cout << "***** Two-way ANOVA *****" << endl << endl;

  unsigned int l = data.size(); // 行数
  if ( l < 2 ) {
    cerr << "At least two rows are necessary." << endl;
    return( false );
  }

  unsigned int k = data[0].size(); // 列数
  if ( k < 2 ) {
    cerr << "At least two columns are necessary." << endl;
    return( false );
  }
  vector<double> mr( l );    // 行毎の平均
  vector<double> mc( k, 0 ); // 列毎の平均
  double m = 0; // 全体の平均

  // 平均値の算出
  for ( unsigned int r = 0 ; r < l ; ++r ) {
    if ( data[r].size() != k ) {
      cerr << "The size of data[" << r << "] is not equal to that of data[0]." << endl;
      return( false );
    }
    mr[r] = sum( data[r] );
    m += mr[r];
    mr[r] /= (double)k;
    for ( unsigned int c = 0 ; c < k ; ++c )
      mc[c] += data[r][c];
  }
  for ( unsigned int c = 0 ; c < k ; ++c )
    mc[c] /= (double)l;
  m /= (double)( k * l );

  cout << "Row size (l) = " << l << " ; Column size (k) = " << k << endl;
  cout << "Mean. (m) = " << m << endl;
  cout << "Mean. by row (mr.) : ";
  for ( unsigned int r = 0 ; r < l - 1 ; ++r )
    cout << mr[r] << ", ";
  cout << mr[l - 1] << endl;
  cout << "Mean. by column (m.c) : ";
  for ( unsigned int c = 0 ; c < k - 1 ; ++c )
    cout << mc[c] << ", ";
  cout << mc[k - 1] << endl << endl;

  // 変動の算出
  Deviation<double> dev( m );
  double sr = sum( mr, dev ) * (double)k; // 行間変動
  double sc = sum( mc, dev ) * (double)l; // 列間変動
  double st = 0; // 全体の変動
  for ( unsigned int r = 0 ; r < l ; ++r )
    st += sum( data[r], dev );
  double se = st - sr - sc; // 誤差変動

  // F値の算出
  double fr = sr * (double)( k - 1 ) / se; // 行間変動に対する F値
  double fc = sc * (double)( l - 1 ) / se; // 列間変動に対する F値

  // 結果出力
  cout << "<Two-way ANOVA Table>" << endl;
  cout << "\t\tSS\tDOF\tUnbiased Var." << endl;
  cout << "Between-row\t" << sr << "\t" << l - 1 << "\t" << sr / (double)( l - 1 ) << "\t\tFr = " << fr << endl;
  cout << "Between-column\t" << sc << "\t" << k - 1 << "\t" << sc / (double)( k - 1 ) << "\t\tFc = " << fc << endl;
  cout << "Error\t\t" << se << "\t" << ( k - 1 ) * ( l - 1 ) << "\t" << se / (double)( ( k - 1 ) * ( l - 1 ) ) << endl;
  cout << "Total.\t\t" << st << "\t" << k * l - 1 << endl << endl;

  // 寄与率の出力
  cout << "Rr^2 = " << 100.0 * sr / ( sr + se ) << " %" << endl;
  cout << "Rc^2 = " << 100.0 * sc / ( sc + se ) << " %" << endl << endl;

  // 検定
  FDistribution fDistR( l - 1, ( k - 1 ) * ( l - 1 ) );
  double fr0 = binSearch( fDistR, 1.0 - a, threshold );
  cout << "Fr0 = " << fr0 << endl;
  double pr = 1.0 - fDistR.lower_p( fr );
  cout << "p-value = " << pr << endl << endl;

  FDistribution fDistC( k - 1, ( k - 1 ) * ( l - 1 ) );
  double fc0 = binSearch( fDistC, 1.0 - a, threshold );
  cout << "Fc0 = " << fc0 << endl;
  double pc = 1.0 - fDistC.lower_p( fc );
  cout << "p-value = " << pc << endl << endl;

  bool res = false; // 検定結果
  if ( fr0 < fr ) {
    cout << "Null hypothesis ( average of all samples by rows are same ) is rejected." << endl;
    res = true;
  } else {
    cout << "Null hypothesis ( average of all samples by rows are same ) is accepted." << endl;
  }

  if ( fc0 < fc ) {
    cout << "Null hypothesis ( average of all samples by columns are same ) is rejected." << endl;
    res = true;
  } else {
    cout << "Null hypothesis ( average of all samples by columns are same ) is accepted." << endl;
  }

  return( res );
}

一元配置分散分析の場合と異なり、今回は各行・各列のデータ数はそれぞれ等しい必要があるため、データのチェック方法は oneWayANOVA の場合とは異なっています。行間・列間・全体の平均をそれぞれ mr, mc, m として求め、さらにそれを利用して行間変動 sr、列間変動 sc、全体の変動 st をそれぞれ求めます。誤差変動 sest, sr, sc から求められるので、あとは行間・列間変動に対する F 値を求め、最後に検定を行います。


二元配置分散分析は、データが変動する要因が二つあって、各要因によってデータに有意な差があるかを検出するときに利用することになります。例えば、ある工具を使って製品を作るとき、工具と作業者の両方の要因による作業時間の変動を調べるような場合、二元配置分散分析が利用できる可能性があります。しかし、一元配置分散分析の場合と同様に、各データは等しい分散を持つことが前提になっているので、作業者が特定の工具を使って作業するときの作業時間のばらつきが個々の条件で異なるような場合は利用できません。
他の利用方法としては、各要因ごとの時系列の変化をデータとして得たとき、要因毎と時期のそれぞれに対して有意な差があるかを検出するような場合が考えられます。例えば、五人の作業者が新しい工具を使って作業を行い、その作業時間の平均を毎日一週間分得たとします。

作業日1234567平均
作業者
A65.360.859.859.450.850.850.656.8
B63.159.659.356.246.348.951.054.9
C65.360.157.252.953.949.648.455.3
D61.359.752.053.147.550.452.853.8
E62.556.154.052.747.748.152.753.4
平均63.559.356.554.949.249.651.154.9
作業時間の時系列推移

上図の左側がその結果を示した表で、右側は時系列推移をグラフで示したものです。グラフを見る限り、最初の頃に比較して後半は作業時間が短縮されていく様子が分かります。各作業者の時系列推移はかなりよく似ており、作業者の間にはほとんど差がないように見えます。実際にサンプル・プログラムで検定した結果は次のようになります。

***** Two-way ANOVA *****

Row size (l) = 5 ; Column size (k) = 7
Mean. (m) = 54.8543
Mean. by row (mr.) : 56.7857, 54.9143, 55.3429, 53.8286, 53.4
Mean. by column (m.c) : 63.5, 59.26, 56.46, 54.86, 49.24, 49.56, 51.1

<Two-way ANOVA Table>
                SS      DOF    Unbiased Var.
Between-row     49.9783 4      12.4946         Fr = 2.68696
Between-column  851.907 6      141.984         Fc = 30.5338
Error           111.602 24     4.65007
Total.          1013.49 34

Rr^2 = 30.931 %
Rc^2 = 88.4172 %

Fr0 = 2.77629
p-value = 0.0555365

Fc0 = 2.50819
p-value = 4.23296e-10

Null hypothesis ( average of all samples by rows are same ) is accepted.
Null hypothesis ( average of all samples by columns are same ) is rejected.

作業者ごとの作業時間の変動を表す行間変動に対する F 値 Fr は約 2.78 で、その p 値は 0.05 より大きいため、各作業者の作業時間平均が等しいという仮説は保留されるのに対し、時系列変動を表す列間変動に対する F 値 Fc は約 2.51 で、その p 値は 10-10 のオーダーと非常に小さく、時系列変動が変わらないという仮説は棄却されます。作業者毎の作業時間の変動に対する寄与率は 30% 程度であり、誤差変動の占める割合の方が高いという傾向が見られるのに対し、時系列変動に対する寄与率は 90% 程度あり、作業日毎の変動が占める割合は非常に高いことも上記結果から得ることができます。


3) 反復のある二元配置分散分析法(Two-way Repeated-Measures ANOVA)

二元配置分散分析法では、二つの要因に対する標本はただ一つのみでした。しかし、計測を繰り返して行うなどして、複数の標本が得られる場合も当然あり得ることであり、二元配置分散分析では平均などの代表値を利用した検定しかできないことになります。複数の標本に対して検定を行う手法として「反復のある二元配置分散分析法(Two-way Repeated-Measures ANOVA)」があります。これと区別するため、先に説明した二元配置分散分析法は「反復のない二元配置分散分析法(Two-way Factorial ANOVA)」と呼んで区別する場合もあります。

二つの要因 Ar ( r = 1, 2, ... l ) と Bc ( c = 1, 2, ... k ) に対し、その組み合わせ ( Ar, Bc ) ごとに n 個の標本 xrci ( i = 1, 2, ... n ) を割り当てるとします。各要因 Ar, Bc ごとの平均は

mr.. = Σc{1→k}( Σi{1→n}( xrci ) ) / kn

m.c. = Σr{1→l}( Σi{1→n}( xrci ) ) / ln

で表し、全体の平均を

m = Σr{1→l}( Σc{1→k}( Σi{1→n}( xrci ) ) ) / kln

とします。また、二つの要因の組み合わせ ( Ar, Bc ) ごとの標本の平均を

mrc. = Σi{1→n}( xrci ) / n

とすれば、(反復のない)二元配置分散分析法は mrc. を使った検定と考えることができます。よって、

ST'=Σr{1→l}( Σc{1→k}( ( mrc. - m )2 ) )
SR'=r{1→l}( ( mr.. - m )2 )
SC'=c{1→k}( ( m.c. - m )2 )
SE'=Σr{1→l}( Σc{1→k}( ( mrc. - mr.. - m.c. + m )2 ) )

とすれば、これは(反復のない)二元配置分散分析法における変動量と同じ式になり、

ST' = SR' + SC' + SE'

が成り立ちます。ここで、

ST' = Σr{1→l}( Σc{1→k}( mrc.2 - m2 ) )

であることを利用して、

ST'=Σr{1→l}( Σc{1→k}( { Σi{1→n}( xrci2 ) / n - m2 } - { Σi{1→n}( xrci2 ) / n - mrc.2 } ) )
=Σr{1→l}( Σc{1→k}( Σi{1→n}( xrci2 ) / n - m2 ) )
 - Σr{1→l}( Σc{1→k}( Σi{1→n}( xrci2 ) / n - mrc.2 ) )
=Σr{1→l}( Σc{1→k}( Σi{1→n}( xrci2 - m2 ) ) ) / n
 - Σr{1→l}( Σc{1→k}( Σi{1→n}( xrci2 - mrc.2 ) ) ) / n
=Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - m )2 ) ) ) / n
 - Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - mrc. )2 ) ) ) / n

のように分解すれば、

ST=Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - m )2 ) ) )
SR=knΣr{1→l}( ( mr.. - m )2 )
SC=lnΣc{1→k}( ( m.c. - m )2 )
SRC=r{1→l}( Σc{1→k}( ( mrc. - mr.. - m.c. + m )2 ) )
SE=Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - mrc. )2 ) ) )

と新たに平方和を定義することによって、

ST = SR + SC + SRC + SE

という関係式が得られます。SRC は、ある行・列の要素の全体との差から行間・列間変動を除外した量を表し、「交互作用」と呼ばれます。SE は誤差変動になり、この値だけが確率変数に従うものとすれば、そのモデルは次のように表されることになります。

xrci = μ + αr + βc + γrc + εrci

ここで、反復のない二元配置分散分析法の場合と同様に、xrci が正規分布 N( μrc., σ2 ) に従い互いに独立であるとしたとき、μ は全体の平均、αr は各行の平均 μr.. と μ との差、βc は各列の平均 μ.c. と μ との差を表すことになります。さらに、γrc は μrc. と μ, αr, βc との差になり、最後の εrci が正規分布 N( 0, σ2 ) に従う確率変数で表される誤差変動となります。これらを式で表すと、

αr = μr.. - μ = Σc{1→k}( μrc. ) / k - μ

βc = μ.c. - μ = Σr{1→l}( μrc. ) / l - μ

γrc = ( μrc. - μ ) - ( μr.. - μ ) - ( μ.c. - μ )

になり、Σr{1→l}( αr ) = Σc{1→k}( βc ) = 0 が成り立つと同時に、

Σr{1→l}( γrc )=Σr{1→l}( μrc. - μ - αr - βc )
=Σr{1→l}( ( μrc. - μ ) - ( μr.. - μ ) - ( μ.c. - μ ) )
=l( μ.c. - μ ) - l( μ - μ ) - l( μ.c. - μ ) = 0
Σr{1→k}( γrc )=Σc{1→k}( μrc. - μ - αr - βc )
=Σc{1→k}( ( μrc. - μ ) - ( μr.. - μ ) - ( μ.c. - μ ) )
=k( μr.. - μ ) - k( μr.. - μ ) - k( μ - μ ) = 0

が成り立つことになります。

各平方和の期待値を求める前に、各平均の平方 mrc.2, mr..2, m.c.2, m2 の期待値を求めておきます。

E[mrc.2]=E[{ Σi{1→n}( xrci ) / n ) }2]
=E[Σi{1→n}( xrci2 ) + 2Σi{1→n}( Σj{i+1→n}( xrcixrcj ) )] / n2
={ Σi{1→n}( E[xrci2] ) + 2Σi{1→n}( Σj{i+1→n}( E[xrci]E[xrcj] ) ) } / n2
={ Σi{1→n}( σ2 + μrc.2 ) + 2Σi{1→n}( Σj{i+1→n}( μrc.2 ) ) } / n2
={ n( σ2 + μrc.2 ) + n( n - 1 )μrc.2} / n2
=σ2 / n + μrc.2
E[mr..2]=E[{ Σc{1→k}( mrc. ) / k }2]
=E[Σc{1→k}( mrc.2 ) + Σc{1→k}( Σd{1→k;d≠c}( mrc.mrd. ) ) )] / k2
={ Σc{1→k}( E[mrc.2] ) + Σc{1→k}( Σd{1→k;d≠c}( E[mrc.]E[mrd.] ) ) } / k2
={ Σc{1→k}( σ2 / n + μrc.2 ) + Σc{1→k}( Σd{1→k;d≠c}( μrc.μrd. ) ) } / k2
=σ2 / kn + Σc{1→k}( μrc.2 + Σd{1→k;d≠c}( μrc.μrd. ) ) / k2
=σ2 / kn + Σc{1→k}( μrc.Σd{1→k}( μrd. ) ) / k2
=σ2 / kn + μr..Σc{1→k}( μrc. ) / k
=σ2 / kn + μr..2
E[m.c.2]=σ2 / ln + μ.c.2
E[m2]=E[{ Σr{1→l}( mr.. ) / l }2]
=E[Σr{1→l}( mr..2 ) + Σr{1→l}( Σs{1→l;s≠r}( mr..ms.. ) )] / l2
={ Σr{1→l}( E[mr..2] ) + Σr{1→l}( Σs{1→l;s≠r}( E[mr..]E[ms..] ) ) } / l2
={ Σr{1→l}( σ2 / kn + μr..2 ) + Σr{1→l}( Σs{1→l;s≠r}( μr..μs.. ) ) } / l2
=σ2 / kln + Σr{1→l}( μr..2 + Σs{1→l;s≠r}( μr..μs.. ) ) / l2
=σ2 / kln + Σr{1→l}( μr..Σs{1→l}( μs.. ) ) / l2
=σ2 / kln + μΣr{1→l}( μr.. ) / l
=σ2 / kln + μ2

解き方は反復のない二元配置分散分析法の時とほとんど変わりませんが、データ数が異なることから、二元配置分散分析法の時と値が多少異なります。次に、mr..mrc., m.c.mrc., mmr.., mm.c., mrc.xrci, mxrci の期待値を求めると

E[mr..mrc.]=E[{ Σd{1→k}( mrd. ) / k }mrc.]
=E[mrc.2 + Σd{1→k;d≠c}( mrc.mrd. )] / k
={ σ2 / n + μrc.2 + Σd{1→k;d≠c}( μrc.μrd. ) } / k
=σ2 / kn + μrc.Σd{1→k}( μrd. ) / k
=σ2 / kn + μr..μrc.
E[m.c.mrc.]=σ2 / ln + μ.c.μrc.
E[mmr..]=E[{ Σs{1→l}( ms.. ) / l }mr..]
=E[mr..2 + Σs{1→l;s≠r}( mr..ms.. )] / l
={ σ2 / kn + μr..2 + Σs{1→l;s≠r}( μr..μs.. ) } / l
=σ2 / kln + μr..Σs{1→l}( μs.. ) / l
=σ2 / kln + μμr..
E[mm.c.]=σ2 / kln + μμ.c.
E[mrc.xrci]=E[xrciΣj{1→n}( xrcj ) / n }]
=E[xrci2 + Σj{1→n;j≠i}( xrcixrcj )] / n
={ ( σ2 + μrc.2 ) + Σj{1→n;j≠i}( μrc.2 ) } / n
=σ2 / n + μrc.2
E[mxrci]=E[xrciΣs{1→l}( Σd{1→k}( msd. ) ) / kl]
=E[mrc.xrci + Σs{1→l}( Σd{1→k;(s,d)≠(r,c)}( msd.xrci ) )] / kl
={ ( σ2 / n + μrc.2 ) + Σs{1→l}( Σd{1→k;(s,d)≠(r,c)}( μsd.μrc. ) ) } / kl
={ σ2 / n + μrc.Σs{1→l}( Σd{1→k}( μsd. ) ) } / kl
=( σ2 / n + μrc.・klμ ) / kl
=σ2 / kln + μμrc.

になります。最後の E[mxrci] の算出において、Σd{1→k;(s,d)≠(r,c)}(...) という表現は s = r かつ d = c の場合だけを和の項に含めないことを意味しています。各要素が互いに独立であることから、二乗和にならない限り E[XY] = E[X]E[Y] として計算することができますが、二乗和の場合は計算方法が異なるためその項だけ特別扱いする必要があります。上記のどの計算においても、基本的には、二乗和を抽出してまた和の中に戻すという操作を行っているだけです。

以上の結果を使い、E[ST], E[SR], E[SC], E[SE] を求めると

E[ST]=E[Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - m )2 ) ) )]
=Σr{1→l}( Σc{1→k}( Σi{1→n}( E[xrci2 - 2mxrci + m2] ) ) )
=Σr{1→l}( Σc{1→k}( Σi{1→n}( ( σ2 + μrc.2 )
 - 2( σ2 / kln + μμrc. ) + ( σ2 / kln + μ2 ) ) ) )
=( kln - 1 )σ2 + nΣr{1→l}( Σc{1→k}( μrc.2 - 2μμrc. + μ2 ) )
=( kln - 1 )σ2 + nΣr{1→l}( Σc{1→k}( ( μrc. - μ )2 ) )
E[SR]=E[knΣr{1→l}( ( mr.. - m )2 )]
=knΣr{1→l}( E[mr..2 - 2mmr.. + m2] )
=knΣr{1→l}( ( σ2 / kn + μr..2 )
 - 2( σ2 / kln + μμr.. ) + ( σ2 / kln + μ2 ) )
=2 - 2σ2 + σ2 + knΣr{1→l}( μr..2 - 2μμr.. + μ2 )
=( l - 1 )σ2 + knΣr{1→l}( ( μr.. - μ )2 )
=( l - 1 )σ2 + knΣr{1→l}( αr2 )
E[SC]=( k - 1 )σ2 + lnΣc{1→k}( ( μ.c. - μ )2 )
=( k - 1 )σ2 + lnΣc{1→k}( βc2 )
E[SE]=E[Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - mrc. )2 ) ) )]
=Σr{1→l}( Σc{1→k}( Σi{1→n}( E[xrci2 - 2mrc.xrci + mrc.2] ) ) )
=Σr{1→l}( Σc{1→k}( Σi{1→n}( ( σ2 + μrc.2 )
 - 2( σ2 / n + μrc.2 ) + ( σ2 / n + μrc.2 ) ) ) )
=kln( σ2 - 2σ2 / n + σ2 / n )
=kl( n - 1 )σ2

最後に、S[SRC] は

E[SRC]=E[ST - SR - SC - SE]
={ ( kln - 1 )σ2 + nΣr{1→l}( Σc{1→k}( ( μrc. - μ )2 ) ) }
 - { ( l - 1 )σ2 + knΣr{1→l}( αr2 ) } - { ( k - 1 )σ2 + lnΣc{1→k}( βc2 ) } - kl( n - 1 )σ2
={ ( kln - 1 ) - ( l - 1 ) - ( k - 1 ) - kl( n - 1 ) }σ2
 + nΣr{1→l}( Σc{1→k}( ( μrc. - μ )2 ) ) - knΣr{1→l}( αr2 ) - lnΣc{1→k}( βc2 )
=( k - 1 )( l - 1 )σ2 + nΣr{1→l}( Σc{1→k}( ( μrc. - μ )2 - αr2 - βc2 ) )

となって、全ての平方和に対する期待値が得られました。以上の結果から、全体の変動 ST は行・列の要素に対する平均 μrc. が全て等しいとき最小値 ( kln - 1 )σ2 となり、各要素の平均のズレが大きくなるほど変動も大きくなることがわかります。同様に、行・列ごとの変動 SR, SC はそれぞれ、行・列ごとの平均 μr.., μ.c. がすべて等しい時に最小値 ( l - 1 )σ2, ( k - 1 )σ2 となって、ズレが大きくに従って変動も大きくなります。また、交互作用の変動 SRC は全体の変動 ST と同様に、μrc. が全て等しければ(行間・列間の変動もなくなるため)最小値 ( k - 1 )( l - 1 )σ2 となり、各要素の平均のズレが大きくなるほど変動も大きくなりますが、行・列間のズレが大きくなると逆に変動は小さくなり、逆に小さければ変動は大きくなる傾向にあります。各行・列の平均には差があまりないのに、要素間の変動は大きいというようなときに SRC は大きくなるため「交互作用」という名が付けられています。

SE は行・列間及び各要素の変動に関係なく一定の値を取ることが期待できることから、

FR = { SR / ( l - 1 ) } / { SE / kl( n - 1 ) }

FC = { SC / ( k - 1 ) } / { SE / kl( n - 1 ) }

FRC = { SRC / ( k - 1 )( l - 1 ) } / { SE / kl( n - 1 ) }

としたとき、変動がなければ各 F 値は 1 となり、変動が大きくなるほどその値は大きくなります。一元配置分散分析にて、mi の平均が全て等しいならば、

SC / σ2 = Σi{1→k}( ni( mi - m )2 ) / σ2

が自由度 k - 1 の χ2-分布に従う(補足1)ことを利用すれば、SR / σ2 と SC / σ2 はそれぞれ ni = kn, ln となった特別な場合と考えられるので、自由度 l - 1, k - 1 の χ2-分布に従うことになります。
これは次のように示すこともできます。mr.. と m.c. はそれぞれ kn 個、ln 個の要素の標本平均であり、各要素が同じ正規分布 N( μ, σ2 ) に従うならば、mr.. と m.c. はそれぞれ正規分布 N( μ, σ2 / kn ), N( μ σ2 / ln ) に従うことになるので (*3-1)、

Σr{1→l}( ( mr.. - m )2 ) / ( σ2 / kn )

Σc{1→k}( ( m.c. - m )2 ) / ( σ2 / ln )

は自由度 l - 1, k - 1 の χ2-分布に従うことになります。

また、反復のない二元配置分散分析にて

Σr{1→l}( Σc{1→k}( ( xrc - mr. - m.c + m )2 ) ) / σ2

が自由度 ( k - 1 )( l - 1 ) の χ2-分布に従うことを示していますが(補足2)、SRC は要素が mrc. となっただけで、これは正規分布 N( μrc, σ2 / n ) に従うことになるので(*3-1)、前要素に対して μrc = μ ならば

Σr{1→l}( Σc{1→k}( ( mrc. - mr.. - m.c. + m )2 ) ) / ( σ2 / n ) = SRC / σ2

はやはり自由度 ( k - 1 )( l - 1 ) の χ2-分布に従うことになります。

最後に、SE については

Σi{1→n}( ( xrci - mrc. )2 )

が自由度 n - 1 の χ2-分布に従い、SE がその kl 個の和で表されることから、χ2-分布の再生性より自由度 kl( n - 1 ) の χ2-分布に従うという結果が得られます(*3-2)。従って、FR は自由度 ( l - 1, kl( n -1 ) ) の、FC は自由度 ( k - 1, kl( n - 1 ) ) の、そして FRC は自由度 ( ( k - 1 )( l - 1 ), kl( n - 1 ) ) の F-分布にそれぞれ従います(*3-3)。以上をまとめ、分散分析表として表すと

平方和自由度不偏分散F値
行間変動SR = knΣr{1→l}( ( mr.. - m )2 )l - 1vR = SR / ( l - 1 )FR = { SR / ( l - 1 ) } / { SE / kl( n - 1 ) }
列間変動SC = lnΣc{1→k}( ( m.c. - m )2 )k - 1vC = SC / ( k - 1 )FC = { SC / ( k - 1 ) } / { SE / kl( n - 1 ) }
交互作用SRC = nΣr{1→l}( Σc{1→k}( ( mrc. - mr.. - m.c. + m )2 ) )( k - 1 )( l - 1 )vRC = SRC / ( k - 1 )( l - 1 )FRC = { SRC / ( k - 1 )( l - 1 ) } / { SE / kl( n - 1 ) }
誤差変動SE = Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - mrc. )2 ) ) )kl( n - 1 )vE = SE / kl( n - 1 ) 
全体ST = Σr{1→l}( Σc{1→k}( Σi{1→n}( ( xrci - m )2 ) ) )kln - 1  

になります。

*3-1)(5) 正規分布」の「3) 大数の法則(Law of Large Numbers)と中心極限定理(Central Limit Theorem)」を参照

*3-2)(6) 標本分布」の「カイ二乗分布に対する性質」を参照

*3-3)(6) 標本分布」の「F-分布に対する性質」を参照


反復のある二元配置分散分析を行うためのサンプル・プログラムを以下に示します。

/*
  twoWay_RM_ANOVA : 反復のある二元配置分散分析(two-way Repeated-Measures ANOVA)

  const vector< vector< vector<double> > >& data : データ列
  double a : 危険率
  double threshold : binSearchで棄却域を求める時のしきい値

  戻り値 : True ... 行間・列間変動・交互作用のいずれかの仮説が棄却 , False ... どの仮説も保留 または データが条件を満たさない
*/
bool twoWay_RM_ANOVA( const vector< vector< vector<double> > >& data, double a, double threshold )
{
  cout << "***** Two-way Repeated-Measures ANOVA *****" << endl << endl;

  unsigned int l = data.size(); // 行数
  if ( l < 2 ) {
    cerr << "At least two rows are necessary." << endl;
    return( false );
  }

  unsigned int k = data[0].size(); // 列数
  if ( k < 2 ) {
    cerr << "At least two columns are necessary." << endl;
    return( false );
  }

  unsigned int n = data[0][0].size(); // 要素数
  if ( n < 2 ) {
    cerr << "At least two data are necessary." << endl;
    return( false );
  }

  vector<double> mr( l, 0 );   // 行毎の平均
  vector<double> mc( k, 0 );   // 列毎の平均
  vector<double> mrc( l * k ); // 要素(行列)毎の平均
  double m = 0; // 全体の平均

  // 平均値の算出
  for ( unsigned int r = 0 ; r < l ; ++r ) {
    if ( data[r].size() != k ) {
      cerr << "The size of data[" << r << "] is not equal to that of data[0]." << endl;
      return( false );
    }
    for ( unsigned int c = 0 ; c < k ; ++c ) {
      if ( data[r][c].size() != n ) {
        cerr << "The size of data[" << r << "][" << c << "] is not equal to that of data[0][0]." << endl;
        return( false );
      }
      mrc[r * k + c] = sum( data[r][c] );
      mr[r] += mrc[r * k + c];
      mc[c] += mrc[r * k + c];
      mrc[r * k + c] /= (double)n;
    }
    m += mr[r];
    mr[r] /= (double)( k * n );
  }
  for ( unsigned int c = 0 ; c < k ; ++c )
    mc[c] /= (double)( l * n );
  m /= (double)( k * l * n );

  cout << "Row size (l) = " << l << " ; Column size (k) = " << k << " ; Data size (n) = " << n << endl;
  cout << "Mean. (m) = " << m << endl;
  cout << "Mean. by row (mr..) : ";
  for ( unsigned int r = 0 ; r < l - 1 ; ++r )
    cout << mr[r] << ", ";
  cout << mr[l - 1] << endl;
  cout << "Mean. by column (m.c.) : ";
  for ( unsigned int c = 0 ; c < k - 1 ; ++c )
    cout << mc[c] << ", ";
  cout << mc[k - 1] << endl;
  cout << "Mean. by data (mrc.) : " << endl;
  for ( unsigned int r = 0 ; r < l ; ++r ) {
    for ( unsigned int c = 0 ; c < k - 1 ; ++c )
      cout << mrc[r * k + c] << ", ";
    cout << mrc[( r + 1 ) * k - 1] << endl;
  }
  cout << endl;

  // 行間・列間変動の算出
  Deviation<double> dev( m );
  double sr = sum( mr, dev ) * (double)( k * n ); // 行間変動
  double sc = sum( mc, dev ) * (double)( l * n ); // 列間変動

  // 全変動・誤差変動の算出
  double st = 0; // 全変動
  double se = 0; // 誤差変動
  for ( unsigned int r = 0 ; r < l ; ++r ) {
    for ( unsigned int c = 0 ; c < k ; ++c ) {
      Deviation<double> dev_rc( mrc[r * k + c] );
      st += sum( data[r][c], dev );
      se += sum( data[r][c], dev_rc );
    }
  }

  double src = st - sr - sc - se; // 交互作用

  // F値の算出
  double fr = sr * (double)( k * l * ( n - 1 ) ) / ( se * (double)( l - 1 ) ); // 行間変動に対する F値
  double fc = sc * (double)( k * l * ( n - 1 ) ) / ( se * (double)( k - 1 ) ); // 列間変動に対する F値
  double frc = src * (double)( k * l * ( n - 1 ) ) / ( se * (double)( ( k - 1 ) * ( l - 1 ) ) ); // 交互作用に対する F値

  // 結果出力
  cout << "<Two-way Repeated-Measures ANOVA Table>" << endl;
  cout << "\t\tSS\tDOF\tUnbiased Var." << endl;
  cout << "Between-row\t" << sr << "\t" << l - 1 << "\t" << sr / (double)( l - 1 ) << "\t\tFr = " << fr << endl;
  cout << "Between-column\t" << sc << "\t" << k - 1 << "\t" << sc / (double)( k - 1 ) << "\t\tFc = " << fc << endl;
  cout << "Interaction\t" << src << "\t" << ( k - 1 ) * ( l - 1 ) << "\t" << src / (double)( ( k - 1 ) * ( l - 1 ) ) << "\t\tFrc = " << frc << endl;
  cout << "Error\t\t" << se << "\t" << k * l * ( n - 1 ) << "\t" << se / (double)( k * l * ( n - 1 ) ) << endl;
  cout << "Total.\t\t" << st << "\t" << k * l * n - 1 << endl << endl;

  // 寄与率の出力
  cout << "Rr^2 = " << 100.0 * sr / ( sr + se ) << " %" << endl;
  cout << "Rc^2 = " << 100.0 * sc / ( sc + se ) << " %" << endl;
  cout << "Rrc^2 = " << 100.0 * src / ( src + se ) << " %" << endl << endl;

  // 検定
  FDistribution fDistR( l - 1, k * l * ( n - 1 ) );
  double fr0 = binSearch( fDistR, 1.0 - a, threshold );
  cout << "Fr0 = " << fr0 << endl;
  double pr = 1.0 - fDistR.lower_p( fr );
  cout << "p-value = " << pr << endl << endl;

  FDistribution fDistC( k - 1, k * l * ( n - 1 ) );
  double fc0 = binSearch( fDistC, 1.0 - a, threshold );
  cout << "Fc0 = " << fc0 << endl;
  double pc = 1.0 - fDistC.lower_p( fc );
  cout << "p-value = " << pc << endl << endl;

  FDistribution fDistRC( ( k - 1 ) * ( l - 1 ), k * l * ( n - 1 ) );
  double frc0 = binSearch( fDistRC, 1.0 - a, threshold );
  cout << "Frc0 = " << frc0 << endl;
  double prc = 1.0 - fDistRC.lower_p( frc );
  cout << "p-value = " << prc << endl << endl;

  bool res = false; // 検定結果
  if ( fr0 < fr ) {
    cout << "Null hypothesis ( average of all samples by rows are same ) is rejected." << endl;
    res = true;
  } else {
    cout << "Null hypothesis ( average of all samples by rows are same ) is accepted." << endl;
  }

  if ( fc0 < fc ) {
    cout << "Null hypothesis ( average of all samples by columns are same ) is rejected." << endl;
    res = true;
  } else {
    cout << "Null hypothesis ( average of all samples by columns are same ) is accepted." << endl;
  }

  if ( frc0 < frc ) {
    cout << "Null hypothesis ( No interaction between every sample ) is rejected." << endl;
    res = true;
  } else {
    cout << "Null hypothesis ( No interaction between every sample ) is accepted." << endl;
  }

  return( res );
}

かなり長いプログラムとなっていますが、一元・二元配置分散分析法の場合と処理の流れはほとんど変わりません。最初に各平均を求め、それを元に変動・F値を求めたら、最後に検定を行うだけです。寄与率は、行間変動・列間変動・交互作用それぞれについて求められ、各変動と誤差変動との和に対する各変動の比率で得られます。


反復のない二元配置分散分析法で示した例について、実際には各要素が三回繰り返して測定されていたとします。その結果を以下に示します(先の例では一回目の測定結果だけを使っていました)。



作業日1234567平均
作業者
A65.360.859.859.450.850.850.656.8
B63.159.659.356.246.348.951.054.9
C65.360.157.252.953.949.648.455.3
D61.359.752.053.147.550.452.853.8
E62.556.154.052.747.748.152.753.4
平均63.559.356.554.949.249.651.154.9


A65.061.161.761.059.051.250.958.6
B64.764.664.260.256.355.756.660.3
C67.264.759.757.556.257.055.059.6
D62.460.357.555.752.451.152.356.0
E60.761.459.359.358.049.549.056.8
平均63.661.059.158.655.050.550.056.8


A65.961.959.260.958.451.351.958.5
B65.464.763.761.154.147.047.357.6
C66.261.659.256.750.645.647.455.3
D64.562.960.861.961.057.856.560.8
E59.460.857.757.554.849.645.055.0
平均64.361.460.059.454.849.149.656.9
作業時間の時系列推移(反復測定)

上図の左側がその結果を示した表で、計三回の計測を実施しています。右側は時系列推移をグラフで示したもので、各測定回数ごとに色分けをしていますす。一回目だけではなく二・三回目の測定でも、最初の頃に比較して後半は作業時間が短縮しています。また、各作業者間の作業時間にはほとんど差がないように見えます。実際にサンプル・プログラムで検定した結果は次のようになります。

***** Two-way Repeated-Measures ANOVA *****

Row size (l) = 5 ; Column size (k) = 7 ; Data size (n) = 3
Mean. (m) = 56.8438
Mean. by row (mr..) : 57.9476, 57.619, 56.7619, 56.8524, 55.0381
Mean. by column (m.c.) : 63.9267, 61.3533, 59.02, 57.74, 53.8, 50.9067, 51.16
Mean. by data (mrc.) : 
65.4, 61.2667, 60.2333, 60.4333, 56.0667, 51.1, 51.1333
64.4, 62.9667, 62.4, 59.1667, 52.2333, 50.5333, 51.6333
66.2333, 62.1333, 58.7, 55.7, 53.5667, 50.7333, 50.2667
62.7333, 60.9667, 56.7667, 56.9, 53.6333, 53.1, 53.8667
60.8667, 59.4333, 57, 56.5, 53.5, 49.0667, 48.9

<Two-way Repeated-Measures ANOVA Table>
                SS      DOF    Unbiased Var.
Between-row     106.822 4      26.7056         Fr = 2.46545
Between-column  2292.93 6      382.154         Fc = 35.2804
Interaction     172.276 24     7.17818         Frc = 0.662689
Error           758.233 70     10.8319
Total.          3330.26 104

Rr^2 = 12.3486 %
Rc^2 = 75.1493 %
Rrc^2 = 18.5142 %

Fr0 = 2.50266
p-value = 0.0528002

Fc0 = 2.23119
p-value = 2.71051e-19

Frc0 = 1.67383
p-value = 0.869777

Null hypothesis ( average of all samples by rows are same ) is accepted.
Null hypothesis ( average of all samples by columns are same ) is rejected.
Null hypothesis ( No interaction between every sample ) is accepted.

列間変動すなわち時系列での変動には差があり、行間変動すなわち個人差と交互作用に対しては差がないという予想通りの結果が得られています。

交互作用については、次のような例を考えると分かりやすいと思います。

作業日1234567平均
作業者
A65.358.859.855.449.849.447.755.2
B50.151.551.152.258.562.561.655.4
C62.360.157.256.953.249.649.255.5
D48.853.452.556.158.960.158.155.4
E49.752.655.357.157.655.360.455.4
平均55.255.355.255.555.655.455.455.4
作業時間の時系列推移(交互作用のある例)

上に示したデータを見ると、行間・列間の平均はほとんど差がありません。しかし、グラフを見ると、作業者によって作業時間が長くなる場合と短くなる場合があり、それらによって変動が相殺されたような形になっています。これを反復のない二元配置分散分析法で検定すると次のような結果になります。

***** Two-way ANOVA *****

Row size (l) = 5 ; Column size (k) = 7
Mean. (m) = 55.3743
Mean. by row (mr.) : 55.1714, 55.3571, 55.5, 55.4143, 55.4286
Mean. by column (m.c) : 55.24, 55.28, 55.18, 55.54, 55.6, 55.38, 55.4

<Two-way ANOVA Table>
                SS            DOF     Unbiased Var.
Between-row     0.432571      4       0.108143               Fr = 0.00346901
Between-column  0.718857      6       0.11981                Fc = 0.00384325
Error           748.175       24      31.174
Total.          749.327       34

Rr^2 = 0.0577834 %
Rc^2 = 0.0959891 %

Fr0 = 2.77629
p-value = 0.999974

Fc0 = 2.50819
p-value = 1

Null hypothesis ( average of all samples by rows are same ) is accepted.
Null hypothesis ( average of all samples by columns are same ) is accepted.

行間・列間の変動はほとんどなく、全変動の大半は誤差変動に含まれています。しかし、グラフを見る限りはこれを測定誤差と見ることは誤っているように思えます。実はこの中に、交互作用による変動が混在している可能性があります。



作業日1234567平均
作業者
A65.358.859.855.449.849.447.755.2
B50.151.551.152.258.562.561.655.4
C62.360.157.256.953.249.649.255.5
D48.853.452.556.158.960.158.155.4
E49.752.655.357.157.655.360.455.4
平均55.255.355.255.555.655.455.455.4


A63.756.763.857.450.952.145.755.8
B51.649.749.553.456.164.963.655.5
C60.458.059.255.055.647.850.755.2
D49.652.853.056.257.160.758.755.4
E49.553.152.854.859.553.362.655.1
平均55.054.155.755.455.855.856.355.4


A63.660.758.956.946.854.447.255.5
B51.450.052.156.157.160.061.955.5
C63.963.458.154.053.546.448.455.4
D50.353.954.353.459.260.657.755.6
E50.156.351.955.257.355.461.355.4
平均55.956.855.055.154.855.355.355.5
作業時間の時系列推移(交互作用のある例の反復測定)

同じような測定を三回繰り返した結果が上記のようになったとします(一回目の測定結果は先の例と同一です)。一回目の測定で作業時間が短くなる傾向にあるのは作業者 A, C で、残り B, D, E は逆に長くなる傾向にありました。そこで、それぞれのグループで色分けをして三回の測定結果全てをプロットしたグラフが上記右側になります。グラフから見ても明らかに、時系列傾向が A, C と B, D, E で異なっていることが分かります。これを反復のある二元配置分散分析法で検定した結果は以下のようになります。

***** Two-way Repeated-Measures ANOVA *****

Row size (l) = 5 ; Column size (k) = 7 ; Data size (n) = 3
Mean. (m) = 55.4219
Mean. by row (mr..) : 55.4762, 55.4714, 55.3762, 55.4952, 55.2905
Mean. by column (m.c.) : 55.3533, 55.4, 55.3, 55.34, 55.4067, 55.5, 55.6533
Mean. by data (mrc.) : 
64.2, 58.7333, 60.8333, 56.5667, 49.1667, 51.9667, 46.8667
51.0333, 50.4, 50.9, 53.9, 57.2333, 62.4667, 62.3667
62.2, 60.5, 58.1667, 55.3, 54.1, 47.9333, 49.4333
49.5667, 53.3667, 53.2667, 55.2333, 58.4, 60.4667, 58.1667
49.7667, 54, 53.3333, 55.7, 58.1333, 54.6667, 61.4333

<Two-way Repeated-Measures ANOVA Table>
               SS             DOF           Unbiased Var.
Between-row    0.632952       4             0.158238       Fr = 0.0704981
Between-column 1.29962        6             0.216603       Fc = 0.0965009
Interaction    2253.37        24            93.8903        Frc = 41.8299
Error          157.12         70            2.24457
Total.         2412.42        104

Rr^2 = 0.40123 %
Rc^2 = 0.820365 %
Rrc^2 = 93.4818 %

Fr0 = 2.50266
p-value = 0.990736

Fc0 = 2.23119
p-value = 0.996517

Frc0 = 1.67383
p-value = 0

Null hypothesis ( average of all samples by rows are same ) is accepted.
Null hypothesis ( average of all samples by columns are same ) is accepted.
Null hypothesis ( No interaction between every sample ) is rejected.

交互作用に対する平方和は 2.25 x 103 で全変動のほとんどを占めています。また、検定結果を見ても、差がないという仮説は棄却されています。

交互作用の計算式は

SRC=r{1→l}( Σc{1→k}( ( mrc. - mr.. - m.c. + m )2 ) )
=r{1→l}( Σc{1→k}( { ( mrc. - m ) - ( mr.. - m ) - ( m.c. - m ) }2 ) )

であり、mr.. - m, m.c. - m が小さく、mrc. - m が大きいほど値が増加する傾向にあります。つまり、各要素の全体の平均とのばらつきがあったとして、行・列ごとのそれぞれの変動が似ていれば行間・列間変動が大きくなり、逆に変動がバラバラの場合は交互作用が大きくなります。


4) クラスカル・ワリス検定(Kruskal-Wallis Test)

ANOVA では、標本として測定値などの実際のデータがあることを前提としていました。しかし、テストの順位や成績(グレード)など、いわゆる順序尺度に対しては ANOVA をそのまま適用することができません。順序尺度に対して(一元配置)分散分析法に相当するものとして「クラスカル・ワリス検定(Kruskal-Wallis Test)」があります。

一元配置分散分析法で示した例に対し、キャベツの質量に関する情報が順位でしか得られなかったと仮定します。このとき、その順位は次のようになります(重いキャベツほど順位は高く(値が小さく)なります)。

肥料なし 肥料A  肥料B  全体 
実測値順位実測値順位実測値順位実測値順位
114911015115.5144911
21365151496915124
315115.51569315027
414361216321138914
514291314978  
6  15922  
合計723255.5929728.558523622381120
平均144611.115504.751463914928

値が等しい場合、適当に順位付けをした後でその順位の平均を求め、その結果をそれぞれに割り当てます。上記の例の場合は、実測値が 1511(g) となる要素が二つあるので、それぞれの順位 5, 6 を使い、その平均 ( 5 + 6 ) / 2 = 5.5 を割り当てています。

もし、それぞれの集団の平均と分散が等しければ、無作為に抽出した要素はランダムに並ぶことが期待できて、その結果としてそれぞれの順位の平均は等しくなることが予想できます。逆に分布が異なれば、ある集団だけ順位の平均が高くなったり低くなることも理解できると思います。上記の例では、肥料 A を使ったキャベツに対する順位の平均が他の二集団に比べて小さく、重いキャベツがこの集団に多く分布しているということが分かります。この差を使った検定方法が「クラスカル・ワリス検定(Kruskal-Wallis Test)」です。

k 個の集団 Ai ( i = 1, 2, ... k ) から ni 個の要素を抽出したとき、同じ母集団から抽出された要素は同一であると見なして並べる場合の順列は次式で表されます。

P( N ; n ) = P( N ; n1, n2, ... nk ) = N! / n1!・n2!・...・nk!

但し n1 + n2 + ... + nk = N

各要素は全てどの順位についても等しい確率でなりうると仮定すると、それぞれの要素の順位の平均 μ は 1 から N までの和(順位和)を N で割った値になるので

μ = N( N + 1 ) / 2N = ( N + 1 ) / 2

と求められます。この値は、1 から N までの数列の中央値にあたる数を表しています。また、「全体の変動」 ST

ST=Σi{1→N}( ( i - μ )2 )
=Σi{1→N}( i2 - 2iμ + μ2 )
=Σi{1→N}( i2 ) - 2μΣi{1→N}( i ) + Nμ2
=Σi{1→N}( i2 ) - Nμ2
=( 2N3 + 3N2 + N ) / 6 - N( N + 1 )2 / 4
=N( N2 - 1 ) / 12

となって、μ と ST は各集団の個数に依存しない値であることが分かります。

例として、k = 3, n1 = 2, n2 = 1, n3 = 1 の場合を考えると、順列は

P( 4 ; 2, 1, 1 ) = 4! / 2!・1!・1! = 12

となり、Ai の要素を i で表せば、並べ方は実際

( 1, 1, 2, 3 ), ( 1, 1, 3, 2 ), ( 1, 2, 1, 3 ), ( 1, 3, 1, 2 ),
( 1, 2, 3, 1 ), ( 1, 3, 2, 1 ), ( 2, 1, 1, 3 ), ( 3, 1, 1, 2 ),
( 2, 1, 3, 1 ), ( 3, 1, 2, 1 ), ( 2, 3, 1, 1 ), ( 3, 2, 1, 1 )

の 12 通りになります。どの要素も順位の平均は μ = ( 4 + 1 ) / 2 = 2.5 であり、全体の変動は ST = 4( 42 - 1 ) / 12 = 5 になります。

この例の中で、A1 の順位に着目すると、

( 1, 2 ), ( 1, 2 ), ( 1, 3 ), ( 1, 3 ),
( 1, 4 ), ( 1, 4 ), ( 2, 3 ), ( 2, 3 ),
( 2, 4 ), ( 2, 4 ), ( 3, 4 ), ( 3, 4 )

となって、1 から 4 までの数字がそれぞれ 6 回ずつ出現していることが分かります。並べ方が 12 通りでそれぞれに二つの要素があるので、要素数は 24 個になり、順位は 1 から 4 まであって、それぞれは等しい回数出現しているので 24 / 4 = 6 という値が得られることになります。従って、順位の和の合計は

( 1 + 2 + 3 + 4 ) x 6 = 60

で、その平均は 60 / 12 = 5 になります。この結果から、全体が N 個、Ai が ni 個の要素であるならば、順位の和 Ti の合計は

Σi{1→N}( i ) x { ni・P( N ; n ) / N } = ni( N + 1 )P( N ; n ) / 2

であり、その期待値 E[Ti] は P( N ; n ) で割った値 ni( N + 1 ) / 2 となることが分かります。これは、全体の平均 μ に要素数 ni を掛けた値と等しくなります。

以前、「(10) 順序統計量」の中で「2) U-分布」を紹介しました。順位和と U 値の間には次のような関係式が成り立ちます(*4-1)。

U = niui + ui( ui + 1 ) / 2 - Tx

但し、ni, ui はそれぞれ、ある母集団 Ai の要素数とそれ以外の要素数を、TxAi 以外の順位和をそれぞれ表します。このとき、Tx + Ti = N( N + 1 ) / 2 が成り立つので、上式は

U=ni( N - ni ) + ( N - ni )( N - ni + 1 ) / 2 - { N( N + 1 ) / 2 - Ti }
=Ti + ( N - ni )( 2ni + N - ni + 1 ) / 2 - N( N + 1 ) / 2
=Ti + [ ( N - ni ){ ( N + 1 ) + ni } - N( N + 1 ) ] / 2
=Ti + { N( N + 1 ) + Nni - ( N + 1 )ni - ni2 - N( N + 1 ) } / 2
=Ti - ni( ni + 1 ) / 2

になります。結局、ある母集団における順位和の確率分布は、その母集団の要素数 ni とそれ以外の母集団の要素数 ui = N - ni における U-分布 Pni,ui(U) に対して変数変換 Ti = U + ni( ni + 1 ) / 2 を行った(つまり分布を単純に ni( ni + 1 ) / 2 だけ正方向へシフトした)ものと等しくなります。U-分布 Pni,ui(U) の期待値は niui / 2 = ni( N - ni ) / 2 であり、この値に ni( ni + 1 ) / 2 を加えれば ni( N + 1 ) / 2 となるので、先ほど求めた順位和の期待値 E[Ti] と一致します。また、分散 V[Ti] は U-分布の値と変わらず、

V[Ti] = niui( ni + ui + 1 ) / 12 = ni( N - ni )( N + 1 ) / 12

となります。

一元配置分散分析法のときと同様に、各集団間の変動 SC と集団内の変動 SE を計算すると

SC=Σi{1→k}( ni( mi - m )2 )
=Σi{1→k}( ni{ Ti / ni - ( N + 1 ) / 2 }2 )
=Σi{1→k}( Ti2 / ni - ( N + 1 )Ti + ni( N + 1 )2 / 4 )
=Σi{1→k}( Ti2 / ni ) - ( N + 1 )Σi{1→k}( Ti ) + ( N + 1 )2Σi{1→k}( ni ) / 4
=Σi{1→k}( Ti2 / ni ) - ( N + 1 )・N( N + 1 ) / 2 + ( N + 1 )2N / 4
=Σi{1→k}( Ti2 / ni ) - N( N + 1 )2 / 4
=Σi{1→k}( Ti2 / ni ) - T2 / N
SE=ST - SC
=N( N2 - 1 ) / 12 - Σi{1→k}( Ti2 / ni ) + N( N + 1 )2 / 4
=N( N + 1 )( 2N + 1 ) / 6 - Σi{1→k}( Ti2 / ni )

となります。但し、T は順位和の総計を表し、T = Σi{1→k}( Ti ) の関係が成り立ちます。この式から、順位和の平方を要素数で割った値の和 Σi{1→k}( Ti2 / ni ) が大きいほど集団間の変動は大きく、集団内の変動は小さくなることが分かります。SC と SE の期待値は

E[SC]=Σi{1→k}( E[Ti2] / ni ) - N( N + 1 )2 / 4
=Σi{1→k}( ( V[Ti] + E[Ti]2 ) / ni ) - N( N + 1 )2 / 4
=Σi{1→k}( ( N - ni )( N + 1 ) / 12 + ni( N + 1 )2 / 4 ) - N( N + 1 )2 / 4
={ kN - Σi{1→k}( ni ) }( N + 1 ) / 12 + N( N + 1 )2 / 4 - N( N + 1 )2 / 4
=( k - 1 )N( N + 1 ) / 12
E[SE]=E[ST] - E[SC]
=N( N2 - 1 ) / 12 - ( k - 1 )N( N + 1 ) / 12
=N( N + 1 ){ ( N - 1 ) - ( k - 1 ) } / 12
=( N - k )N( N + 1 ) / 12

と求められます。ここまでの内容は、全要素の順位に対する組み合わせが全て等しい確率で発生することを前提としていることに注意してください。もし、ある集団の値が他と比べて差があり、順位に偏りが生ずるようなときは、この結果は成り立たなくなります。

もし、N 個の要素の順位が全て互いに独立だと仮定すれば、順位の不偏分散 σ2

σ2 = ST / ( N - 1 ) = N( N + 1 ) / 12

なので、全要素の順位に対する組み合わせが全て等しい確率で発生するならば、E[SC] = ( k - 1 )σ2 になります。一元配置分散分析法では SC / σ2 が自由度 k - 1 の χ2-分布に従うことを示しました。各要素の順位が互いに独立で、正規分布 N( ( N + 1 ) / 2, N( N + 1 ) / 12 ) に従うと仮定すれば、順位和に対しても同様に

SC / σ2={ Σi{1→k}( Ti2 / ni ) - T2 / N } / { N( N + 1 ) / 12 }
={ 12 / N( N + 1 ) }Σi{1→k}( Ti2 / ni ) - 3( N + 1 )

が自由度 k - 1 の χ2-分布に従うことになります。この値は通常 H で表して検定を行うため、「クラスカル・ワリスの H 検定(Kruskal-Wallis H Test)」とも呼ばれることがあります。


実際には、各要素の順位は互いに独立ではないので上に示した仮定は成り立ちません。また、全要素の順位に対する組み合わせが全て等しい確率で発生するならば、各要素の順位は 1 から N までのどの値も等しい確率で発生することになるので、各要素の順位は実際には正規分布ではなく離散的な一様分布に従うと考えられます。しかし、前にも述べたとおり、各群の順位和 Ti は U-分布 Pni,ui(U) を ni( ni + 1 ) / 2 だけ正方向へシフトしたものと等しく、N が非常に大きければ、中心極限定理から正規分布 ( ni( N + 1 ) / 2, niui( N + 1 ) / 12 ) に近似することができるので、各集団の順位平均 mi は正規分布 N( ( N + 1 ) / 2, ui( N + 1 ) / 12ni ) に近似できることになります。よって、

zi = √ni{ mi - ( N + 1 ) / 2 } / { ui( N + 1 ) / 12 }1/2

は標準正規分布 N( 0 , 1 ) に従いますが、mi は分散が異なるのでそれを補正するため yi = ( ui / N )1/2zi とすると、この二乗和が自由度 k - 1 の χ2-分布に従うことになって(補足3)、

Σi{1→k}( yi2 )=Σi{1→k}( ni{ mi - ( N + 1 ) / 2 }2 / { N( N + 1 ) / 12 } )
=Σi{1→k}( { Ti2 / ni - ( N + 1 )Ti + ni( N + 1 )2 / 4 } / { N( N + 1 ) / 12 } )
={ Σi{1→k}( Ti2 / ni ) - N( N + 1 )2 / 2 + N( N + 1 )2 / 4 } / { N( N + 1 ) / 12 }
={ 12 / N( N + 1 ) }Σi{1→k}( Ti2 / ni ) - { N( N + 1 )2 / 4 } / { N( N + 1 ) / 12 }
={ 12 / N( N + 1 ) }Σi{1→k}( Ti2 / ni ) - 3( N + 1 )

より、H = SC / σ2 の値と一致します。なお、各要素の順位が互いに独立でなく、正規分布にも従わないことから分かるように、一元配置分散分析法の「SE / σ2 が自由度 N - k の χ2-分布に従う」という性質はこの場合は成り立ちません。


正確な SC / σ2 の確率密度を求めるサンプル・プログラムを以下に示します。

const double DEFAULT_THRESHOLD = 1.0E-6; // デフォルトのしきい値

/*
  浮動小数点数に対して大小比較をするための関数オブジェクト
  差分の絶対値があるしきい値より小さければ、等しいと判断されるようにする
*/
class CmpDbl
{
  double _threshold; // しきい値

public:

  // コンストラクタ
  CmpDbl( double threshold = DEFAULT_THRESHOLD )
    : _threshold( threshold ) {}

  // 比較用関数オブジェクト
  bool operator()( const double& d1, const double& d2 ) {
    if ( fabs( d1 - d2 ) < fabs( _threshold * d1 ) )
      return( false );

    return( d1 < d2 );
  }
};

// 浮動小数点数をキーとする Map
typedef map<double,double,CmpDbl> MapProbDist;

/*
  calcRankSqSum : 順位の二乗の平均の和を求める

  const vector<unsigned int>& perm : 対象の順列
  const vector<unsigned int>& n : 各集合の要素数

  戻り値 : 順位の二乗の平均の和
*/
double calcRankSqSum( const vector<unsigned int>& perm, const vector<unsigned int>& n )
{
  double ans = 0;       // 求める二乗和
  unsigned int idx = 0; // 各要素の読み込み開始位置

  for ( unsigned int i = 0 ; i < n.size() ; ++i ) {
    double t = perm[idx];
    for ( unsigned int j = 1 ; j < n[i] ; ++j )
      t += perm[idx + j];
    ans += t * t / (double)n[i];
    idx += n[i];
  }

  return( ans );
}

/*
  rankSqSumDist : 順位の二乗の平均の和に対する確率分布を求める

  const vector<unsigned int>& n : 各集合の要素数
  MapProbDist& p : 確率変数とその確率密度

  戻り値 : true ... 計算成功 ; false ... 計算に失敗
*/
bool rankSqSumDist( const vector<unsigned int>& n, MapProbDist& p )
{
  typedef MapProbDist::iterator MIt;

  if ( n.size() == 0 ) return( false );
  for ( unsigned int i = 0 ; i < n.size() ; ++i )
    if ( n[i] < 1 ) return( false );

  p.clear(); // 確率分布は初期化

  unsigned int s = sum( n );      // 要素数の総和
  vector<unsigned int> perm( s ); // 順列
  for ( unsigned int i = 1 ; i <= s ; ++i )
    perm[i - 1] = i;
  double permCnt = fact( s );     // 順列の総数

  // 各要素の順位の平均と分散
  double ave = (double)( s * ( s + 1 ) * ( s + 1 ) ) / 4.0;
  double var = (double)( s * ( s + 1 ) ) / 12.0;

  // 順位の二乗の平均の和に対する組み合わせ総数を求める
  do {
    // Sc / var(確率変数)
    double x = ( calcRankSqSum( perm, n ) - ave ) / var;
    MIt it = p.find( x );
    if ( it == p.end() )
      p[x] = 1;
    else
      it->second += 1;
  } while ( std::next_permutation( perm.begin(), perm.end() ) );

  // 和を確率密度に変換
  for ( MIt it = p.begin() ; it != p.end() ; ++it )
    it->second /= permCnt;

  return( true );
}

処理の内容は非常に単純で、順列の先頭から番号順に各集団を割り当てたときの SC / σ2 を求める操作を全順列について繰り返し、それぞれの和が何回出現したかを数えます。順位の二乗和の平均の和は calcRankSqSum で求め、各要素の順位の平均 (ave) と分散 (var) を使って SC / σ2 を計算してそれを map<double,double> 型のキーに利用しています。最後に、和を順列の総数 permCnt で割れば確率密度に変換することができます。
double 型の変数をキーに使う場合、計算時の誤差によって一致するはずのキーが異なるものと判断される場合があります。そのため、ある差の絶対値があるしきい値より小さければ等しいとみなすように、比較用の関数オブジェクト CmpDbl を用意しています。しきい値はコンストラクタで変更できるようにしてあります。

順列を求めるため、Standard Template Library(STL) にある next_permutation を利用しています。最初に 1 から N まで並べた配列を用意しておけば、呼び出すたびに新たな順列を作成し、すべての順列を作成したら戻り値として false を返してくれる優れものです。また、要素の和を求める関数 sum と、階乗を求める関数 fact はその実装については示していませんが、他の章にて紹介していますし、それほど難しい処理でもないので内容は割愛しています。

サンプル・プログラムを使って、いくつかの条件で確率密度と累積確率を求めた結果を以下に示します。

図 4-1. 要素数の変化に対する確率密度の推移
( n1, n2, n3 ) = ( 2, 2, 1 )( n1, n2, n3 ) = ( 2, 2, 2 )
(2,2,1)(2,2,2)
( n1, n2, n3 ) = ( 2, 2, 3 )( n1, n2, n3 ) = ( 2, 2, 4 )
(2,2,3)(2,2,4)
( n1, n2, n3 ) = ( 2, 3, 4 )( n1, n2, n3 ) = ( 3, 3, 3 )
(2,3,4)(3,3,3)
( n1, n2, n3 ) = ( 4, 4, 4 )
(4,4,4)

三つの集団に対して要素数を変化させながら確率密度を求めると、上に示したような結果になりました。要素数が増えるに従って、SC / σ2 の増加により確率密度は小さくなる傾向が見られます。これは、各要素が全てどの順位に対しても等しい確率でなりうる場合、集団間の平均に対するバラツキが大きくなる確率は低くなるということを意味しています。累積密度の分布では χ2-分布と比較できるように両プロットを表示しています。各集団の要素数がある程度大きくなれば、χ2-分布に近似できる様子がこの結果からも分かります。

図 4-2. 集団数の変化に対する確率密度の推移
( n1, n2, n3 ) = ( 2, 2, 2 )( n1, n2, n3, n4 ) = ( 2, 2, 2, 2 )
(2,2,2)(2,2,2,2)
( n1, n2, n3, n4, n5 ) = ( 2, 2, 2, 2, 2 )( n1, n2, n3, n4 ) = ( 3, 3, 3, 3 )
(2,2,2,2,2)(2,2,4)

上図は、集団の数を変化させたときの確率密度の変化を示したグラフです。集団の数が増えても、集団内の数が小さければ χ2-分布との差は大きいままになります。( n1, n2, n3, n4, n5 ) = ( 2, 2, 2, 2, 2 ) の場合、SC / σ2 の確率密度は中央を最大としてほぼ左右対称に分布しています。それに対して ( n1, n2, n3, n4 ) = ( 3, 3, 3, 3 ) の場合は SC / σ2 が増加するに従い確率密度が減少する様子が確認できます。

要素数が極端に少ない場合、正確な確率密度を計算した方が精度は上がりますが、計算できる要素数には限界があります。今回のサンプル・プログラムでは、要素数 10 個を少し超える程度までであれば実用範囲内ですが、15 個で試したときはしばらく待っても処理が完了しませんでした。通常、サンプルとして利用するデータは比較的大きな要素数となるので、χ2-分布で近似しても充分正確な結果が得られます。しかし、要素数の小さなサンプルについても対応できるようにするには、要素数があるしきい値より小さい場合、正確な確率密度を計算してそれを利用するか、検定だけに利用するのであれば、危険率に対する棄却域をあらかじめ配列などの形で保持しておく必要があります。


同順位の要素があった場合、順位の平均を各要素に割り当てるため、同順位の要素がない場合と比較して平均 μ に変化はありませんが、分散 σ2 は一致しなくなります。例えば、五つの要素が全て異なる値ならば、順位は ( 1, 2, 3, 4, 5 ) と付けられるので、その平均は

μ = 5 x ( 5 + 1 ) / 2 = 3

であり、分散は

σ2 = 5 x ( 5 + 1 ) / 12 = 5 / 2

になります。しかし、三番目と四番目の順位にあたる要素が等しい場合、順位は ( 1, 2, 3.5, 3.5, 5 ) になります。このとき、平均は

μ = ( 1 + 2 + 3.5 + 3.5 + 5 ) / 5 = 3

であり、同順位がない場合と等しくなりますが、分散は

σ2 = { ( 1 - 3 )2 + ( 2 - 3 )2 + ( 3.5 - 3 )2 + ( 3.5 - 3 )2 + ( 5 - 3 )2 } / ( 5 - 1 ) = 19 / 8

となって、若干小さな値になります。H は SC を σ2 で割った値なので、同順位が多くなるほど、σ2 が小さくなる分だけ H の値は大きくなるはずです。よって、同順位がある場合はその分 H を補正する必要があります。一般に、m 個の要素が同順位であった場合、その順位が l から l + m - 1 までであれば、平均順位は

Σi{l→l+m-1}( i ) / m = l + Σi{0→m-1}( i ) / m = l + ( m - 1 ) / 2

で求められます。平均が μ ならば、平均との差の平方和は

m{ l + ( m - 1 ) / 2 - μ }2=m{ l + ( m - 1 ) / 2 }2 - 2m{ l + ( m - 1 ) / 2 }μ + mμ2
={ l2m + lm2 - lm + m( m - 1 )2 / 4 } - 2m{ l + ( m - 1 ) / 2 }μ + mμ2

になります。それに対して順位が等しくない場合、平均との差の平方和は

Σi{l→l+m-1}( ( i - μ )2 )=Σi{l→l+m-1}( i2 ) - 2μΣi{l→l+m-1}( i ) + mμ2
=( l + m - 1 ){ ( l + m - 1 ) + 1 }{ 2( l + m - 1 ) + 1 } / 6
 - ( l - 1 ){ ( l - 1 ) + 1 }{ 2( l - 1 ) + 1 } / 6 - 2m{ l + ( m - 1 ) / 2 }μ + mμ2
=( l + m )( l + m - 1 )( 2l + 2m - 1 ) / 6
 - l( l - 1 )( 2l - 1 ) / 6 - 2m{ l + ( m - 1 ) / 2 }μ + mμ2
={ 2( l + m )3 - 3( l + m )2 + ( l + m ) - 2l3 + 3l2 - l } / 6
 - 2m{ l + ( m - 1 ) / 2 }μ + mμ2
=( l2m + lm2 + m3 / 3 - lm - m2 / 2 + m / 6 ) - 2m{ l + ( m - 1 ) / 2 }μ + mμ2

となります。下式の結果と上式の結果の差は

( l2m + lm2 + m3 / 3 - lm - m2 / 2 + m / 6 ) - { l2m + lm2 - lm + m( m - 1 )2 / 4 }
=( m3 / 3 - m2 / 2 + m / 6 ) - ( m3 / 4 - m2 / 2 + m / 4 )
=( m3 - m ) / 12

となって、不偏分散の値としては同順位がある方が ( m3 - m ) / 12( N - 1 ) 小さくなるという結果が得られます。ここで、差の値は順位 l に依存せず、同順位となった要素の数だけで決まることに注意してください。上記のようにまじめに計算しなくても、順位の等しいグループの中央から等間隔で左右に離れた値どうしが、符号が異なるだけで絶対値が等しいという事と、その値が中央値からの距離のみで決まることを考慮すれば、順位の初期値 l には依存しないことが簡単に理解できます。

同順位の要素が複数ある場合も同様の結果になるため、そのようなグループが p 個存在したとすれば

Σi{1→p}( mi3 - mi ) / 12( N - 1 )

但し、mi は i 番目の同順位グループの要素数

だけ値が小さくなることになります。順位に重複のないときの分散は N( N + 1 ) / 12 で求められるので、同順位が存在した場合の分散との比 R は

R={ N( N + 1 ) / 12 - Σi{1→p}( mi3 - mi ) / 12( N - 1 ) } / { N( N + 1 ) / 12 }
=1 - Σi{1→p}( mi3 - mi ) / ( N3 - N )

となり、H を R で割れば補正を行うことができます。


クラスカル・ワリス検定を行うためのサンプル・プログラムを以下に示します。

// 浮動小数点数をキーとする Map
typedef map<double,vector<unsigned int>,CmpDbl> MapRank;
typedef pair<double,vector<unsigned int> > PairRank;

/*
  getRank : データから順位を求める

  const vector< vector<double> >& data : 各集団の要素
  vector< vector<double> >& rank : 求める順位
  double& r : 同順位が存在した場合の補正値
  unsigned int& totalCnt : 全要素数

  戻り値 : true ... 計算成功 ; false ... データが一つもない
*/
bool getRank( const vector< vector<double> >& data, vector< vector<double> >& rank, double& r, unsigned int& totalCnt )
{
  MapRank mapRank; // 昇順にデータを並べるためのバッファ
  totalCnt = 0;    // 全要素数を初期化

  // バッファ mapRank にデータをコピー(同時に並べ替え)
  for ( unsigned int i = 0 ; i < data.size() ; ++i ) {
    totalCnt += data[i].size();
    for ( unsigned int j = 0 ; j < data[i].size() ; ++j ) {
      MapRank::iterator it = mapRank.find( data[i][j] );
      if ( it == mapRank.end() )
        mapRank.insert( PairRank( data[i][j], vector<unsigned int>( 1, i ) ) );
      else
        ( it->second ).push_back( i );
    }
  }

  // 順位の初期化
  rank.resize( data.size() );
  for ( unsigned int i = 0 ; i < rank.size() ; ++i )
    rank[i].clear();

  if ( totalCnt == 0 ) return( false );

  double curRank = 1; // 現在の末尾の順位
  r = 0; // 補正値 r の初期化
  for ( MapRank::const_iterator cit = mapRank.begin() ; cit != mapRank.end() ; ++cit ) {
    const vector<unsigned int>& v = cit->second;
    unsigned int cnt = v.size();
    double avgRank = curRank + (double)( cnt - 1 ) / 2.0; // 平均順位

    // rankに平均順位を登録
    for ( unsigned int i = 0 ; i < cnt ; ++i )
      rank[v[i]].push_back( avgRank );

    // 末尾の順位と補正値を更新
    curRank += v.size();
    r += cnt * ( cnt - 1 ) * ( cnt + 1 );
  }

  // 補正値 r の計算
  double d = totalCnt * ( totalCnt - 1 ) * ( totalCnt + 1 );
  r = ( d - r ) / d;

  return( true );
}

/*
  kruskalWallisTest : Kruskal-Wallis検定

  const vector< vector<double> >& rank : 各集団の順位
  double a : 危険率
  double r : 補正値
  unsigned int totalCnt : 全要素数
  double threshold : binSearchで棄却域を求める時のしきい値

  戻り値 : true ... 帰無仮説が棄却された ; false ... 帰無仮説が保留された
*/
bool kruskalWallisTest( const vector< vector<double> >& rank, double a, double r, unsigned int totalCnt, double threshold )
{
  cout << "***** Kruskal-Wallis Test *****" << endl << endl;

  double h = 0; // H = Sc / var
  unsigned int k = rank.size(); // 集団の数(自由度)
  if ( k < 2 ) {
    cerr << "The size of group must be greater than 1." << endl;
    return( false );
  }

  // Σ( Ti^2 / ni ) の計算
  for ( unsigned int i = 0 ; i < k ; ++i )
    h += pow( sum( rank[i] ), 2 ) / (double)( rank[i].size() );

  // H の計算
  h -= totalCnt * pow( totalCnt + 1, 2 ) / 4.0;
  h /= ( (double)( totalCnt * ( totalCnt + 1 ) ) / 12.0 );
  h /= r;

  cout << "#Class(k) = " << k << " ; #Data(N) = " << totalCnt << endl;
  cout << "H : " << h << " ( R : " << r << " )" << endl << endl;

  ChiSquareDistribution chiDist( k - 1 );
  double rej = binSearch( chiDist, 1.0 - a, threshold ); // 棄却域の定義域
  cout << "Critical Region : >" << rej << endl;
  double p = 1.0 - chiDist.p( 0, h ); // p値
  cout << "p-value : " << p << endl << endl;

  if ( rej < h ) {
    cout << "Null hypothesis is rejected." << endl;
    return( true );
  } else {
    cout << "Null hypothesis is accepted." << endl;
    return( false );
  }
}

データを元に順位を求める処理と検定を行う処理は別々に分けてあり、それぞれ getRankkruskalWallisTest で行います。このようにしたのは、順位のみで実際のデータがない場合にも対応できるようにするためです。

getRank では、最初に map<double,vector<unsigned int>> 型のコンテナ mapRank へ要素と集団の番号を登録します。このとき、等しいと判断されたものは同じキーに登録されるように、その値は vector<unsigned int> としています。map クラスはキーに対して昇順に並んだ形でデータが保持されるので、先頭の要素から順番に値(すなわち各集団の番号からなる配列)を取得すれば、そのまま順位を決定することができます。
各集団の番号からなる配列への参照を v に代入し、その要素数を使って平均順位 avgRank を計算した後、順位を登録する配列 rank にその値を登録します。この処理方法から理解できるように、各集団の要素とその順位は一対一に対応しているわけではなく、順位は昇順に登録されることに注意してください。各要素に対応して順位が並ぶようにする場合は、処理方法を見直す必要があります。
順位に重複があるかどうかは、配列への参照 v の要素数から分かるので、これを利用して同時に補正値 r を求めています。また、要素数の合計も検定時に必要になるので、この処理の中でいっしょに求めて totalCnt へ登録しています。

kruskalWallisTest へは、getRank で求めた各データをそのまま渡すことを想定して引数を決めています。ここで、順位だけがデータとして分かっていたとして、その順位を使って検定を行う場合、順位に重複がないのであれば r に 1 を代入しておけば正しく処理ができます(重複があるのであれば、他の場所で補正値を計算しておく必要があります)。また、要素数はあらかじめ計算しておく必要があります。処理自体は非常に単純で、H = SC / σ2 を求めて変数 h へ代入し、自由度 k - 1 の χ2-分布から棄却域を求めて両者を比較して、その結果を返します。χ2-分布用のクラスは「(6) 標本分布」のカイ二乗分布用サンプル・プログラムにて紹介しています。また、検定に利用している binSearch は「(8) 推定」の中の「平均値の区間推定用サンプル・プログラム」をご覧ください。


最初に示した、肥料に対するキャベツの質量のデータを使って、危険率 5% でクラスカル・ワリス検定を行った結果を以下に示します。"#Class(k)" は群の数、"#Data(N)" が全要素数で、"H" が H 値 ( = SC / σ2 )、"R" が同順位があった場合の補正値、"Critical Region" は棄却域で、"p-value" が p 値をそれぞれ表しています。

***** Kruskal-Wallis Test *****

#Class(k) = 3 ; #Data(N) = 15
H : 5.78157 ( R : 0.998214 )

Critical Region : >5.99146
p-value : 0.0555325

Null hypothesis is accepted.

一元配置分散分析での結果と異なり、同じ危険率でもここでは仮説が保留となっています。といっても、保留となっているだけで仮説が正しいというわけではありません。p 値を見ると 5.55 % 程度と低い値なので、データに偏りがないとは言えないと判断することは充分にできます。

図 4-3. 実測値と順位の分布比較
実測値順位
キャベツの質量に対する分布 キャベツの質量に対する順位の分布

上の図は、キャベツの質量に対する度数分布を集団別に色分けしたものを、実測値とその順位のそれぞれで表したものです。比較しやすいように、順位は質量の小さいものから番号付け(すなわち昇順と)しています。「肥料なし」と「肥料 A」を比較すると、実測値も順位も分布に偏りが見られます。

データに正規性が見られない場合や、分散が等しくないような場合は「一元配置分散分析法」ではなく「クラスカル・ワリス検定」を利用するのが一般的な考え方のようです。「マン・ホイットニーの U 検定(Mann-Whitney U Test)」を「分布が等しいか」判定する目的で利用する場合、二つの分布の分散が等しくなくても「等しい」と判断される場合がありましたが、「順位平均値が等しいか」判定するのであれば分散が等しいかどうかを気にする必要はなくなります。それは「クラスカル・ワリス検定」でも当てはまり、各分布の分散を気にすることは不要になります。なお、処理の内容を見ても明らかなように、平均値の差の有意性を判定する場合、「マン・ホイットニーの U 検定」は「クラスカル・ワリス検定」で集合が二つだった場合と同じものになります。
ここで気を付けなければならないのは、判定対象は「順位平均」であって「データの平均」ではないという点です。順位平均を対象としていることから、実際の値の分布に対しては「平均値」ではなく「中央値」が一致するかどうかを検定していることになります。従って、分布が対称でない場合は平均値と中央値は一致しないことから、検定結果が差の有意性を示したとしても、平均が一致しないと判断できるとは限らないことになります。例えば、値の大きい側に裾が広がった分布と小さい側に裾が広がった分布とで検定を行った結果、差に有意性があると判断されたとしても、平均値は中央値よりも裾の広がった側にシフトするため有意性がなくなる可能性があります。データから順位を求めてクラスカル・ワリス検定を行う場合は、あくまでも中央値に対する検定であって、平均値に対しては分布の対称性が示されない限り判断できないことを理解しておく必要があります。

*4-1)(10) 順序統計量」の「3) マン・ホイットニーの U 検定(Mann-Whitney U Test)」を参照


5) フリードマン検定(Friedman Test)

クラスカル・ワリス検定では、全ての集合をまとめて順位付けを行いましたが、ここで各集合の中で順位付けを行うことを考えます。例として、「二元配置分散分析法」で示した作業者毎の作業時間の時系列推移を使って具体的に説明します。各作業者ごとに、作業時間の短いものから順位付けをすると、以下のような結果になります。

表 5.1 作業者ごとの作業時間の順位
作業日1234567平均
作業者実測値順位実測値順位実測値順位実測値順位実測値順位実測値順位実測値順位実測値順位
A65.3760.8659.8559.4450.82.550.82.550.6156.84
B63.1759.6659.3556.2446.3148.9251.0354.94
C65.3760.1657.2552.9353.9449.6248.4155.34
D61.3759.7652.0353.1547.5150.4252.8453.84
E62.5756.1654.0552.73.547.7148.1252.73.553.44
平均63.5759.3656.54.654.93.949.21.949.62.151.12.554.94

作業者ごとに作業時間の順位付けをしても、1, 2 日目の作業時間は下位の順位になり、日が経つにつれて順位が上がる様子は実績値の場合と変わりません。作業者ごとに順位付けしているので、作業者ごとの順位平均は全て 4 になる反面、日毎の順位平均はそれぞれバラバラになっていることから、作業者が工具に慣れて作業時間が短縮していると予想することができます。これを定量的に表すための検定法が「フリードマン検定 (Friedman Test)」です。

「(反復のない)二元配置分散分析法」では、全体の変動を行間・列間の変動と誤差変動に分解してそれぞれの変動を調べる手法でした。「フリードマン検定」は、行間(または列間)の変動をわざとゼロにした上で、全体の変動を列間(または行間)変動と誤差変動に分解する手法と考えることができます。「(反復のない)二元配置分散分析法」における、全体の変動 ST、行間・列間変動 SR, SC、誤差変動 SE は次式で求められるのでした。

ST=Σr{1→l}( Σc{1→k}( ( xrc - m )2 ) )
SR=r{1→l}( ( mr. - m )2 )
SC=c{1→k}( ( m.c - m )2 )
SE=Σr{1→l}( Σc{1→k}( { ( xrc - m ) - ( mr. - m ) - ( m.c - m ) }2 ) )
=Σr{1→l}( Σc{1→k}( ( xrc - mr. - m.c + m )2 ) )

「フリードマン検定」において、行間(または列間)変動を SG で表すとすれば、もう一方の(列間または行間の)変動はゼロなので、上式は次のように表すことができます。

ST=Σg{1→k}( Σi{1→n}( ( Rgi - m )2 ) )
SG=g{1→k}( ( mg - m )2 )
SE=Σg{1→k}( Σi{1→n}( { ( Rgi - m ) - ( mg - m ) }2 ) )
=Σg{1→k}( Σi{1→n}( ( Rgi - mg )2 ) )

但し、m は全体の順位平均、mg は行または列ごとの順位平均、Rgi は列間(または行間)の順位をそれぞれ表します。ここでも ST = SG + SE は成り立ち、mg と m の差が小さい、すなわち各行または各列間で順位に差がほとんどなければ誤差のほとんどは SE に含まれ、逆に各行または各列間での順位差が大きくなるほど SG が大きくなります。

各行ごとに 1 から k までの順位を付けたと仮定します。このとき、ある行の順位和は 1 から k までの和 k( k + 1 ) / 2 なので、Rgi の期待値 E[Rgi] はこの順位和を k で割った値と等しく

E[Rgi] = ( k + 1 ) / 2

となります。また、Rgi は U-分布 P1,k-1(U) に対して変数変換 Rgi = U + 1・( 1 + 1 ) / 2 = U + 1 を行った分布と等しいと考えられるので ( 4) クラスカル・ワリス検定(Kruskal-Wallis Test) 参照 )、分散 V[Rgi] は

V[Rgi] = ( k + 1 )( k - 1 ) / 12

になります。各行の分布は互いに独立であると見なせるので、ある列の順位和を Ti としたとき、その期待値と分散 E[Ti], V[Ti] はそれぞれ E[Rgi], V[Rgi] の和で表され

E[Ti] = Σg{1→n}( E[Rgi] ) = n( k + 1 ) / 2

V[Ti] = Σg{1→n}( V[Rgi] ) = n( k + 1 )( k - 1 ) / 12

となって、

zi = { Ti - n( k + 1 ) / 2 } / { n( k + 1 )( k - 1 ) / 12 }1/2

としたとき、クラスカル・ワリス検定の場合と同様な考え方によって、{ ( k - 1 ) / k }1/2 で zi を補正した上で二乗和を計算した結果は自由度 k - 1 の χ2-分布に従います(補足4)。これを計算すると、

Σi{1→k}( { ( k - 1 ) / k }zi2 )
={ ( k - 1 ) / k }Σi{1→k}( { Ti - n( k + 1 ) / 2 }2 / { n( k + 1 )( k - 1 ) / 12 } )
=Σi{1→k}( ( nmi - nm )2 ) / { nk( k + 1 ) / 12 }
=i{1→k}( ( mi - m )2 ) / { k( k + 1 ) / 12 }
=SG / { k( k + 1 ) / 12 }

という結果が得られます。但し、m は全体の平均を表し、kn 個の要素の順位和が nk( k + 1 ) / 2 となることから m = ( k + 1 ) / 2 という結果になることを利用しています。クラスカル・ワリス検定の場合と同じように、この式は別の形で表すことができて、

SG / { k( k + 1 ) / 12 }
=Σi{1→k}( { Ti - n( k + 1 ) / 2 }2 ) / { nk( k + 1 ) / 12 }
=Σi{1→k}( Ti2 - n( k + 1 )Tg + { n( k + 1 ) / 2 }2 ) / { nk( k + 1 ) / 12 }
=[ Σi{1→k}( Ti2 ) - n( k + 1 )・nk( k + 1 ) / 2 + k{ n( k + 1 ) / 2 }2 ] / { nk( k + 1 ) / 12 }
={ Σi{1→k}( Ti2 ) - n2k( k + 1 )2 / 4 } / { nk( k + 1 ) / 12 }
={ 12 / nk( k + 1 ) }Σi{1→k}( Ti2 ) - 3n( k + 1 )

で計算することもできます。ところで、一つの群(一行)の中での不偏分散は

Σi{1→k}( { i - ( k + 1 ) / 2 }2 ) / ( k - 1 )
={ Σi{1→k}( i2 ) - k( k + 1 )2 / 4 } / ( k - 1 )
={ k( k + 1 )( 2k + 1 ) / 6 - k( k + 1 )2 / 4 } / ( k - 1 )
=k( k + 1 ){ 2( 2k + 1 ) - 3( k + 1 ) } / 12( k - 1 )
=k( k + 1 )( k - 1 ) / 12( k - 1 ) = k( k + 1 ) / 12

と計算できるので、σ2 = k( k + 1 ) / 12 として、以下は SG / σ2 で表します。ここで、今回求めた分散と、先ほど U-分布を利用して求めた V[Rgi] が異なるのは、σ2 は各要素が独立であると仮定した場合の値であるからです。実際には、最後の要素の順位は他の要素の順位が確定した時点で自動的に決まってしまうので、自由度を一つ落とすことで k( k + 1 ) / 12 から ( k - 1 )( k + 1 ) / 12 に変化したと考えれば理解しやすいのではないかと思います。また、全ての群をまとめたときの不偏分散は、各群の不偏分散の和で表されるので nk( k + 1 ) / 12 となります。


SG / σ2 の正確な確率分布を求めるためのサンプル・プログラムを以下に示します。

/*
  calcRankSqSum : permList内のgNum[i]で指定された順列の各要素毎の和を計算し、その二乗和を求める

  permListの第一添字を行、第二添字を列としたとき
   gNum[i] は i 番目の行の順列の番号を表し
   permList[gNum[i]] は i 番目の行の順列を表し
   permList[gNum[i]][j] は 対象の順列の j 番目の列の要素を表す

  const vector< vector<unsigned int> >& permList : 順列のリスト
  const vector<unsigned int>& gNum : 各行の順列の番号

  戻り値 : 求めた二乗和
*/
double calcRankSqSum( const vector< vector<unsigned int> >& permList, const vector<unsigned int>& gNum )
{
  unsigned int k = permList[0].size(); // 群内の要素数(列数)
  unsigned int n = gNum.size();        // 群の数(行数)

  double sg = 0; // 求める二乗和 (Sg)
  for ( unsigned int j = 0 ; j < k ; ++j ) {
    unsigned int d = 0; // 列 j の要素の和
    for ( unsigned int i = 0 ; i < n ; ++i )
      d += permList[gNum[i]][j];
    sg += pow( d, 2 );
  }

  return( sg );
}

/*
  next_gNum : 順列の次の組み合わせを求める

  const vector<unsigned int>& gNum : 各行の順列の番号
  unsigned int index : 現在着目している最小添字
  max_gNum : 順列の個数(指定できる順列の番号の最大値)

  戻り値 : true ... 次の組み合わせがあった ; false ... これ以上組み合わせは存在しない
*/
bool next_gNum( vector<unsigned int>& gNum, unsigned int index, unsigned int max_gNum )
{
  // 最小添字が順列の番号リストのサイズを超えたらこれ以上組み合わせはない
  if ( index >= gNum.size() ) return( false );

  // index が指す番号が最大値より小さい場合
  if ( gNum[index] < max_gNum - 1 ) {
    ++( gNum[index] );
    return( true );
  }

  // index が指す番号が最大値に達した場合
  // 番号をリセットして次の添字の値を更新
  gNum[index] = 0;

  return( next_gNum( gNum, index + 1, max_gNum ) );
}

/*
  friedmanDist : フリードマン検定の正確な確率分布を求める

  unsigned int k : 群内の要素数(列数)
  unsigned int n : 群の数(行数)
  MapProbDist& p : 確率変数とその確率密度

  戻り値 : true ... 計算成功 ; false ... 計算に失敗
*/
bool friedmanDist( unsigned int k, unsigned int n, MapProbDist& p )
{
  typedef MapProbDist::iterator MIt;

  if ( k == 0 || n == 0 ) return( false );

  // k 個の順列を初期化
  vector<unsigned int> perm( k );
  for ( unsigned int i = 1 ; i <= perm.size() ; ++i )
    perm[i - 1] = i;
  // 順列のリストを作成
  vector< vector<unsigned int> > permList;
  do {
    permList.push_back( perm );
  } while ( next_permutation( perm.begin(), perm.end() ) );

  vector<unsigned int> gNum( n, 0 ); // 各行の順列番号

  p.clear(); // 確率分布は初期化

  double offset = (double)( n * ( k + 1 ) ) * 3.0;   // オフセット値( 3n( k + 1 ) )
  double var = (double)( n * k * ( k + 1 ) ) / 12.0; // 分散値( nk( k + 1 ) / 12 )

  // 計算開始
  do {
    // Sg / var - offset (確率変数)
    double x = calcRankSqSum( permList, gNum ) / var - offset;
    MIt it = p.find( x );
    if ( it == p.end() )
      p[x] = 1;
    else
      it->second += 1;
  } while ( next_gNum( gNum, 0, permList.size() ) );

  // 和を確率密度に変換(組み合わせ総数は ( k! )^n )
  for ( MIt it = p.begin() ; it != p.end() ; ++it )
    it->second /= pow( fact( k ), n );

  return( true );
}

処理の方法はかなり強引で、順列の全ての組み合わせを求めた上で、各群(ここでは各行で表しています)に対する順列の組み合わせ全てに対して確率変数 SG / σ2 を計算しながら、クラスカル・ワリス検定の確率分布を求めるときに利用した MapProbDist 型のインスタンスに登録します。最終的には各確率変数の個数が得られるので、これを全組み合わせ総数 ( これは 群の数を n、各群の要素数を k としたとき ( k! )n で求められます ) で割って確率密度としています。


k を群内の要素数、n を群の数として、いくつかの条件で SG / σ2 の確率分布を求めた結果のグラフを以下に示します。

図 5-1. 行列数の変化に対する確率密度の推移
( k, n ) = ( 3, 2 )( k, n ) = ( 3, 3 )
(3,2)(3,3)
( k, n ) = ( 3, 5 )( k, n ) = ( 4, 2 )
(3,5)(4,2)
( k, n ) = ( 4, 3 )( k, n ) = ( 4, 5 )
(4,5)(4,5)
( k, n ) = ( 5, 2 )( k, n ) = ( 5, 3 )
(5,2)(5,3)
( k, n ) = ( 5, 5 )
(5,5)

行・列数のいずれかが 5 より小さい場合は正確に求めたほうが精度は上がりますが、どちらも 5 以上であれば χ2-分布で近似しても問題はなさそうです。統計学の書籍などに載せてある検定表を見ると、たいていは群の数が 5 から 6 個、群内の要素数が 10 個程度を限度として棄却域を用意しているようなので、それ以上の場合は χ2-分布で近似するのが一般的であると考えられます。


各群の中で同順位があった場合、クラスカル・ウォリス検定の場合と同様に σ2 を補正する必要があります。考え方はクラスカル・ウォリス検定の場合と全く同じで、それぞれの群の中で同順位になるグループが全部合わせて mi ( i = 1, 2, ... p ) 個あったとき、同順位がない場合と比べたときの平方和の差の総和は

Σi{1→p}( mi3 - mi ) / 12

になります。順位に重複がないときの不偏分散は nk( k + 1 ) / 12 なので、補正値 R は

R={ nk( k + 1 ) / 12 - Σi{1→p}( ( mi3 - mi ) / 12( k - 1 ) ) } / { nk( k + 1 ) / 12 }
=1 - Σi{1→p}( mi3 - mi ) / n( k3 - k )

という結果になります。


フリードマン検定のサンプル・プログラムを以下に示します。

/*
  getRank : データから順位を求める(Friedman検定用)

  const vector< vector<double> >& data : 各集団の要素
  vector< vector<double> >& rank : 求める順位
  double& r : 同順位が存在した場合の補正値

  戻り値 : true ... 計算成功 ; false ... データが一つもない ; 要素数が等しくない群が存在する
*/
bool getRank( const vector< vector<double> >& data, vector< vector<double> >& rank, double& r )
{
  r = 0; // 補正値 r の初期化

  unsigned int n = data.size();
  if ( n < 1 ) {
    cerr << "The group count must be greater than 0." << endl;
    return( false );
  }
  unsigned int k = data[0].size();
  for ( unsigned int i = 1 ; i < n ; ++i ) {
    if ( data[i].size() != k ) {
      cerr << "The size of all groups must be equal." << endl;
      return( false );
    }
  }

  // 順位の初期化
  rank.resize( n );
  for ( unsigned int i = 0 ; i < rank.size() ; ++i )
    rank[i].resize( k );

  // バッファ mapRank にデータをコピー(同時に並べ替え)
  for ( unsigned int i = 0 ; i < n ; ++i ) {

    MapRank mapRank; // 昇順にデータを並べるためのバッファ
    for ( unsigned int j = 0 ; j < k ; ++j ) {
      MapRank::iterator it = mapRank.find( data[i][j] );
      if ( it == mapRank.end() )
        mapRank.insert( PairRank( data[i][j], vector<unsigned int>( 1, j ) ) );
      else
        ( it->second ).push_back( j );
    }

    double curRank = 1; // 現在の末尾の順位
    for ( MapRank::const_iterator cit = mapRank.begin() ; cit != mapRank.end() ; ++cit ) {
      const vector<unsigned int>& v = cit->second;
      unsigned int cnt = v.size();
      double avgRank = curRank + (double)( cnt - 1 ) / 2.0; // 平均順位

      // rankに平均順位を登録
      for ( unsigned int j = 0 ; j < cnt ; ++j )
        rank[i][v[j]] = avgRank;

      // 末尾の順位と補正値を更新
      curRank += v.size();
      r += cnt * ( cnt - 1 ) * ( cnt + 1 );
    }
  }

  // 補正値 r の計算
  double d = n * k * ( k + 1 ) * ( k - 1 );
  r = ( d - r ) / d;

  return( true );
}

/*
  friedmanTest : Friedman検定

  const vector< vector<double> >& rank : 各集団の順位
  double a : 危険率
  double r : 補正値
  double threshold : binSearchで棄却域を求める時のしきい値

  戻り値 : true ... 帰無仮説が棄却された ; false ... 帰無仮説が保留された
*/
bool friedmanTest( const vector< vector<double> >& rank, double a, double r, double threshold )
{
  cout << "***** Friedman Test *****" << endl << endl;

  unsigned int n = rank.size();
  if ( n < 1 ) {
    cerr << "The group count must be greater than 0." << endl;
    return( false );
  }
  unsigned int k = rank[0].size();
  for ( unsigned int i = 1 ; i < n ; ++i ) {
    if ( rank[i].size() != k ) {
      cerr << "The size of all groups must be equal." << endl;
      return( false );
    }
  }

  // Σ( Tg^2 / var ) の計算
  double chi = 0;
  for ( unsigned int i = 0 ; i < k ; ++i ) {
    double t = 0;
    for ( unsigned int j = 0 ; j < n ; ++j )
      t += rank[j][i];
    chi += pow( t, 2 );
  }

  // chi の計算
  chi /= (double)( n * k * ( k + 1 ) ) / 12.0;
  chi -= (double)( n * ( k + 1 ) ) * 3.0;
  chi /= r;

  cout << "#Class(k) = " << k << " ; #Group(n) = " << n << endl;
  cout << "chi^2 : " << chi << " ( R : " << r << " )" << endl << endl;

  ChiSquareDistribution chiDist( k - 1 );
  double rej = binSearch( chiDist, 1.0 - a, threshold ); // 棄却域の定義域
  cout << "Critical Region : >" << rej << endl;
  double p = 1.0 - chiDist.p( 0, chi ); // p値
  cout << "p-value : " << p << endl << endl;

  if ( rej < chi ) {
    cout << "Null hypothesis is rejected." << endl;
    return( true );
  } else {
    cout << "Null hypothesis is accepted." << endl;
    return( false );
  }
}

getRank はデータから順位を求めるための関数ですが、クラスカル・ワリス検定の場合と異なり各群ごとに順位を求める必要があるので、前に作成したものとは別に用意しています。しかし、異なる部分は、順位の求め方と補正値の計算式くらいなので、クラスカル・ワリス検定のときに用意した関数をある程度流用して作成することができます。friedmanTest がフリードマン検定用のメイン・ルーチンで、カイ二乗値を変数 chi で求めたら、それを使って検定を行います。


最初に示した、「二元配置分散分析法」での作業者毎の作業時間の時系列推移を、フリードマン検定を使って危険率 5% で検定した結果を以下に示します。"#Class(k)" は群内の要素数、"#Group(n)" が群の数で、"chi^2" が χ2 値 ( = SG / σ2 )、"R" が同順位があった場合の補正値、"Critical Region" は棄却域で、"p-value" が p 値をそれぞれ表しています。

***** Friedman Test *****

#Class(k) = 7 ; #Group(n) = 5
chi^2 : 25.5108 ( R : 0.992857 )

Critical Region : >12.5916
p-value : 0.00027454

Null hypothesis is rejected.

順位に置き換えた場合でも、仮説は棄却されて「時系列変動が見られる」と判定されています。下のグラフは、順位に置き換えたときの時系列推移を表していますが、順位に変換しても、日が経つにつれて作業時間が短くなる傾向を見ることができます。

順位の時系列変動

今回は、ANOVA をテーマとして取り上げてみましたが、予想以上に内容が多くなり、かなりのボリュームとなってしまいました。統計解析用のソフトウェアには必ずと言っていいほど利用できるくらい非常に有名な検定法なので、一般的にもよく知られており利用頻度も非常に高いと思います。しかし、制約が割と多く、実際には利用できない場面もあるので、順序統計量を利用したクラスカル・ワリス検定やフリードマン検定といった手法も知っておくと、さらに活用範囲は広くなるのではないでしょうか。


補足1) SC / σ2 が自由度 k - 1 の χ2-分布に従うことの証明

SC / σ2 は次のように表されるのでした。

SC / σ2 = Σi{1→k}( ni( mi - m )2 ) / σ2

但し、

m = Σi{1→k}( nimi ) / N

N = Σi{1→k}( ni )

になります。これによく似たものとして、xi が正規分布 N( μ, σ2 ) に従うときに、y = Σi{1→N}( ( xi - m )2 ) / σ2 = Ns2 / σ2 が自由度 N - 1 の χ2-分布に従うことを前に示しましたが (「(6) 標本分布」の「補足2)」を参照 )、このときと同じやり方で証明することができます。

まず、SC の式を次のように変形します。

SC / σ2=Σi{1→k}( ( √nimi - √nim )2 ) / σ2
=Σi{1→k}( ( √ni( mi - μ ) / σ - √ni( m - μ ) / σ )2 )

ここで、μ は各集団に共通な平均値を表します。つまり、どの母集団も平均は等しい(そして最初の仮定から分散も等しい)としていることに注意してください。このとき、mi は正規分布 N( μ, σ2 / ni ) に従うことから、それを正規化した √ni( mi - μ ) / σ は標準正規分布 ( 0, 1 ) に従うことになります。この正規化した値を yi で表し、( m - μ ) / σ = my とすれば、

SC / σ2 = Σi{1→k}( yi - √nimy )2 )

になります。y = ( y1, y2, ... yk )T を任意の直交行列 Qz = Qy と変換した変数 z = ( z1, z2, ... zk )T は、やはり標準正規分布 N( 0, 1 ) に従います。Q の第 k 行を uk = ( ( n1 / N )1/2, ( n2 / N )1/2, ... ( nk / N )1/2 ) としたとき、

||uk||2 = Σi{1→k}( ni / N ) = 1

であり、Q は直交行列として成り立ちます。このとき、

zk=Σi{1→k}( √niyi / √N )
=Σi{1→k}( ni( mi - μ ) ) / √Nσ
=N( m - μ ) / √Nσ = √Nmy

となるので、

Σi{1→k}( ( yi - √nimy )2 )=Σi{1→k}( yi2 - 2√nimyyi + nimy2 )
=Σi{1→k}( yi2 - 2myni( mi - μ ) / σ + nimy2 )
=Σi{1→k}( yi2 ) - 2myΣi{1→k}( ni( mi - μ ) ) / σ + Nmy2
=Σi{1→k}( zi2 ) - 2my・N( m - μ ) / σ + Nmy2
=Σi{1→k-1}( zi2 ) + zk2 - 2Nmy2 + Nmy2
=Σi{1→k-1}( zi2 ) + Nmy2 - Nmy2
=Σi{1→k-1}( zi2 )

と計算できます。直交行列による変換ではノルムは変化しないので、Σi{1→k}( yi2 ) = Σi{1→k}( zi2 ) が成り立つことに注意してください (「固有値問題 (1) 対称行列の固有値」の「3) 直交行列(Orthogonal Matrix)」参照)。よって、SC は標準正規分布に従う k - 1 個の独立な確率変数の二乗和で表され、自由度 k - 1 の χ2-分布に従うことになります。


補足2) SE / σ2 が自由度 ( k - 1 )( l - 1 ) の χ2-分布に従うことの証明

二元配置分散分析法のモデルにおいて、SE / σ2 は次のように表されるのでした。

SE / σ2 = Σr{1→l}( Σc{1→k}( ( xrc - mr. - m.c + m )2 ) ) / σ2

まず、yk(r-1)+c = ( xrc - mr. - m.c + m ) / σ とすると、E[yk(r-1)+c] は

E[yk(r-1)+c]=E[( xrc - mr. - m.c + m ) / σ]
=( μrc - μr. - μ.c + μ ) / σ
=[ { μ + ( μr. - μ ) + ( μ.c - μ ) } - μr. - μ.c + μ ] / σ
=0

になり、V[yk(r-1)+c] は

V[yk(r-1)+c]=E[( xrc - mr. - m.c + m )2 / σ2]
=E[xrc2 + mr.2 + m.c2 + m2 -2mr.xrc - 2m.cxrc + 2mxrc + 2mr.m.c + -2mmr. - 2mm.c] / σ2

と展開できるので、各項を計算すると、

E[xrc2] = σ2 + μrc2

E[mr.2] = σ2 / k + μr.2

E[m.c2] = σ2 / l + μ.c2

E[m2] = σ2 / kl + μ2

はすでに二元配置分散分析法の説明の中で求めてあり、E[mr.xrc] は

E[mr.xrc]=E[xrcΣi{1→k}( xri ) / k]
=E[xrc2 / k + xrcΣi{1→k;i≠c}( xri ) / k]
=E[xrc2] / k + E[xrci{1→k;i≠c}( E[xri] ) / k
=( σ2 + μrc2 ) / k + μrcΣi{1→k;i≠c}( μri ) / k
=σ2 / k + μrcΣi{1→k}( μri ) / k
=σ2 / k + μrcμr.

と計算することができます。mr. は xri ( i = 1, 2, ... k ) の平均であり、xrc との積を考えると i = c の場合だけ互いに独立ではなく期待値は E[xrc2] になるので分離して計算をしています。それ以外 ( i ≠ c の場合 ) は互いに独立なので E[xrc]E[xri] = μrcμri となって、上記のような結果になります。同様な方法で

E[m.cxrc] = σ2 / l + μrcμ.c

E[mxrc] = σ2 / kl + μrcμ

と求めることができます。他の項に関しても、

E[mr.m.c]=E[Σi{1→k}( xrij{1→l}( xjc ) / kl]
=E[Σi{1→k}( xriΣj{1→l}( xjc ) ) / kl]
={ σ2 + Σi{1→k}( μriΣj{1→l}( μjc ) ) } / kl
=σ2 / kl + Σi{1→k}( μriμ.c ) / k
=σ2 / kl + μr.μ.c
E[mmr.]=E[mr.Σj{1→l}( mj. ) / l]
={ E[mr.2] + E[mr.j{1→l;j≠r}( E[mj.] ) } / l
={ ( σ2 / k + μr.2 ) + μr.Σj{1→l;j≠r}( μj. ) } / l
=σ2 / kl + μr.Σj{1→l}( μj. ) / l
=σ2 / kl + μμr.
E[mm.c]=σ2 / kl + μμ.c

になり、V[yk(r-1)+c] は

V[yk(r-1)+c]=( σ2 + μrc2 ) + ( σ2 / k + μr.2 ) + ( σ2 / l + μ.c2 ) + ( σ2 / kl + μ2 )
 - 2( σ2 / k + μrcμr. ) - 2( σ2 / l + μrcμ.c ) + 2( σ2 / kl + μrcμ )
 + 2( σ2 / kl + μr.μ.c ) - 2( σ2 / kl + μμr. ) - 2( σ2 / kl + μμ.c )
=( 1 - 1 / k - 1 / l + 1 / kl )σ2 + μrc2 + μr.2 + μ.c2 + μ2
 - 2μrcμr. - 2μrcμ.c + 2μrcμ + 2μr.μ.c - 2μμr. - 2μμ.c
=( k - 1 )( l - 1 )σ2 / kl + ( μrc - μr. - μ.c + μ )2
=( k - 1 )( l - 1 )σ2 / kl

という結果になります (この結果から、yk(r-1)+c は正規分布 N( 0, ( k - 1 )( l - 1 )σ2 / kl ) に従うことになります)。yk(r-1)+c と yk(s-1)+d の共分散 E[yk(r-1)+cyk(s-1)+d] は

E[yk(r-1)+cyk(s-1)+d]=E[( xrc - mr. - m.c + m )( xsd - ms. - m.d + m ) / σ2]
=E[xrcxsd + mr.ms. + m.cm.d + m2 - mr.xsd - ms.xrc - m.cxsd - m.dxrc
 + mxrc + mxsd + mr.m.d + ms.m.c - mmr. - mms. - mm.c - mm.d] / σ2

と展開できて、

E[xrcxsd]=E[xrc]E[xsd] = μrcμsd
E[mr.ms.]=E[Σc{1→k}( xrcd{1→k}( xsd )] / k2
=E[Σc{1→k}( Σd{1→k}( xrcxsd ) )] / k2
=Σc{1→k}( Σd{1→k}( E[xrc]E[xsd] ) ) / k2
=Σc{1→k}( Σd{1→k}( μrcμsd ) ) / k2
=Σc{1→k}( μrcμs. ) / k
=μr.μs.
E[m.cm.d]=μ.cμ.d
E[mr.xsd]=E[xsdΣc{1→k}( xrc )] / k
=E[xsdc{1→k}( E[xrc] ) / k
=μsdΣc{1→k}( μrc ) / k
=μr.μsd
E[ms.xrc]=μs.μrc
E[m.cxsd]=μ.cμsd
E[m.dxrc]=μ.dμrc

と各項を求めることができるので、

E[yk(r-1)+cyk(s-1)+d]={ μrcμsd + μr.μs. + μ.cμ.d + ( σ2 / kl + μ2 ) - μr.μsd - μs.μrc - μ.cμsd - μ.dμrc
 + ( σ2 / kl + μrcμ ) + ( σ2 / kl + μsdμ ) + ( σ2 / kl + μr.μ.d ) + ( σ2 / kl + μs.μ.c )
 - ( σ2 / kl + μμr. ) - ( σ2 / kl + μμs. ) - ( σ2 / kl + μμ.c ) - ( σ2 / kl + μμ.d ) } / σ2
=σ2 / kl + ( μrcμsd + μr.μs. + μ.cμ.d + μ2 - μr.μsd - μs.μrc - μ.cμsd - μ.dμrc
 + μrcμ + μsdμ + μr.μ.d + μs.μ.c - μμr. - μμs. - μμ.c - μμ.d ) / σ2
=σ2 / kl + ( μrc - μr. - μ.c + μ )( μsd - μs. - μ.d + μ ) / σ2
=1 / kl

という結果が得られます。しかし、これは r ≠ s かつ c ≠ d の場合であり、r = s かつ c ≠ d ならば

E[mr.ms.] = E[mr.2] = σ2 / k + μr.2

E[mr.xsd] = E[mr.xrd] = σ2 / k + μrdμr.

E[ms.xrc] = E[mr.xrc] = σ2 / k + μrcμr.

なので、

E[yk(r-1)+cyrk+d]={ μrcμrd + ( σ2 / k + μr.2 ) + μ.cμ.d + ( σ2 / kl + μ2 )
 - ( σ2 / k + μrdμr. ) - ( σ2 / k + μr.μrc ) - μ.cμrd - μ.dμrc
 + ( σ2 / kl + μrcμ ) + ( σ2 / kl + μrdμ ) + ( σ2 / kl + μr.μ.d ) + ( σ2 / kl + μr.μ.c )
 - ( σ2 / kl + μμr. ) - ( σ2 / kl + μμr. ) - ( σ2 / kl + μμ.c ) - ( σ2 / kl + μμ.d ) } / σ2
=σ2( 1 / kl - 1 / k ) + ( μrcμrd + μr.2 + μ.cμ.d + μ2 - μr.μrd - μr.μrc - μ.cμrd - μ.dμrc
 + μrcμ + μrdμ + μr.μ.d + μr.μ.c - μμr. - μμr. - μμ.c - μμ.d ) / σ2
=σ2( 1 - l ) / kl + ( μrc - μr. - μ.c + μ )( μrd - μr. - μ.d + μ ) / σ2
=-( l - 1 ) / kl

であり、r ≠ s かつ c = d も同様な計算で

E[yk(r-1)+cysk+c] = -( k - 1 ) / kl

になります。そこで、r 行 c 列目の要素 arc を、

arc=( k - 1 )( l - 1 ) / kl=M( r = c )
=-( l - 1 ) / kl=L( r ≠ c かつ | r - c | < k )
=-( k - 1 ) / kl=K( r ≠ c かつ ( r - c ) ≡ 0 ( mod k ) )
=1 / kl=I( それ以外 )

とする行列 A を考えると、これは yk(r-1)+c を確率変数とする kl 変数の多変量正規分布の共分散行列を表します。どんな行列になるのか上記の条件ではイメージしにくいかもしれませんが、以下のように考えると理解しやすいと思います。
A は kl 行 kl 列の対称行列であり、対角成分は r = c なので全て M になります。全体を k 行 k 列の行列 (これを Auv ( u, v = 1, 2, ... l ) で表します) に区切ったとき、対角成分にあたる行列 Auu は、対角成分は M であり、その他の成分は | r - c | < k となる要素であることから全て L です。
( r - c ) ≡ 0 ( mod k ) にあたる要素 (これは合同式による表現で、r - c が k で割り切れることを表しています) とは、例えば r = 1 ならば c = 1, k + 1, 2k + 1, ... k( l - 1 ) + 1 であり、r ≠ c という条件を加えれば ( r, c ) = ( 1, 1 ) は除外されるので、これは Auv ( u ≠ v ) における 1 行 1 列目の要素を表します。これを一般化すると、r 行目に対してこれを k で割った余り r' に対し、c = k + r', ... k( l - 1 ) + r' 列目の要素、すなわち Auv ( u ≠ v ) における対角成分を意味し、これらが K であることになります。それ以外は I であることから、

Auv=|M,L,...L|( u = v )
|L,M,...L|
|::...:|
|LL...M|
=|K,I,...I|( u ≠ v )
|I,K,...I|
|::...:|
|II...K|

と表され、u = v の場合を Ad、それ以外を And とすれば、

A=|Ad,And,...And|
|And,Ad,...And|
|::...:|
|And,And,...Ad|

になります。一般に、ある行列から任意の行・列を削除したものを「部分行列(Submatrix)」といい、上記のように行列を垂直・水平に区切って分割した部分行列を「ブロック(Block)」といいます。また、ブロックに分割し、表現される行列は「分割行列(Partitioned Matrix)」と呼ばれます。

A2 の r 行 c 列目の要素を a2rc とすると、r = c のとき (つまり、A2 の対角成分の解) は

a2rr=M2 + ( k - 1 )L2 + ( l - 1 )K2 + { kl - ( k - 1 ) - ( l - 1 ) - 1 }I2
={ ( k - 1 )( l - 1 ) / kl }2 + ( k - 1 )( l - 1 )2 / (kl)2
 + ( l - 1 )( k - 1 )2 / (kl)2 + { kl - ( k - 1 ) - ( l - 1 ) - 1 } / (kl)2
={ ( k - 1 )( l - 1 ) / kl }2 + ( k - 1 )( l - 1 )2 / (kl)2
 + ( l - 1 )( k - 1 )2 / (kl)2 + ( k - 1 )( l - 1 ) / (kl)2
={ ( k - 1 )( l - 1 ) / (kl)2 }{ ( k - 1 )( l - 1 ) + ( l - 1 ) + ( k - 1 ) + 1 }
={ ( k - 1 )( l - 1 ) / (kl)2 }kl
=( k - 1 )( l - 1 ) / kl = M ( = arr )

また、r ≠ c かつ | r - c | < k となる場合 (つまり、A2 に対し Ad と等しい行・列のブロックの中で対角成分以外の要素) は、

a2rc=2ML + ( k - 2 )L2 + 2( l - 1 )KI + ( k - 2 )( l - 1 )I2
=-2( k - 1 )( l - 1 )2 / (kl)2 + ( k - 2 )( l - 1 )2 / (kl)2
 - 2( l - 1 )( k - 1 ) / (kl)2 + ( k - 2 )( l - 1 ) / (kl)2
={ ( l - 1 ) / (kl)2 }{ -2( k - 1 )( l - 1 ) + ( k - 2 )( l - 1 ) - 2( k - 1 ) + ( k - 2 ) }
={ ( l - 1 ) / (kl)2 }{ -2( kl - k - l + 1 ) + ( kl - k - 2l + 2 ) - 2( k - 1 ) + ( k - 2 ) }
={ ( l - 1 ) / (kl)2 }(-kl)
=-( l - 1 ) / kl = L ( = arc )

r ≠ c かつ ( r - c ) ≡ 0 ( mod k ) となる場合 (つまり、A2 に対し And と等しい行・列のブロックの中の対角要素) は、

a2rc=2MK + ( l - 2 )K2 + 2( k - 1 )LI + ( k - 1 )( l - 2 )I2
=-2( k - 1 )2( l - 1 ) / (kl)2 + ( l - 2 )( k - 1 )2 / (kl)2
 - 2( k - 1 )( l - 1 ) / (kl)2 + ( k - 1 )( l - 2 ) / (kl)2
={ ( k - 1 ) / (kl)2 }{ -2( k - 1 )( l - 1 ) + ( l - 2 )( k - 1 ) - 2( l - 1 ) + ( l - 2 ) }
={ ( k - 1 ) / (kl)2 }{ -2( kl - k - l + 1 ) + ( kl - l - 2k + 2 ) - 2( l - 1 ) + ( l - 2 ) }
={ ( k - 1 ) / (kl)2 }(-kl)
=-( k - 1 ) / kl = K ( = arc )

その他の場合、

a2rc=2MI + 2KL + 2( k - 2 )LI + 2( l - 2 )KI + ( k - 2 )( l - 2 )I2
=2( k - 1 )( l - 1 ) / (kl)2 + 2( k - 1 )( l - 1 ) / (kl)2 - 2( k - 2 )( l - 1 ) / (kl)2 - 2( l - 2 )( k - 1 ) / (kl)2 + ( k - 2 )( l - 2 ) / (kl)2
={ 4( kl - k - l + 1 ) - 2( kl - k - 2l + 2 ) - 2( kl - 2k - l + 2 ) + ( kl - 2k - 2l + 4 ) } / (kl)2
=kl / (kl)2
=1 / kl = I ( = arc )

となるので、A は「べき等行列(Idempotent Matrix)」です。A は、その固有ベクトルからなる直交行列 Q によって QDQT の形に表されるのでした。ここで、D は固有値を対角成分とする対角行列です(「固有値分解 - (2) カルーネン・レーベ展開」の「1) 対称行列と二次形式」参照 )。べき等行列ということは、A2 = A であることを意味するので、

A2 = QDQTQDQT = QD2QT = QDQT ( = A )

より D2 = D であり、これが成り立つためには固有値すなわち分散が 0 か 1 のいずれかでなければなりません。べき等行列の階数が対角成分の和に等しい(「(12) 二標本の解析 -2-」の「補足2) べき等行列の階数」参照) ことから、A の対角成分の和を求めると

Σi{1→kl}( M ) = kl・( k - 1 )( l - 1 ) / kl = ( k - 1 )( l - 1 )

であり、これは固有値が 1 の対角成分が ( k - 1 )( l - 1 ) 個あって、残りは全てゼロであることを意味します。直交変換によって、確率変数 yk(r-1)+c は ( k - 1 )( l - 1 ) 個の独立な変数 zi ( i = 1, 2, ... (k-1)(l-1) ) に変換され、残り kl - ( k - 1 )( l - 1 ) 個の変数は消失します。直交変換によってノルムは変わらないので、Σr{1→l}( Σc{1→k}( yk(r-1)+c2 ) ) = Σi{1→(k-1)(l-1)}( zi2 ) であり、これは自由度 ( k -1 )( l - 1 ) の χ2-分布に従うことになります。

もし、全ての r, c に対して αr = 0 かつ βc = 0 が成り立てば、ST / σ2 = Σr{1→l}( Σc{1→k}( ( xrc - m )2 ) ) / σ2 は明らかに自由度 kl - 1 の χ2-分布に従います。また、同じ条件下では SR / σ2 = kΣr{1→l}( ( mr. - m )2 ) / σ2 は自由度 l - 1 の、SC / σ2 = lΣc{1→l}( ( m.c - m )2 ) / σ2 は自由度 k - 1 の χ2-分布にそれぞれ従うのでした。ST = SR + SC + SE より、χ2-分布の再生性から ST の自由度は SR, SC, SE それぞれの自由度の和と等しくなることから、SE の自由度を X とすると

kl - 1 = ( l - 1 ) + ( k - 1 ) + X

となるので、X = kl - l - k + 1 = ( k - 1 )( l - 1 ) となり、先ほど得られた結果と一致します。行間・列間変動が大きくなれば、ST, SR, SC は分布からかけ離れた値を取るようになりますが、SE はその影響を受けません。


補足3) 順位和の平方和

順位和 Ti は、次のように正規化することで近似的に標準正規分布に従うのでした。

zi = { Ti - ni( N + 1 ) / 2 } / { niui( N + 1 ) / 12 }1/2

Ti の期待値 E[Ti] は ni( N + 1 ) / 2、分散 V[Ti] は niui( N + 1 ) / 12 なので上記結果が得られるわけですが、Ti と Tj (i≠j) の共分散 γij = E[( Ti - E[Ti] )( Tj - E[Tj] )] はどうなるかというと、

γij=E[( Ti - E[Ti] )( Tj - E[Tj] )]
=E[TiTj] - E[Ti]E[Tj]
=E[TiTj] - ninj( N + 1 )2 / 4

で求めることができます。全体の順位和を T ( = N( N + 1 ) / 2 ) として、Ti がある定数だったときの Tj の値を考えると、それは Ti に含まれる要素を除いた N - ni 個の要素で、その和が T - Ti となるものの中から、nj 個の要素を抽出したときの和となります。N - ni 個の要素から nj 個の要素を抽出する組み合わせは N-niCnj 通りで、各組み合わせには nj 個の要素があります。選択可能な要素は N - ni 個で、全ての組み合わせの中には njN-niCnj 個の要素があり、しかも各要素は等しい回数だけ選ばれているので、その回数は

njN-niCnj / ( N - ni )

で求められ、各要素の和が T - Ti であることから総計は

njN-niCnj( T - Ti ) / ( N - ni )

になります。組み合わせの数は N-niCnj ですから、順位和の平均は

nj( T - Ti ) / ( N - ni )

ということになります。例として、5 個の順位要素から { 1, 3 } が Ti に含まれるとして、{ 2, 4, 5 } の中の二つの要素で Tj を求めることを考えてみます。このとき、N = 5, ni = nj = 2, T = 15, Ti = 4 となります。

{ 2, 4, 5 } から二つの要素を抽出する組み合わせは

{ 2, 4 } { 2, 5 } { 4, 5 }

の 3 通りです。各組み合わせには二個の要素があるので、要素の総数は 2 x 3 = 6 個で、選択可能な要素数は 3 個なので、各要素はそれぞれ 6 / 3 = 2 回ずつ組み合わせの中に出現していることになります。各要素の和は 15 - 4 = 11 ( = 2 + 4 + 5 ) ですから、総計は 11 x 2 = 22 となり、順位和の平均は 22 / 3 と求めることができます。

Ti を固定したとき、Tj の期待値は nj( T - Ti ) / ( N - ni ) となるので、E[TiTj] は

E[TiTj]=E[Tinj( T - Ti ) / ( N - ni )]
={ nj / ( N - ni ) }E[Ti( T - Ti )]
={ nj / ( N - ni ) }{ N( N + 1 )E[Ti] / 2 - E[Ti2] }
={ nj / ( N - ni ) }{ niN( N + 1 )2 / 4 - V[Ti] - E[Ti]2 }
={ nj / ( N - ni ) }{ niN( N + 1 )2 / 4 - ni( N - ni )( N + 1 ) / 12 - ni2( N + 1 )2 / 4 }
={ nj / ( N - ni ) }{ ni( N - ni )( N + 1 )2 / 4 - ni( N - ni )( N + 1 ) / 12 }
=ninj{ ( N + 1 )2 / 4 - ( N + 1 ) / 12 }

になります。従って、

γij=ninj{ ( N + 1 )2 / 4 - ( N + 1 ) / 12 } - ninj( N + 1 )2 / 4
=-ninj( N + 1 ) / 12

であり、相関係数 ρij は ( V[Ti]V[Tj] )1/2 で割ることで得られて

ρij = -( ninj / uiuj )1/2

となります。群が二つしかなければ、N = ni + nj より ui = nj, uj = ni なので、ρij = -1 になります。この結果は、片側の順位和が決まればもう片側の順位和も一意な値になるので、分布は直線上に並ぶことからも容易に理解できます。また、N が ni, nj に対して大きくなるほど相関係数は小さくなることも、この結果からわかります。


この結果から、T = ( T1, T2, ... Tk ) は分散共分散行列 VT

VT= { ( N + 1 ) / 12 }|n1u1,-n1n2,...-n1ni,...n1nk|
|-n2n1,n2u2,...-n2ni,...-n2nk|
|::::|
|-nin1,-nin2,...niui,...-nink|
|::::|
|-nkn1,-nkn2,...-nkni,...nkuk|

としたとき、多変量正規分布

N( T, VT ) = { 1 / ( 2π )k/2|VT|1/2 } exp( -( T - μ, VT-1( T - μ ) ) )

但し、μ = ( n1( N + 1 ) / 2, n2( N + 1 ) / 2, ... nk( N + 1 ) / 2 )

に従い、z = ( z1, z2, ... zk ) は分散共分散行列 VZ

VZ=|1,-( n1n2 / u1u2 )1/2,...-( n1ni / u1ui )1/2,...-( n1nk / u1uk )1/2|
|-( n2n1 / u2u1 )1/2,1,...-( n2ni / u2ui )1/2,...-( n2nk / u2uk )1/2|
|::...:...:|
|-( nin1 / uiu1 )1/2,-( nin2 / uiu2 )1/2,...1,...-( nink / uiuk )1/2|
|::...:...:|
|-( nkn1 / uku1 )1/2,-( nkn2 / uku2 )1/2,...-( nkni / ukui )1/2,...1|

としたとき、多変量正規分布

N( z, VZ ) = { 1 / ( 2π )k/2|VZ|1/2 } exp( -( z, VZ-1z ) )

に従うことになります。VZ が「べき等行列(Idempotent Matrix)」であれば、補足2で示したように固有値が 1 か 0 のいずれかになり、階数から固有値 1 の個数がわかるので、二乗和が階数を自由度とする χ2-分布に従うことが証明できます。しかし、VZ はべき等行列ではありません。そこで、次のような行列の積に分解します。

VZ=|(N/u1)1/2,0,...0,...0||u1/N,-(n1n2)1/2/N,...-(n1ni)1/2/N,...-(n1nk)1/2/N||(N/u1)1/2,0,...0,...0|
|0,(N/u2)1/2,...0,...0||-(n2n1)1/2/N,u2/N,...-(n2ni)1/2/N,...-(n2nk)1/2/N||0,(N/u2)1/2,...0,...0|
|::::||::...:...:||::::|
|0,0,...(N/ui)1/2,...0||-(nin1)1/2/N,-(nin2)1/2/N,...ui/N,...-(nink)1/2/N||0,0,...(N/ui)1/2,...0|
|::::||::...:...:||::::|
|0,0,...0,...(N/uk)1/2||-(nkn1)1/2/N,-(nkn2)1/2/N,...-(nkni)1/2/N,...uk/N||0,0,...0,...(N/uk)1/2|

中央の行列を A としたとき、A のべき乗 A2 の i 行 i 列の対角要素は

nin1 / N2 + nin2 / N2 + ... + ui2 / N2 + ... + nink / N2
={ niΣj{1→k;j≠i}( nj ) + ( N - ni )2 } / N2
={ niΣj{1→k;j≠i}( nj ) + N2 - 2Nni + ni2 } / N2
={ niΣj{1→k}( nj ) + N2 - 2Nni } / N2
=( Nni + N2 - 2Nni ) / N2
=( N - ni ) / N = ui / N

となり、r 行 c 列 ( r ≠ c ) の要素は

( nrn1・n1nc )1/2 / N2 + ( nrn2・n2nc )1/2 / N2 + ... - ur( nrnc )1/2 / N2 + ...
 - ( nrnc )1/2uc / N2 + ... + ( nrnk・nknc )1/2 / N2
{ ( nrnc )1/2Σi{1→k;i≠r,i≠c}( ni ) - ( N - nr )( nrnc )1/2 - ( N - nc )( nrnc )1/2 } / N2
{ ( nrnc )1/2Σi{1→k}( ni ) - 2N( nrnc )1/2 } / N2
-( nrnc )1/2 / N

なので A2 = A が成り立ち、A はべき等行列になります。べき等行列の階数は対角要素の和と等しいことから、

Σi{1→k}( ui / N )=Σi{1→k}( ( N - ni ) / N )
=( kN - N ) / N = k - 1

と求められるので、

yi = ( N / ui )1/2zi

としたとき、yi が標準正規分布に従えば、yi の二乗和は自由度 k - 1 の χ2-分布に従うことになります。ところが、標準正規分布に従うのは zi の方だったので、yi が標準正規分布に従うように zi に補正を掛ける必要があります。そこで、zi に ( ui / N )1/2 を掛ければ、yi は標準正規分布に従い、yi の二乗和、すなわち ( ui / N )1/2zi の二乗和は自由度 k - 1 の χ2-分布に従うことになります。この値は、H = SC / σ2 そのものです。


補足4) SG / σ2 が自由度 k - 1 の χ2-分布に従うことの証明

フリードマン検定において、SG / σ2 が自由度 k - 1 の χ2-分布に従うことを証明するには、前に示した「補足4) 順位和の平方和」の内容を流用すれば簡単に証明できます。まず、要素数 k である群の中の要素 Rgi に対する期待値 E[Rgi] と分散 V[Rgi] は

E[Rgi] = ( k + 1 ) / 2

V[Rgi] = ( k + 1 )( k - 1 ) / 12

となるのでした。また、同じ群の中の二つの要素 Rgi, Rgi の共分散は、「補足4) 順位和の平方和」で求めた順位和の共分散 γij = -ninj( N + 1 ) / 12 において、ni = nj = 1, N = k とした場合に等しいので

γij = -( k + 1 ) / 12

であり、相関係数 ρij

ρij = -1 / ( k - 1 )

になります。群が n 個あって、それぞれ独立であるとしたとき、各群の i 番目の要素の和を Ti と表せば、この期待値 E[Ti] と分散 V[Ti] は各群の値の和となって

E[Ti] = n( k + 1 ) / 2

V[Ti] = n( k + 1 )( k - 1 ) / 12

という結果が得られます。また、E[TiTj] は

E[TiTj]=E[Σg{1→n}( Rgig{1→n}( Rgj )]
=Σg{1→n}( Σh{1→n;h≠g}( E[Rgi]E[Rhj] ) ) + Σg{1→n}( E[RgiRgj] )
=Σg{1→n}( Σh{1→n;h≠g}( E[Rgi]E[Rhj] ) ) + Σg{1→n}( γij + E[Rgi]E[Rgj] )
=Σg{1→n}( Σh{1→n}( E[Rgi]E[Rhj] ) ) + nγij
=E[Ti]E[Tj] + nγij

なので、Ti, Tj の共分散と相関係数を γTij, ρTij とすれば、

γTij = nγij = -n( k + 1 ) / 12

ρTij = -1 / ( k - 1 )

となります。よって、

zi = { Ti - n( k + 1 ) / 2 } / { n( k + 1 )( k - 1 ) / 12 }1/2

としたとき、z = ( z1, z2, ... zk ) は分散共分散行列 VZ

VZ=|1,-1 / ( k - 1 ),...-1 / ( k - 1 ),...-1 / ( k - 1 )|
|-1 / ( k - 1 ),1,...-1 / ( k - 1 ),...-1 / ( k - 1 )|
|::...:...:|
|-1 / ( k - 1 ),-1 / ( k - 1 ),...1,...-1 / ( k - 1 )|
|::...:...:|
|-1 / ( k - 1 ),-1 / ( k - 1 ),...-1 / ( k - 1 ),...1|

とする k 変量正規分布

N( z, VZ ) = { 1 / ( 2π )k/2|VZ|1/2 } exp( -( z, VZ-1z ) )

に従うことになります。ここで、VZ を次のような行列の積に分解します。

VZ=|{k/(k-1)}1/2,0,...0,...0||(k-1)/k,-1/k,...-1/k,...-1/k||{k/(k-1)}1/2,0,...0,...0|
|0,{k/(k-1)}1/2,...0,...0||-1/k,(k-1)/k,...-1/k,...-1/k||0,{k/(k-1)}1/2,...0,...0|
|::::||::...:...:||::::|
|0,0,...{k/(k-1)}1/2,...0||-1/k,-1/k,...(k-1)/k,...-1/k||0,0,...{k/(k-1)}1/2,...0|
|::::||::...:...:||::::|
|0,0,...0,...{k/(k-1)}1/2||-1/k,-1/k,...-1/k,...(k-1)/k||0,0,...0,...{k/(k-1)}1/2|

中央の行列を A とすると、そのべき乗行列 A2 の対角要素 a2ii

a2ii={ ( k - 1 ) / k }2 + ( k - 1 ) / k2
={ ( k2 - 2k + 1 ) + ( k - 1 ) } / k2
=k( k - 1 ) / k2
=( k - 1 ) / k

であり、非対角要素 a2rc ( r ≠ c ) は

a2rc=-2( k - 1 ) / k2 + ( k - 2 ) / k2
={ ( -2k + 2 ) + ( k - 2 ) } / k2
=-k / k2
=-1 / k

なので、A はべき等行列ということになります。べき等行列の階数は対角要素の和と等しいことから、

Σi{1→k}( ( k - 1 ) / k ) = k - 1

と求められるので、

yi = ( k / ( k - 1 ) )1/2zi

としたとき、yi が標準正規分布に従えば、yi の二乗和は自由度 k - 1 の χ2-分布に従うことになります。ところが、標準正規分布に従うのは zi の方だったので、yi が標準正規分布に従うように zi に補正を掛ける必要があります。そこで、zi に ( ( k - 1 ) / k )1/2 を掛ければ、yi は標準正規分布に従い、yi の二乗和、すなわち ( k / ( k - 1 ) )1/2zi の二乗和は自由度 k - 1 の χ2-分布に従うことになります。この値は、SG / σ2 そのものです。


<参考文献>
1. 「統計数学入門」 本間鶴千代著 (森北出版)
2. 「統計のための行列代数(上)」 D.A.ハーヴィル著 (シュプリンガー・ジャパン)
3. 「Vassar College」-「Concepts & Applications of Inferential Statistics
3-1. The Kruskal-Wallis Test
3-2. The Friedman Test
4. 我楽多頓陳館 --- 統計学入門
5. Wikipedia

◆◇◆更新履歴◆◇◆

◎ 「5) フリードマン検定(Friedman Test)」のサンプル・プログラムを修正しました (2015-08-02)
計算時の一時変数の型に誤りがあり、正しいカイ二乗値が計算できていませんでした。出力結果も併せて修正してあります。
◎ 「3) 反復のある二元配置分散分析法(Two-way Repeated-Measures ANOVA)」のサンプル・プログラムを修正しました (2015-12-13)
単純なコーディング・ミスで交互作用に対する p 値が正しく計算されていませんでいた。出力結果も併せて修正してあります。

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system