分散分析の基礎

この記事では、分散分析と呼ばれる検定手法について解説します。特に一元配置分散分析の解説となります。
分散分析の理論と、ソフトウェアを使った実行方法を解説します。

まず分散分析の基本事項を整理してから、分散分析で用いられる検定統計量であるF比について解説します。F比の計算、F比の解釈を通して「分散分析が何を目指しているのか」を理解してほしいです。
続いて統計モデルやパラメトリックブートストラップ検定の解説を通して、p値について解説します。

RとPythonでの実装方法についても解説します。
データやコードはGitHubで公開しています。

この記事では、統計学の初歩をすでに勉強したことがあるという人を対象にしています。平均や分散、あるいは統計的仮説検定という言葉を初めて聞いたという方は、当ブログの入門記事「統計学基礎」などもあわせて参照してください。

 


スポンサードリンク

 

目次

1. 分散分析の使い道
2. 分散分析の学び方
3. 分散分析の検定統計量（F比）
   1. 差があるかどうかを判断することの難しさ
   2. F比の直観的な説明
   3. F比を計算する手順
   4. 補足：F比を計算する手順を数式で表現する
4. 統計モデルとp値
   1. 統計モデルと確率分布
   2. 分散分析モデル
   3. 正規分布の初歩
   4. 検定統計量の標本分布とp値
5. 分散分析の実行
   1. Rによる実行
   2. Pythonによる実行
6. パラメトリックブートストラップ検定（PB検定）
   1. PB検定の考え方
   2. 帰無仮説が正しいことを想定した分散分析モデル
   3. 乱数の生成
   4. PB検定の実行

 

1．分散分析の使い道

分散分析の仕組みに移る前に、分散分析の使い道を紹介します。

分散分析は、多群において平均値に差があるかどうかを判断する時に使います。多群というのは2群でもよいし、3群以上でも良いです。
例えば「東京都と大阪府（2群）で、小学6年生の平均身長に差があるか」や「群馬県と山梨県と茨城県（3群）で、40代男性の体重に差があるか」などを判定する際に使われます。
なお、差があるかを判断したい興味の対象を因子と呼びます。因子のカテゴリーを水準と呼びます。「東京都と大阪府（2群）で、小学6年生の平均身長に差があるか」という問題の因子は場所で、水準は東京都と大阪府です。「群馬県と山梨県と茨城県（3群）で、40代男性の体重に差があるか」という問題の水準は群馬県と山梨県と茨城県です。

平均値に差があるかどうかを知りたい場面はしばしばあります。生物の研究をしている人でも「食べ物をA,B,Cに変えた時、家畜の体重の平均値に差があるか」を調べるならば分散分析が登場するでしょう。利用場面は広いです。

当サイトでは以前にt検定を取り上げました。t検定も平均値に差があるかどうかを判断する際に使用されます。しかしt検定は2つの水準での平均値の差を調べることしかできません。3つ以上の水準を相手にする場合は、例えば「AとBの差の検定」「AとCの差の検定」「BとCの差の検定」と3回繰り返す必要が出てきます。詳しくは説明しませんが、これは検定の多重性の問題が発生して好ましくありません。
水準の数が3つ以上の場合は、分散分析を使用することを検討しましょう。

 

2．分散分析の学び方

分散分析を学ぶ際には、以下の手順を踏むのをお勧めします。

1. 検定統計量について学ぶ
   1. 検定統計量の解釈を学ぶ
   2. 検定統計量の計算方法を学ぶ
2. p値の計算方法について学ぶ
   1. 統計モデルについて学ぶ
   2. 正規分布について学ぶ
   3. 統計検定量の標本分布の考え方を学ぶ
   4. p値を計算する方法とp値の解釈を学ぶ

ざっくりいうと、統計的仮説検定は、以下の手順で計算を進めます。

1. 統計モデルを作る
2. 検定統計量を計算する
3. 検定統計量からp値を計算する
4. p値の結果からなにがしかの判断を下す(あるいは判断を保留する)

最後の「p値に基づき判断を下す」ところだけピックアップすると、計算の意味がとらえにくくなります。
前半の統計モデル、そして検定統計量が何者か理解できれば、統計的仮説検定の理解度は大きく増すはずです。

 

3．分散分析の検定統計量（F比）

ここから本格的に分散分析の解説をします。まずは分散分析における検定統計量であるF比の解説をします。

 

差があるかどうかを判断することの難しさ

今回は「食べ物をA,B,Cに変えた時、家畜の体重の平均値に差があるか」を調べる目的で分散分析を使うことにします。因子は餌の種類であり、水準は餌A、餌B、餌Cの3つです。
なお、ここでの分散分析は一元配置分散分析と呼ばれます。

ところで、分散分析などという面倒な手続きをなぜする必要があるのでしょうか。例えば「家畜の体重が5㎏以上異なっていれば、差があるとみなす」というような単純なルールを設定すれば良いように思えます。けれども、このやり方ではうまくいきません。

例えば以下のデータがあったとします。
結果①
餌A：平均体重50㎏
餌B：平均体重51㎏
餌C：平均体重49㎏

このとき「重さが1㎏ほどしか離れていないから差がない」と言いたくなります。

同じテーマで以下の結果が得られたとします。
結果②
餌A：平均体重2㎏
餌B：平均体重3㎏
餌C：平均体重1㎏

元の家畜の体重が軽いので、1㎏の差が重要な意味を持つように思えます。単なる平均値の差の大きさで判断を下すのは問題がありそうです。

分散分析では検定統計量を計算し、p値を計算し……と少し面倒な手続きを踏みます。ぱっと見で差があるかどうかわからない場合には、いろいろな観点から評価をして判断を下したいですね。判断をするための一助として、少し面倒ではありますが、分散分析という検定の手続きを行うわけです。
なお、仮説検定はあくまでも「判断をサポートする」ためのものです。仮説検定の結果だけを見て機械的に判断を下すことはできないので注意しましょう。最終的には、検定結果に加えて、グラフや生物学的な解釈など、いろいろな観点から総合的に判断を下すことになります。

 

F比の直観的な説明

検定統計量という言葉は少し難しいですが、分散分析の場合「検定統計量が大きな値を取れば、平均値に差があるとみなせるだろう」と考えられるような値だと思うと良いでしょう。後ほど登場するp値とも密接な関わりがあります。分散分析の検定統計量はF比と呼ばれます。
分散分析の目的は、水準ごとに平均値に差があるかどうかを判断することでした。検定統計量であるF比を計算して、この値が大きいかどうかを判断することで、平均値に差があるかどうかを判断する流れです。

分散分析は、名前の通り分散が大きな役割を果たします。F比は以下のように計算されます。

$$F比=\frac{水準間のばらつき}{水準内のばらつき}$$

