R言語でポアソン検定を実行する関数とその実行例を紹介します。
あるイベント数に関するデータがあるときに、イベントの単位時間当たりの発生頻度を検定する際に用いるポアソン検定について見ていきます。
この記事で扱うRスクリプトは以下からダウンロードできます。
ポアソン検定を詳しく見たい方は次の記事を参照してください。
-
【統計学】ポアソン検定 ポアソン分布に関する正確検定
ポアソン検定について解説する。 ポアソン過程に関するデータの平均時間当たりのイベントの発生数の検定であるポアソン検定の手順や検定の例を紹介する。 検定の例では、ポアソン検定のp値の具体的な計算例を紹介 ...
続きを見る
ポアソン検定
ポアソン検定の手順と検定を実行する関数を以下にまとめました。
検定の概要
ポアソン検定は以下のポアソン過程の事象に関する正確検定となります。
観測開始から終了までのイベント発生数に関して、単位人当たりにどれだけイベントが発生するかを検定するときに次のポアソン検定を用います。
ポアソン検定
\(x\)をパラメータ\(\lambda\)のポアソン分布\(Pois(\lambda)\)からの観測値とする。次の仮説検定を考える。
有意水準\(\alpha\)の仮説検定を行うために、次のポアソン分布の両側確率から計算されるp値を用いる。
ここに
上の集合は、\(x\)に対応する確率変数\(X\)が観測値\(x\)よりも極端な値をとる確率を意味する。
関数poisson.test
次の関数poisson.testでポアソン検定を実行することができます。
poisson.testの引数は以下の通り。
| x | numeric型。イベントの発生数。長さが1または2のベクトルで、長さが2のときは二標本検定となる。 |
| T | numeric型。イベントの発生数についての単位時間。長さが1または2のベクトル。 |
| r | numeric型。帰無仮説の下での単位時間当たりのイベントの発生頻度、または頻度の比。 |
| alternative | character型。"two.sided"で両側検定、"greater"で右片側検定、"less"で左片側検定を実行できる。 |
| conf.level | numeric型。信頼区間を計算する際の信頼水準。 |
実行例
ポアソン検定の実行例をいくつかみていきます。
簡単な例からデータセットを用いた実行例も紹介します。
簡単な例
まずは簡単な例について見ていきます。ポアソン検定がどのような検定なのかを簡単に解説します。
次の10個のイベントが単位時間当たり（例えば1時間）で得られたとします。
1 2 | #ポアソン検定 numberOfEvents <- 10 |
「単位時間当たりのイベントの発生数は8であるか」という仮説検定を実行するには、次のように関数poisson.testの引数x, T, rにイベントの発生数、単位時間、帰無仮説の下でのイベントの発生頻度を渡します。
1 2 | lambda0 <- 8 testResult <- poisson.test(numberOfEvents, 1, lambda0) |
testResultを参照すると、ポアソン検定の結果がコンソールに表示されます。
1 2 3 4 5 6 7 8 9 10 11 12 | > testResult Exact Poisson test data: numberOfEvents time base: 1 number of events = 10, time base = 1, p-value = 0.4746 alternative hypothesis: true event rate is not equal to 8 95 percent confidence interval: 4.795389 18.390356 sample estimates: event rate 10 |
「p-value = 0.4746」であることから、帰無仮説は採択されて単位時間当たりのイベント発生数は8であるという結論が得られました。
また、ポアソン検定のp値やイベントの発生数の信頼区間は、次のようにtestResultから取得することができます。
1 2 3 4 5 | events <- testResult$statistic #イベント数 timeBase <- testResult$parameter #イベント数に対する時間の基準 pValue <- testResult$p.value #p値 ci <- testResult$conf.int #平均時間当たりのイベント数の信頼区間 nullValue <- testResult$null.value #帰無仮説の下での平均時間当たりのイベント数 |
この検定の各種結果は次の通りです。
1 2 3 4 5 | > resultTable <- data.frame(No.events = events, time.base = timeBase, p.value = pValue, + c.i.lower = ci[1], c.i.upper = ci[2], null.value = nullValue, row.names = NULL) > resultTable No.events time.base p.value c.i.lower c.i.upper null.value 1 10 1 0.4746118 4.795389 18.39036 8 |
データセットを用いた例
次にデータセットを用いた例を紹介します。
次のデータセットを用いたポアソン検定の実行例をみていきます。
1 2 3 | #データセットを用いた例 dataset <- data.frame(Jan = 1, Feb = 6, Mar = 8, Apr = 5, May = 9, Jun = 3, Jul = 5, Aug = 8, Sep = 2, Oct = 4, Nov = 4, Dec = 5, row.names = "no.accidents") |
次のように1月から12月までの交通事故の発生件数に関するデータとなります。
1 2 3 | > dataset Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec no.accidents 1 6 8 5 9 3 5 8 2 4 4 5 |
「1か月あたりの事故の発生回数は10回であるか>という仮説検定をポアソン検定により実行したいと思います。
1か月あたりの事故件数に関する検定を行うには、poisson.testの引数を次のように設定します。xには1月から12月までの事故件数を渡し、xが12カ月の間に起こったイベントであるため、Tには12を渡します。rには帰無仮説の下での1か月あたりの事故発生数を与えてあげます。
1 2 | lambda0 <- 10 testResult <- poisson.test(sum(dataset), ncol(dataset), lambda0) |
testResultを参照すると、この検定のp値や1カ月当たりの事故発生数の信頼区間が出力されます。
1 2 3 4 5 6 7 8 9 10 11 12 | > testResult Exact Poisson test data: sum(dataset) time base: ncol(dataset) number of events = 60, time base = 12, p-value = 1.922e-09 alternative hypothesis: true event rate is not equal to 10 95 percent confidence interval: 3.815527 6.435991 sample estimates: event rate 5 |
「p-value = 1.992e-09 < 0.05」より、 有意水準0.05の下で帰無仮説は棄却され、1カ月当たりの事故発生数は10であるとはいえないという結論が得られました。
上の例で紹介したように、検定結果は次のようにtestResultの$を付けることで取得することができます。
1 2 3 4 5 | events <- testResult$statistic #イベント数 timeBase <- testResult$parameter #イベント数に対する時間の基準 pValue <- testResult$p.value #p値 ci <- testResult$conf.int #平均時間当たりのイベント数の信頼区間 nullValue <- testResult$null.value #帰無仮説の下での平均時間当たりのイベント数 |
最後に、検定結果をcsvファイルなどの他のファイルに保存する例を紹介します。
csvファイルへ保存するには、次のようにこれらの検定結果をデータフレームに格納した後に、関数write.csvで検定結果のデータフレームを出力すればおｋです。下の画像のような検定結果のファイルが出力されます。
1 2 3 | resultTable <- data.frame(No.events = events, time.base = timeBase, p.value = pValue, c.i.lower = ci[1], c.i.upper = ci[2], null.value = nullValue, row.names = NULL) write.csv(resultTable, "ポアソン検定.csv", row.names = FALSE) |
ポアソン検定結果
まとめ
R言語でポアソン検定を実行する関数と実行例について見てきました。
関数poisson.testによってポアソン検定を行うことができます。
イベントの発生確率が小さく観測値数が非常に多い場合には、関数binom.testではなくpoisson.testを用いましょう。