定義通りですと理解しにくいかもしれません。F比のイメージは以下のようなところになります。

$$F比=\frac{効果の大きさ}{誤差の大きさ}$$

効果というのは、今回の事例では餌の違い（異なる水準の間にある違い）となります。餌が変わることによって体重が大きく変わった（異なる水準の間で体重が大きく変わった）ならば「効果が大きい」と言えます。

餌によって体重が変わるかもしれません。けれども生き物には個性があります。同じ餌を食べている個体（同じ水準内にいる個体）であっても、ある個体はたまたま体重が重めで、ある個体はたまたま体重が軽いということもあります。これが「誤差の大きさ」です。誤差が大きいと「餌の違いによって体重が変わったように見えるが、これはたまたまで、実際は生き物の個体差にすぎない」ということになるかもしれません。

そこでF比の登場です。
F比では効果の大きさ(餌の違いがもたらす影響)と、誤差の大きさ(生き物の個体差など)の比を取ります。そして「誤差があることを加味しても、それでもなお、大きな効果がある」とみなせるかどうかを評価します。
F比の大きさは、餌によって体重が変わるかどうかを判断するための根拠の1つになるはずです。

 

F比を計算する手順

実際にF比を計算してみます。数値例は馬場(2015)『平均・分散から始める一般化線形モデル入門』と同じものを使います。

データ
餌A：2kg　　4kg
餌B：10kg　12kg
餌C：6kg　　8kg

 

効果の大きさの評価のイメージ

まずはF比の分子である「効果の大きさ」を評価します。
餌の影響を無視して、6個体すべての体重の平均値を求めます。これを総平均\( \bar{x} \)と呼ぶことにします。結果は7になります。

$$ \bar{x} = \frac{2+4+10+12+6+8}{6} = 7 $$

続いて、餌の種類ごとに体重の平均値を求めます。各々を\( \bar{x}_1,\bar{x}_2,\bar{x}_3 \)と表記することにします。
餌A：\( \bar{x}_1 = \frac{2+4}{2}=3 \)
餌B：\( \bar{x}_2 = \frac{10+12}{2}=11 \)
餌C：\( \bar{x}_3 = \frac{6+8}{2}=7 \)

ここで「餌の種類を無視した総平均(7)」と、「個別の餌ごとの体重の平均値(3，11，7)」が大きく異なっていれば、餌には体重を変化させる効果があったと言えそうです。これは重要なポイントなので覚えておいてください。

 

元データの平方和と不偏分散

先述の通り、分散分析では、ばらつき、すなわち分散を使って誤差や効果の大きさを表現します。まずは分散や平方和といった用語の復習から進めます。
不偏分散の計算式は以下の通りです。ただし\( N \)はサンプルサイズ、\( x_i \)は個別のデータ、\( \mu \)は総平均です(この計算式の意味が分からなければ、記述統計の基礎を参照してください)。伝統的に分散は\( \sigma^2 \)と表記します。

$$ \sigma^2 = \frac{1}{N-1} \sum_{i=1}^{N} (x_i – \mu)^2 $$

ここで\( N-1 \)で割る前の以下の結果を偏差平方和と呼びます。Sum of Squaredの略でSSと略します

$$ SS = \sum_{i=1}^{N} (x_i – \mu)^2 $$

元データに対して不偏分散を計算してみます。総平均\( \mu = 7 \)であることに注意してください。
偏差平方和SSは以下の通りです。

$$ \begin{eqnarray}
SS &=& (2-7)^2 + (4-7)^2 + (10-7)^2 + (12-7)^2 + (6-7)^2 + (8-7)^2 \\
&=& 70
\end{eqnarray}$$

偏差平方和の計算結果を使って、不偏分散を求めます。

$$ \sigma^2 = \frac{70}{6-1} = 14 $$

 

効果の大きさ(水準間の分散)

効果の大きさを評価する際は、餌の種類ごとの平均値を使います。すなわち、以下のデータが得られたと想定してばらつきを計算します。
\( \bar{x}_1=3, \bar{x}_2=11, \bar{x}_3=7 \)であることを思い出してください。
餌A：3kg　　3kg
餌B：11 kg　11kg
餌C：7kg　　7kg

本来は餌Aの個体は2kgと4㎏でした。しかし、餌の種類がもたらした効果を測定したいので、個体差は無視してしまいます。
餌の種類ごとの平均値そのもののデータが得られたと考えるわけです。
このようにして、個体差を無視した「餌の効果の大きさ」を評価します。
「餌の効果の大きさ」を「餌の種類の違いがもたらす、体重のばらつき」で表現するわけです。

上記の数値を使って偏差平方和を求めます。水準間の偏差平方和をBetweenのBを使って\( SS_B \)と表記します。総平均\( μ=7 \)であることに注意してください。
見やすさのために中カッコを入れています。

$$ \begin{eqnarray}
SS_B &=& \{ (3-7)^2 + (3-7)^2 \} + \{ (11-7)^2 + (11-7)^2 \} + \{ (7-7)^2 + (7-7)^2 \} \\
&=& 64
\end{eqnarray}$$

続いて、水準間のばらつきを求めます。なお、水準間や水準内のばらつきは分散と呼ばずに平均平方と呼ぶことが多いので、本記事でもそれにならいます。
水準間の平均平方を\( \sigma_{B}^{2} \)と表記することにします。上記の計算結果を使って\( SS_B \div (6-1) \)として分散を求めたいところです。
しかし、餌の種類が3種類しかないのに、サンプルサイズ（6）を使って不偏分散を計算するのには違和感があります。
そのため、餌の種類（3種類）を使って以下のように「3-1」で割ります。分母の値は自由度とも呼ばれます。これを\( \nu_B \)と表記することにします。\( \nu_B = 3-1 = 2 \)です。

$$ \sigma_{B}^{2} = \frac{SS_B}{\nu_B} = \frac{64}{3-1} = 32 $$

 

誤差の大きさ(水準内の分散)

続いて誤差の大きさを評価します。これは「餌の違いを考慮した後に、それでも残ったばらつき」として評価されます。

データは元データをそのまま使います。
餌A：2kg　　4kg
餌B：10 kg　12kg
餌C：6kg　　8kg

ここで偏差平方和を求める際、総平均\( \mu \)ではなく「餌の違いを考慮した平均値」を使うことに注意が必要です。すなわち以下の結果を使います。
餌Aの平均値：\( \bar{x}_1=3 \)
餌Bの平均値：\( \bar{x}_2=11 \)
餌Cの平均値：\( \bar{x}_3=7 \)

上記の数値を使って偏差平方和を求めます。水準内の偏差平方和をWithinのWを使って\( SS_W \)と表記します。見やすさのために中カッコを入れています。

$$ \begin{eqnarray}
SS_W &=& \{(2-3)^2 + (4-3)^2 \} + \{(10-11)^2 + (12-11)^2 \} + \{(6-7)^2 + (8-7)^2 \} \\
&=& 6
\end{eqnarray}$$

続いて、水準内の平均平方を求めます。水準内の平均平方を\( \sigma_{W}^{2} \)と表記することにします。
分母には注意が必要です。餌ごとにサンプルサイズが決まっています。まずは薬Aのサンプルサイズから1を引きます。さらに薬B、Cのサンプルサイズからも1を引きます。その結果の合計値が分母になります。
水準内の平均平方の分母の値（自由度）を\( \nu_W \)と表記することにします。\( \nu_W = (2-1)+(2-1)+(2-1) = 3 \)です。

$$ \sigma_{W}^{2} = \frac{SS_W}{\nu_W} = \frac{6}{(2-1)+(2-1)+(2-1)} = 2 $$

 

F比の計算

最後にF比を求めます。

$$ F比=\frac{\sigma_{B}^{2}}{\sigma_{W}^{2}} = 16 $$

誤差のばらつきよりも、効果のばらつきの方が16倍も大きいようです。
F比が大きければ、誤差よりも効果の方が大きい、すなわち餌の効果によって体重が変化したとみなすことができそうです。
ただしF比が15を超えれば「F比が大きい」とみなせるのか、あるいは20を超えなければ大きいとみなせないのかを判断しなければなりません。
統計的仮説検定ではp値と呼ばれる指標を使ってF比が大きいかどうかを判断します。p値は統計モデルの解説をした後に導入します。

 

補足：平方和の分解

今までの計算では3つの平方和が出てきました。元のデータに対する平方和\( SS \)と、水準間の平方和\( SS_B \)、そして水準内の平方和\( SS_W \)です。

ここで、以下の関係があることが知られています。

$$ SS = SS_B + SS_W $$

これを平方和の分解と呼びます。
今回のデータでは\( SS=70, SS_B=64, SS_W=6 \)でしたので、上記の関係が成り立っているのがわかります。

平方和の分解は「データのもともとのばらつきは、餌の効果がもたらしたばらつきと、誤差がもたらしたばらつき、の2つで表現できる」と解釈できます。
これは後ほど統計モデルとしての分散分析を理解する際に重要となってくるので、頭の片隅に置いておいてください。

 

補足：F比を計算する手順を数式で表現する

統計学の教科書を眺めると、Σ記号がたくさん並んだ複雑な数式をお目にかけるかと思います。計算の流れは先ほど説明した通りなのですが、中級以上の教科書を読み進めるための訓練になることを意図して、ここでは数式を使ったF比の計算手続きを解説します。難しいと感じたら、飛ばし読みしても大丈夫です。

まず、データを以下のように表記します。

$$
\begin{eqnarray}
\begin{array}{|c|c|c|c|c|c|c|}
\hline
& Level_1 & Level_2 & \ldots & Level_j & \ldots & Level_J \\
\hline
1 & x_{ 11 } & x_{ 12 } & \ldots & x_{ 1j } & \ldots & x_{ 1J } \\
\hline
2 & x_{ 21 } & x_{ 22 } & \ldots & x_{ 2j } & \ldots & x_{ 2J } \\
\hline
\vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\
\hline
i & x_{ i1 } & x_{ i2 } & \ldots & x_{ ij } & \ldots & x_{ iJ } \\
\hline
\vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\
\hline
n_j & x_{ n_j1 } & x_{ n_j2 } & \ldots & x_{ n_jj } & \ldots & x_{ n_jJ } \\
\hline
\end{array}
\end{eqnarray}
$$

水準の数が\( J \)であるデータがあるとします。水準間で平均値に差があるかどうかを調べます。
ある水準\( j \)におけるサンプルサイズを\( n_j \)と表記します。\( n_j \)は繰り返し数と呼ぶこともあります。水準によって\( n_j \)が変わっても大丈夫です。なお、\( n_j \)をすべての水準で足し上げると総サンプルサイズ\( N \)になります。すなわち下記が成り立ちます。

$$ N = \sum_{j=1}^{J}n_j $$

ある水準\( j \)における第\( i \)番目のデータを\( x_{ij} \)とします。ただし\( i \leq n_j \)です。
ある水準\( j \)におけるデータの平均値を\( \bar{x}_j \)とします。\( \bar{x}_j\)は以下のように計算されます。

$$ \bar{x}_j = \frac{1}{n_j} \sum_{i=1}^{n_j}x_{ij} $$

すべてのデータにおける総平均を\( \bar{x} \)とします。\( \bar{x} \)は以下のように計算されます

$$ \bar{x}_j = \frac{1}{N} \sum_{j=1}^{J} \sum_{i=1}^{n_j}x_{ij} $$

シグマ記号が2つ並ぶと少し読みづらさを感じるかもしれません。
まずは\( \sum_{i=1}^{n_j}x_{ij} \)で水準ごとに合計値を計算します。そしてそれをすべての水準で足し合わせます。
「水準の中での足し算」と「水準での合算」の2回足し算が行われているのでシグマ記号が2つ並びます。

水準間の平方和\( SS_B \)は以下のように計算されます。

$$ SS_B=\sum_{j=1}^{J}{n_j\left({\bar{x}}_j-\bar{x}\right)}^2 $$

水準内の平方和\( SS_W \)は以下のように計算されます。

$$ SS_W=\sum_{j=1}^{J}\sum_{i=1}^{n_j}\left(x_{ij}-{\bar{x}}_j\right)^2 $$

水準間の平均平方和\( \sigma_{B}^{2} \)は以下のように計算されます。

$$ \sigma_B^2=\frac{SS_B}{J-1} $$

水準内の平均平方和\( \sigma_{W}^{2} \)は以下のように計算されます。

$$ \sigma_W^2=\frac{SS_W}{\sum_{j=1}^{J}\left(n_j-1\right)} $$

 

数値例との対応を見ていきます。
数値例では\( A,B,C \)の3水準あったので\( j=1,2,3 \)であり、\( J=3 \)となります。ただし\( j=1 \)のときを水準\( A \)と、\( j=2 \)のときを水準\( B \)と、\( j=3 \)のときを水準\( C \)と定めます。
数値例では全てサンプルサイズは2ずつだったので、\( n_1=n_2=n_3=2 \)です。すべての水準を合算したサンプルサイズは6です。
総平均\( \bar{x} \)は7です。

\( x_{ij}\)と\( {\bar{x}}_j\)は以下のようになります。

$$
\begin{eqnarray}
\begin{array}{|c|c|c|c|}
\hline
& A(j=1) & B(j=2) & C(j=3)\\
\hline
i=1 & x_{ 11 }=2 & x_{ 12 }=10 & x_{ 13 }=6 \\
\hline
i=2 & x_{ 21 }=4 & x_{ 22 }=12 & x_{ 23 }=8 \\
\hline
平均値 & \bar{x}_1=3 & \bar{x}_2=11 & \bar{x}_3=7 \\
\hline
\end{array}
\end{eqnarray}
$$

\( x_{ij} \)のような記号を使うと、水準の数\( J \)や水準ごとの繰り返し数\( n_j \)が変わっても、計算式が変化しません。
記号が出てくると読みづらく感じるかもしれませんが、慣れれば便利です。

 

4．統計モデルとp値

ここでは統計モデルの初歩を解説します。ここの議論は少し抽象的でわかりにくいかもしれません。私が統計モデルの解説をするときは、ほぼ毎回シミュレーションを使って説明します。具体的に数値を見た方が理解しやすいと思うからです。
最初はどうしても抽象的な話が続きますが、後ほどシミュレーションによる解説をします。まずは用語をざっと眺める程度でいいので読み進めていってください。

 

統計モデルと確率分布

統計モデルは「データが得られるプロセス」として用いられます。Upton,Cook(2010)『統計学辞典』ではモデルを『観測したデータを生み出す確率的な過程を簡潔に記述したもの』と定義しています。

統計モデルは確率分布を使って表現されます。例えば「平均3、分散4の正規分布」などはシンプルな統計モデルです。

統計モデルは「データが得られるプロセス」であり、「確率分布を使って表現されるもの」です。この2つの意味のつながりをおさえるのがとても大切です。

赤いボールを渡されたとします。次は何色のボールが来るでしょうか。
普通はわかりません。でも、「赤いボールだけが入った箱からボールを取り出したのだ」という「データが得られるプロセス」がわかっていれば、次も赤いボールが出てくると予測できます。
「赤いボールが1000個、青いボールが1000個入った箱からランダムにボールを取り出したのだ」という「データが得られるプロセス」がわかっていれば、赤いボールと青いボールが半々の確率で出てくると予測できます。

統計学ではデータが得られるプロセスを確率分布で表現します。
「赤いボールだけが入った箱からボールを取り出したのだ」というデータ生成プロセスであれば、「赤いボールが出る確率が1」と解釈できます。
「赤いボールが1000個、青いボールが1000個入った箱からランダムにボールを取り出したのだ」というデータ生成プロセスであれば、「赤いボールが出る確率が50%、青いボールが出る確率が50%」と解釈されます。

統計学では「確率分布」と呼ばれる言葉がよく出てきます。平たく言うと、確率的に変化する値と、その時の確率をセットにしたものが、確率分布です。
「赤いボールが出る確率が50%、青いボールが出る確率が50%」というのも、1つの確率分布ですね。

確率分布は、データが得られるプロセスとして使用されます。ここはとても大切なことなので、覚えておいてください。

分散分析においても、「分散分析モデル」という統計モデルが使われています。このモデルでは正規分布と呼ばれる確率分布が使われています。

 

分散分析モデル

 

分散分析モデルのモデル式

分散分析モデルは正規分布を使うことで以下のように定式化されます。

$$ x_{ij}=\mu_0+\delta_j+\varepsilon_{ij},\ 　\varepsilon_{ij} \sim N(0,\sigma^2) $$

日本語で書き下すと以下のようなモデル式になります。

$$ データ=総平均+処理効果+誤差 $$

ただし\( N\left(\mu,\sigma^2\right) \)は平均\( \mu \)で分散\( \sigma^2 \)の正規分布という意味です。\( \varepsilon_{ij} \sim N\left(0,\sigma^2\right) \)は\( \varepsilon_{ij} \)が平均0で分散\( \sigma^2 \)の正規分布に従うということを表しています。
\( \mu_0 \)は総平均\( \bar{x} \)です。\( \delta_j \)は処理の効果です。効果\( \delta_j \)には添え字\( j \)がついているので、水準ごとに効果が変わることを想定しているのが分散分析モデルです。

 

数値例

モデル式のままだと少し難しいので数値例を当てはめながら見ていきます。
数値例では\( A,B,C \)の3水準ありました。そのため\( j=1,2,3 \)です。ただし\( j=1 \)のときを水準\( A \)と、\( j=2 \)のときを水準\( B \)と、\( j=3 \)のときを水準\( C \)と定めます。
総平均は\( \mu_0 = \bar{x} = 7 \)でした。水準ごとの平均値は\( {\bar{x}}_1=3,{\bar{x}}_2=11,{\bar{x}}_3=7 \)です。
処理の効果\( \delta_j \)は、各々\( \delta_1=-4,\delta_2=4,\delta_3=0 \)です。\( {\bar{x}}_j=\mu_0+\delta_j \)となっていることに注意してください。
誤差\( \varepsilon_{ij} \)は、モデル式を変形して\( \varepsilon_{ij}=x_{ij}-\mu_0-\delta_j \)で計算されます。

表で整理しておきます。まずはデータです。
$$
\begin{eqnarray}
\begin{array}{|c|c|c|c|}
\hline
& A(j=1) & B(j=2) & C(j=3)\\
\hline
i=1 & x_{ 11 }=2 & x_{ 12 }=10 & x_{ 13 }=6 \\
\hline
i=2 & x_{ 21 }=4 & x_{ 22 }=12 & x_{ 23 }=8 \\
\hline
\end{array}
\end{eqnarray}
$$

効果\( \delta_j \)です。
$$
\begin{eqnarray}
\begin{array}{|c|c|c|c|}
\hline
& A(j=1) & B(j=2) & C(j=3)\\
\hline
i=1 & \delta_1=-4 & \delta_2=4 & \delta_3=0 \\
\hline
i=2 & \delta_1=-4 & \delta_2=4 & \delta_3=0 \\
\hline
\end{array}
\end{eqnarray}
$$

誤差\( \varepsilon_{ij} \)は以下のようになります。例えば\( \varepsilon_{11}=x_{11}-\mu_0-\delta_1=2-7-\left(-4\right)=-1 \)となります。

$$
\begin{eqnarray}
\begin{array}{|c|c|c|c|}
\hline
& A(j=1) & B(j=2) & C(j=3)\\
\hline
i=1 & \varepsilon_{ 11 }=-1 & \varepsilon_{ 12 }=-1 & \varepsilon_{ 13 }=-1 \\
\hline
i=2 & \varepsilon_{ 21 }=1 & \varepsilon_{ 22 }=1 & \varepsilon_{ 23 }=1 \\
\hline
\end{array}
\end{eqnarray}
$$

これで\( x_{ij}=\mu_0+\delta_j+\varepsilon_{ij} \)となりますね。

 

モデル式と今までの計算の比較

モデルは現象の理解を促進するのに役立ちます。下記の式は直観的に理解がしやすい形式になっているのではないかなと思います。

$$ x_{ij}=\mu_0+\delta_j+\varepsilon_{ij},\ 　\varepsilon_{ij} \sim N(0,\sigma^2) $$

$$ データ=総平均+処理効果+誤差 $$

ここで、モデルの構成要素である\( \delta_j \)や\( \varepsilon_{ij} \)を使って、F比の計算方法を再確認しておきます。

まずは効果の大きさとして、水準間の平方和\( SS_B \)を計算します。これはシンプルに効果\( \delta_j \)の表の値を全て2乗して合計すれば求まります。
$$ $$

$$ \begin{eqnarray}
SS_B &=& \{ \left(-4\right)^2+\left(-4\right)^2 \} + \{ \left(4\right)^2+\left(4\right)^2 \} + \{ \left(0\right)^2+\left(0\right)^2 \} \\
&=& 64
\end{eqnarray}$$

続いて誤差の大きさとして、水準内の平方和\( SS_W \)を計算します。これはシンプルに誤差\( \varepsilon_{ij} \)の表の値を全て2乗して合計すれば求まります
$$ \begin{eqnarray}
SS_W &=& \left\{\left(-1\right)^2+\left(1\right)^2\right\}+\left\{\left(-1\right)^2+\left(1\right)^2\right\}+\left\{\left(-1\right)^2+\left(1\right)^2\right\} \\
&=& 6
\end{eqnarray}$$

後はこの値から平均平方を求めるだけですね。以下の計算手順は変わらないので省略します。
モデル式を使うと、データの解釈がしやすくなります。もちろんF比の計算も支障なく行えます。

 

モデル式の別の表現

今までは下記のモデル式を使っていました。

$$ x_{ij}=\mu_0+\delta_j+\varepsilon_{ij},\ 　\varepsilon_{ij} \sim N(0, \sigma^2) $$

多くの教科書ではこの形式になっていると思います。ところで、これは下記のように表記してもまったく同じです。

$$ x_{ij} \sim N\left(\mu_0+\delta_j,\sigma^2\right) $$

データ\( x_{ij} \)は「平均値が\( \mu_0+\delta_j \)で、分散が\( \sigma^2 \)の正規分布」に従って得られると読めます。\( \varepsilon_{ij} \)を消した表記の方がシンプルですね。モデルとは『観測したデータを生み出す確率的な過程を簡潔に記述したもの』です。「データは、ある確率分布から生成されている」というイメージを持ちやすいのは後者の表記だと思います。こちらの表記にも慣れておきましょう。

 

正規分布の初歩

 

正規分布の確率密度関数

正規分布の確率密度関数は以下のようになります。
ただし\( N\left(\mu,\sigma^2\right) \)は平均\( \mu \)で分散\( \sigma^2 \)の正規分布という意味です。
確率変数\( X \)を明示して\( N\left(X\middle|\mu,\sigma^2\right) \)と表記することもあります。

$$ N\left(X\middle|\mu,\sigma^2\right)=\frac{1}{\sqrt[]{2\pi\sigma^2}} e^{\{-\frac{(x-μ)^2}{2\sigma^2}\}} $$

上記の数式を初めて見たという方は「確率密度関数と正規分布」などもあわせて参照してください。
正規分布の確率密度関数を使うと、「平均値(期待値)、分散」が求まれば、確率変数\( X \)の確率密度がすくに計算できますね。

 

正規分布と分散分析モデル

正規分布が重要である理由は「中心極限定理」と呼ばれる定理にあります。
大雑把に言うと、中心極限定理は以下のことを示します。
「分散が無限大ではない確率分布から、単純ランダムサンプリングによって得られた標本の期待値、あるいは合計値は、サンプルサイズを大きくすれば、正規分布に近づく」
例えばサイコロを少し多めに1万回振ったとします。そして、1万回のサイコロの目の合計値を計算します。すると、このサイコロの目の合計値は、正規分布に従います。

正規分布は誤差分布とも呼ばれます。たとえば平均すると0になるような小さな誤差が、同一の確率分布に従って独立にたくさん生成されたとします。その誤差の合計値を取ると、これは平均0の正規分布に従うはずです。
データ\( x_{ij} \)は「平均値が\( \mu_0+\delta_j \)で、分散が\( \sigma^2 \)の正規分布」に従って得られると想定するのが分散分析モデルでした。かなり単純化した議論ですが、これは「データは\( \mu_0+\delta_j \)という平均的な体重からなるが、ここに小さな誤差がたくさん積み重なることによって体重がばらつく。そしてその結果としてデータが観測される」と想定しているのが分散分析モデルだと言えるかと思います。

 

検定統計量の標本分布とp値

 

F分布

モデルが特定され、\( \mu_0,\ \delta_j,\sigma^2 \)が推定されたとします。正規分布の確率密度関数を使うことで、例えば水準Bのデータにおいて10という結果になる確率密度や30という結果になる確率密度などは即座に計算されます。今回のデータでは、30よりも大きなデータが得られる確率はとても低くなると計算されるでしょう。（水準Bのデータは平均値が11の正規分布に従って得られると想定されるから）。

ところで、分散分析という仮説検定では検定統計量であるF比が重要な役割を果たすのでした。ここで、途中経過は省きますが（大変に長い計算が要求されます）、分散分析モデルが特定されたならば、F比の従う確率分布が理論的に導出されます。この確率分布をF分布と呼びます。
F分布を使うことで、検定統計量であるF比が、例えば10を超える確率や、20を超える確率などを計算できるようになります。

正規分布の確率密度関数は、平均\( \mu \)と分散\( \sigma^2 \)の2つのパラメータで構成されていました。
F分布の確率密度関数は、「効果の大きさ(水準間の分散)の分母(今回は2)」と「誤差の大きさ(水準内の分散)の分母(今回は3)」の2つのパラメータで構成されています。
これをまとめて\( F\left(\nu_B,\nu_W\right) \)と、あるいは確率変数Xを明示して\( F\left(X\middle|\nu_B,\nu_W\right) \)と表記することにします。

F分布の確率密度関数は大変に複雑ですが下記のようになります。
F分布のパラメータは、\( \nu_B,\nu_W \)という記号を使うととても見づらいので\( m,n \)とします。この2つのパラメータは0よりも大きな整数しかとらないことに注意してください。

$$ F\left(X \middle| m, n \right) = \frac{\Gamma(\frac{m+n}{2}) m^{\frac{m}{2}} n^{\frac{n}{2}} x^{\frac{m}{2}-1} }
{\Gamma(\frac{m}{2}) \Gamma(\frac{n}{2}) (mx+n)^{\frac{m+n}{2}} } $$

ただし\( x \gt 0 \)です。ちなみに\( \Gamma(y) \)はガンマ関数であり\( \Gamma(y) = \int_0^\infty t^{y-1}e^{-t} dt \)となります。
なんだか目がチカチカしますが、この数式を覚える必要性はありません。
分散分析モデルが特定されたならば、F比の従う確率分布が理論的に導出されるということを理解しておけば、実用上は十分です。

今回はm,nが各々2，3となっていました。
このとき、F分布を使うことで、F比が16以上となる確率はおよそ0.025であると計算されます。
この数値がp値となります。

 

p値の解釈

統計的仮説検定では、p値という指標を使って、帰無仮説を受け入れるかどうかの判断を行います。
Wasserstein and Lazar(2017)『統計的有意性とP値に関するASA声明』を引用すると、p値は以下のように解釈されます。

おおざっぱにいうと、P値とは特定の統計モデルのもとで、データの統計的要約（たとえば、2グループ比較での標本平均の差）が観察された値と等しいか、それよりも極端な値をとる確率である。

データの統計的要約というのは、主に検定統計量であり、分散分析の場合はF比となります。すなわち分散分析の場合は「特定の統計モデル(分散分析モデル)のもとで、検定統計量(F比)が観察された値と等しいか、それよりも極端な値をとる確率」と解釈できます。
今回の観察データではF比が16であると計算されました。そうしたら「F比が16と等しいか、それよりも極端な値をとる確率」としてp値が計算されます。

F比の計算のときに紹介したように、F比が大きければ平均値に差があると考えられそうです。けれどもF比が大きいかそうでないかを判断する方法が問題でした。15を超えていれば大きいと言ってよいのか、20を超えなければいけないのか、決めるのは簡単ではありません。
そこで、統計的仮説検定ではp値を使います。
分散分析モデルを推定したうえで「F比が今回観察された値と等しいか、それよりも極端な値をとる確率」としてp値を計算します。この確率が小さければ、餌の効果によって体重が変化したと判断します。
「p値が十分に小さいとみなせる閾値」をどのように定めるかはやはり問題として残ります。伝統的に「p値が0.05を下回ったならば、p値を十分小さいとみなす」ことが多いです。本記事でも0.05という数値を使うことにします。

今回はp値が0.025と計算されたのでした。p値が0.5を下回っているので「餌の種類によって、体重が有意に変化した」と判断することになります。

 

5．分散分析の実行

F比やp値の計算は、コンピュータで実行するのが簡単です。
分散分析の実行方法として、R言語を使う方法とPythonを使う方法を解説します。

データやコードはGitHubで公開しています。

 

Rによる実行

Rによる実行方法を解説します。
Rの初歩的な使い方は「R言語ではじめるプログラミングとデータ分析」や「Rの簡単な使い方」を参照してください

 

分散分析の実行

まずはデータを読み込みます。read.csvという関数を使ってCSVファイルを読み込みます。読み込んだデータはdataという名前で格納しておきます。
データは「data.csv」という名前でGitHubにおいてあります。

> # データの読み込み
> data <- read.csv("data.csv")
> data
  weight food
1      2    A
2      4    A
3     10    B
4     12    B
5      6    C
6      8    C

 
分散分析モデルを推定したうえで、分散分析表を出力します。
R言語ではlm関数を使ってモデルを推定します。「weight ~ food」でモデルの構造を指定します。これをformulaと呼びます。今回は「weightはfoodによって変化する」と想定した分散分析モデルとなります。
推定されたモデルに対してanova関数を適用することで分散分析表が得られます。

> # モデル化
> anova_mod <- lm(weight ~ food, data = data)
> 
> # 分散分析表
> anova_table <- anova(anova_mod)
> anova_table
Analysis of Variance Table

Response: weight
          Df Sum Sq Mean Sq F value  Pr(>F)  
food       2     64      32      16 0.02509 *
Residuals  3      6       2                  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Analysis of Variance Tableは分散分析表と呼ばれます。
餌の効果の行と、誤差の行に分かれて、各々、自由度\( \nu_B,\nu_W \)・偏差平方和\( SS_B,SS_W \)・平均平方和\( \sigma^{2}_{B},\sigma^{2}_{W} \)・F比・p値が出力されています。
手計算した結果と一致していることを確認してください。

 
ほんの数行ですべての計算が終わってしまうので、逆に理解がし辛いかもしれません。
F分布からp値を求める計算方法も紹介しておきます（これも1行で終わりますが）。

> # F分布を使ってp値を計算する
> 1 - pf(16, 2, 3)
[1] 0.02509457

pf関数を使うことで、自由度が2と3であるF分布において、F比が16以下である確率を計算できます。1からこの確率を引くことで「F比が16以上である確率」が計算できます。
0.025なので0.05を下回っていますね。そのため今回は餌の種類によって、体重が有意に変化するとみなすことになります。

 

分散分析の結果の取り扱い

ところで分散分析表の結果はdata.frameと呼ばれる標準的なデータ格納方式が使われています。
そのためF比など個別の結果を取得するのは容易です。

> # 分散分析表はdata.frameなので、結果の取り出しは容易
> class(anova_table)
[1] "anova"      "data.frame"

 
F比を取得してみます。角カッコと列名・行名を指定することで、任意の要素を取得できます。

> # F比
> anova_table["F value"]
          F value
food           16
Residuals        
> anova_table["food", "F value"]
[1] 16

 

Pythonによる実行

続いてPythonを使います。一昔前は統計分析といえばRで機械学習はPythonと使い分けがされることもありました。
しかし、単純な統計分析であれば、PythonでもRとほとんど同様に実行ができます。せいぜいライブラリの読み込みのためのコードが増える程度です。
分散分析に関しても簡単に実行ができます。

Pythonの初歩的な使い方は「Pythonで学ぶあたらしい統計学の教科書」や「Pythonの簡単な使い方」を参照してください。

 

前準備

Pythonは汎用的なプログラミング言語です。統計分析を行う場合は特別なライブラリを別途読み込んで使うことになります。
ライブラリの読み込みコードは以下のようになります。分散分析を実行するだけならばもっと少なくても良いのですが、後ほどPB検定を実行するために少し増やしました。

# ライブラリの読み込み
import numpy as np
import pandas as pd
import scipy as sp
import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf

参考までにライブラリのバージョンを載せておきます。

# ライブラリのバージョン
print('np ',np.__version__)
print('pd ',pd.__version__)
print('sp ',sp.__version__)
print('sm ',statsmodels.__version__)

結果は以下の通りです。

np  1.18.5
pd  1.0.5
sp  1.5.0
sm  0.11.1

 

分散分析の実行

データを読み込みます。pd.read_csv関数を使います。
データは「data.csv」という名前でGitHubにおいてあります。
実装コードと出力を分けるのが手間なので点線で区切ってあります。

# データの読み込み
data = pd.read_csv('data.csv')
print(data)
-------------------------------------------------
   weight food
0       2    A
1       4    A
2      10    B
3      12    B
4       6    C
5       8    C

 
分散分析モデルを推定したうえで、分散分析表を出力します。
モデルの推定にはsmf.ols関数を使います。最後の「.fit()」を忘れないように注意しましょう。
分散分析表の出力にはsm.stats.anova_lm関数を使います。

# モデル化
anova_mod = smf.ols('weight ~ food', data=data).fit()

# 分散分析表
anova_table = sm.stats.anova_lm(anova_mod, type=2)
print(anova_table)
-------------------------------------------------
           df  sum_sq  mean_sq     F    PR(>F)
food      2.0    64.0     32.0  16.0  0.025095
Residual  3.0     6.0      2.0   NaN       NaN

 
F分布を使ってp値を計算してみます。

# F分布を使ってp値を計算する
1 - sp.stats.f.cdf(16, 2, 3)
-------------------------------------------------
0.02509457330439091

ほぼRと同じように実装できます。もちろん結果も同じです。

 

分散分析の結果の取り扱い

Pythonにおいても、分散分析表はDataFrameと呼ばれる標準的なデータ形式が使われているので、結果の取り出しは容易です。

# 分散分析表はDataFrameなので結果の取り出しは容易
type(anova_table)
-------------------------------------------------
pandas.core.frame.DataFrame

 
F比を取り出します。

# F比
anova_table['F']
-------------------------------------------------
food        16.0
Residual     NaN
Name: F, dtype: float64

行名と列名を指定するのはRと同じですがloc関数をかませることを忘れないように注意しましょう。

anova_table.loc['food', 'F']
-------------------------------------------------
16.0

わざわざ分散分析表を作らなくても、推定されたモデルから直接F比を取得することもできます。

# モデルから直接取り出してもよい
anova_mod.fvalue
-------------------------------------------------
16.0

 

6．パラメトリックブートストラップ検定（PB検定）

統計的仮説検定についてより深く理解するために、分散分析という検定方法とは少し異なるのですが、パラメトリックブートストラップ検定(PB検定)という検定方法を解説します。
先ほどの計算ではF分布という便利な確率密度関数を使うことでp値を簡単に計算できました。プログラムも短くて済むのですが、あまりにも簡単に計算が終わってしまって理解が追い付かないことがあるかもしれません。ここではあえてシミュレーションという面倒くさい作業を通して、p値についての理解を深めていただこうと思います。

PB検定では、シミュレーションを行ってF比の標本分布を得て、それに基づいてp値を計算します。F分布という複雑な確率分布を使わなくても済むので、直観的にも理解がしやすいと思います。また統計モデルに従ってデータを生成するという作業を通して、「データが得られる過程」に関する理解も深めてほしいです。

なお、ここではR言語による実装を載せていますが、Pythonでも実行できます。Pythonの結果はこちらから確認できます。

 

餌の種類によって体重が変わらないと想定したモデル

p値の解釈として、分散分析の場合は「特定の統計モデル(分散分析モデル)のもとで、検定統計量(F比)が観察された値と等しいか、それよりも極端な値をとる確率」であると紹介しました。
この言葉をすんなり理解できれば苦労はないのですが、なかなか難しいと感じる人も多いかと思います。
そこであえてPB検定というシミュレーションを使った検定手法を使い、p値の意味を学んでいただくのがこの章のテーマです。

まずは「餌の種類によって体重が変わらないと想定したモデル」を推定してみましょう。これをNullモデルと呼びます。
以下のようにして推定できます。

# foodによって体重が変わらないことを仮定したモデル
anova_null <- lm(weight ~ 1, data = data)

通常の分散分析モデルのコードと比較します。1行目が通常の分散分析モデルで、2行目がNullモデルです。

anova_mod  <- lm(weight ~ food, data = data) # 通常
anova_null <- lm(weight ~ 1   , data = data) # Nullモデル

モデルのformulaにおいて、Nullモデルではfoodが使われていません。
これでfoodによって体重が変わらないことを仮定したモデルが推定できます。

Nullモデルのモデル式は以下のようになります。餌の種類の違いがもたらす効果を表す項がなくなりました。

$$ x_{ij}=\mu_0+\varepsilon_{ij},\ 　\varepsilon_{ij} \sim N(0,\sigma^2) $$

推定された\( \mu_0 \)は以下のようにして取得できます。

> # 推定されたパラメータ
> coef(anova_null)
(Intercept) 
          7 

これは体重の総平均と一致します。

> # 総平均
> mean(data$weight)
[1] 7

推定された標準偏差\( \sigma \)は以下のようになります。

> # モデルの標準偏差
> summary(anova_null)$sigma
[1] 3.741657

これは、体重の標準偏差と一致します。

> # 体重の標準偏差
> sd(data$weight)
[1] 3.741657

Nullモデルのモデル式を変形すると以下のようになります。

$$ x_{ij} \sim N\left(\mu_0,\sigma^2\right) $$

データは平均が\( \mu_0 \)で分散が\( \sigma^2 \)である正規分布から生成されたとみなしているのがNullモデルです。
このため、モデルによって推定されたパラメータは、元の体重データの総平均および標準偏差と一致します。

 

シミュレーションによる乱数の生成

続いて、このNullモデルを使って、疑似的に体重データを生成してみます。
simulate関数を使うことで、疑似データをシミュレーションで生成できます。シミュレーションで生成されたデータを乱数と呼ぶことにします。
引数のn_simは、生成するデータセットの数です。今回はサンプルサイズが6であるデータに対してNullモデルを推定したので、n_sim=1とすると6つの乱数が生成されます。

> # 帰無仮説が正しいと仮定して、データを生成する
> simulate(anova_null, n_sim = 1)
       sim_1
1  1.7894269
2  0.7065733
3  4.4979244
4 12.6246366
5  5.3855722
6  4.4790944
> simulate(anova_null, n_sim = 1)
      sim_1
1 11.807013
2  2.514968
3  2.832782
4  8.404657
5  1.549460
6  7.873122

これからNullモデルによって生成された乱数を使って、p値の意味を確認していきます。
その前に、乱数生成において少し工夫をします。
同じ確率分布に従う乱数であっても、乱数を生成するたびにその結果は変化します。例えばサイコロを何度も投げたとき、2の目が出たり5の目が出たり、いろいろと変化するような感じです。
毎回結果が変わると不便なので、生成される乱数を固定します。それがseedという引数です。指定する数値は何でも良いですが、とりあえず1にしておきます。同じ数値のseedを指定すると、同じ乱数が生成されます。

> # 乱数を固定する
> simulate(anova_null, n_sim = 1, seed = 1)
      sim_1
1  4.656024
2  7.687130
3  3.873364
4 12.968994
5  8.232905
6  3.930088
> simulate(anova_null, n_sim = 1, seed = 1)
      sim_1
1  4.656024
2  7.687130
3  3.873364
4 12.968994
5  8.232905
6  3.930088

 
simulate関数を使わない方法も紹介します。
正規分布に従う乱数を生成する関数であるrnorm関数を使います。

anova_nullでは、総平均7と標準偏差3.741657がパラメータとして推定されていました。
このNullモデルに従う乱数を生成することは、「平均7、標準偏差3.741657の正規分布」に従う確率変数を生成しているのと同じです。
それを以下のコードで確認します。
rnorm関数は、引数nで生成する乱数の数を、meanで正規分布の平均値を、sdで標準偏差を指定します。
set.seed関数で乱数の結果を固定しています。

> # 平均7、標準偏差3.741657の正規分布に従う確率変数を生成しているのと同じ
> set.seed(1)
> mu <- mean(data$weight)
> sd <- sd(data$weight)
> size <- nrow(data)
> rnorm(n = size, mean = mu, sd = sd)
[1]  4.656024  7.687130  3.873364 12.968994  8.232905  3.930088

simulate関数の結果と一致していることを確認してください。

 
ここで、Nullモデルから生成された乱数は、餌の影響を一切受けていません(そのように想定してモデルを作ったからです)。
そのため、Nullモデルから生成された乱数に対して分散分析を適用すると、小さなF比が得られることが予想されます。
実際にやってみます。

> # このシミュレーションデータに対してF比を計算すると、
> # 小さなF比が得られやすい
> set.seed(1)
> res_rnorm <- rnorm(n = size, mean = mu, sd = sd)
> mod_sim <- lm(res_rnorm ~ data$food)
> anova(mod_sim)
Analysis of Variance Table

Response: res_rnorm
          Df Sum Sq Mean Sq F value Pr(>F)
data$food  2  7.029  3.5143  0.1909 0.8355
Residuals  3 55.216 18.4054      

予想通り、F比は0.1909と小さな値になりました。

とはいえ、乱数は（set.seed関数を使って乱数を固定しない限り）変化します。
無数に乱数を作れば、時々は大きめのF比が観測されるはずです。

さて、ここで今回のNullモデルから生成された乱数に対してF比を計算し、その結果が16を超える確率はいくらになるでしょうか。
この確率こそがp値です。

 

PB検定の実行

PB検定を実行してp値を計算してみましょう。
まずはシミュレーションの設定をします。

シミュレーションにより乱数を50000セット(1セットにつきサンプルサイズは6)作ることにします。
多ければ多い方が安定した結果が得られますが、回数を多くすると時間がかかります。
今回のPB検定の実装コードがあまり効率が良くないこともありますが、当方のPCでは1分以上かかりました。
あとはコメントの通りですが、f_ratio_vecで5万個のF比を保管する容れ物を用意します。
水準に関してはexp_vecという名前で用意しておきます。

n_sim <- 50000                # シミュレーションする回数
f_ratio_vec <- numeric(n_sim) # F比を保管する容れ物
exp_vec <- data$food          # 水準

 
PB検定の実行コードです。やや時間がかかるので注意してください。

set.seed(1)
for (i in 1:n_sim) {
  # シミュレーションにより体重データを生成
  # このデータは、体重が餌によって変化しないことを想定している
  simlated_weight <- rnorm(n = size, mean = mu, sd = sd)
  # モデル化と分散分析表の出力
  mod_sim <- lm(simlated_weight ~ exp_vec)
  table_sim <-  anova(mod_sim)
  # F比を保管
  f_ratio_vec[i] <- table_sim["exp_vec", "F value"]
}

1行目で乱数を固定しています。乱数の固定は最初の1回だけ行います。forループの中に入れないように注意してください。
2行目でforループと呼ばれる繰り返し構文を利用しています。for (i in 1:n_sim)で、iを1からn_simまで変化させながら、中カッコ内の処理を繰り返し実行します。
「simlated_weight <- rnorm(n = size, mean = mu, sd = sd)」で乱数を生成しています。これはNullモデルからの乱数生成とみなせます。
そしてmod_simとして分散分析モデルを推定し、table_simとして分散分析表を得ます。
最終的に「table_sim[“exp_vec”, “F value”]」として得られるF比をf_ratio_vecに格納します。

5万個あるf_ratio_vecの中で、F比が16以上であった割合を計算します。

> # F比が16以上である割合
> sum(f_ratio_vec >= 16) / n_sim
[1] 0.02478

0.02478となりました。これがp値です。F分布を用いた方法(およそ0.025)とほぼ同じ結果になりました。

PB検定は、やや手間がかかるものの、p値の意味が比較的理解しやすい手法ではないかと思います。

 
ところで、分散分析という検定手法を適用する際にしばしばいわれるのが「等分散正規分布の仮定」です。
これは母集団分布が正規分布に従っていなくてはならないという仮定と、水準ごとで（平均値が変わっていても）分散は変わらないことを仮定していることをまとめて「等分散正規分布の仮定」と呼ぶようです。
このような仮定の名前を丸暗記することはぜひ避けてほしいなと思います。
分散分析という検定手法の背後にある「統計モデル」について自覚的になれば「等分散正規分布の仮定」といった標語を覚える必要はありませんね。
PB検定を実行すると、モデルの仮定についてより理解が深まるかと思います。

 

参考文献

平均・分散から始める一般化線形モデル入門   このブログの管理人が書いた本です。Rを使った分析方法を解説しています。 一般化線形モデルというのは、分散分析をはじめとしたさまざまな分析手法を体系立てたものです。 サポートページはこちら。

Pythonで学ぶあたらしい統計学の教科書   このブログの管理人が書いた本です。Pythonを使った分析方法を解説しています。 一般化線形モデルの解説も含まれます。 サポートページはこちら。

生物統計学入門   さまざまな確率分布の解説や、数理的な側面の丁寧な解説が載っている、統計学の本です。 入門と銘打っていますが、やや難易度は高めです。

一般化線形モデル入門 原著第2版   一般化線形モデルの数理に興味がある人のための教科書です。

統計学辞典   名前の通り辞典です。わからない用語があった時はこれを参照すると便利です。

 
Wasserstein RL, Lazar NA．(佐藤俊哉 訳)．(2017)．統計的有意性とP値に関するASA声明．
→統計的仮説検定の使い方を学ぶ資料として個人的にお勧めしています。

 

 
更新履歴
2020年12月19日：新規作成