からだにいいもの Rで解析：計数データの過分散に対応する「hermite」パッケージ
発生件数のような計数データを扱うとき、どの分布を当てはめるかが解析の出発点になります。しかし、事象が1件ずつだけでなく複数まとまっても発生する場合、ポアソン分布では分散を説明しきれず、過分散への対処に手間がかかります。
「hermite」パッケージは、一般化エルミート分布にもとづく確率計算と回帰分析を簡単におこなえるパッケージです。確率質量関数、累積分布関数、分位点関数、乱数生成の4つのコマンドが収録されており、分布の形状の確認からシミュレーションまでを一貫して扱えます。
また、説明変数を含む回帰モデルの最尤推定や、ポアソン分布を帰無仮説とする尤度比検定も可能です。本パッケージの利用で、まとまって発生する計数データの特性を反映した解析ができるのではないかと考えます。
パッケージバージョンは1.2.1。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。
＜おすすめのRに関する書籍です＞
![Rユーザのためのtidymodels[実践]入門〜モダンな統計・機械学習モデリングの世界](https://m.media-amazon.com/images/I/51t4FWJKkuL._SL500_.jpg)
www.amazon.co.jp
パッケージのインストール
下記コマンドを実行してください。依存パッケージとしてmaxLikが同時にインストールされます。
# パッケージのインストール
install.packages("hermite")
# パッケージの読み込み
library("hermite")コマンド例
詳細はコメント、パッケージのヘルプを確認してください。
本記事では、架空の動物病院における30日分の1日あたりのペット数を作成し、記事全体を通して使用します。病院には1頭で来院するケースと、多頭飼いのため2頭同時に来院するケースが混在していることを考えました。1件ずつの事象と2件まとまった事象が重なる構造は、次数m = 2の一般化エルミート分布そのものです。
# 架空の動物病院における30日分の1日あたり来院ペット数を作成
raiin <- c(1, 3, 0, 7, 2, 9, 4, 0, 5, 11,
2, 1, 6, 13, 3, 0, 8, 2, 5, 15,
4, 1, 7, 0, 3, 10, 6, 2, 9, 4)
# 各日が休診日の翌日かどうかを示すフラグを作成
yokujitsu <- c(0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
0, 0, 1, 1, 0, 0, 1, 0, 0, 1,
0, 0, 1, 0, 0, 1, 1, 0, 1, 0)
# 平均と分散を比較して過分散の有無を確認
c(平均 = mean(raiin), 分散 = var(raiin), 分散指標 = var(raiin) / mean(raiin))
平均 分散 分散指標
4.766667 16.322989 3.424403 分散が平均の約3.4倍となっており、ポアソン分布が前提とする「分散＝平均」から大きく外れています。以降、このデータに一般化エルミート分布を当てはめます。
一般化エルミート分布の確率質量関数を計算：dhermiteコマンド
| オプション | 意味 | 初期値 |
|---|---|---|
| x | 確率を求める非負整数のベクトル | なし |
| a | エルミート分布の第一パラメーター | なし |
| b | エルミート分布の第二パラメーター | なし |
| m | 一般化エルミート分布の次数、2が標準的なエルミート分布に対応 | 2 |
指定した値がちょうど発生する確率を返します。a、b、mのいずれかが負の値、またはmが2未満や非整数の場合は警告とともにNaNが返ります。
# 単独来院が平均2件、ペア来院が平均1.4件のときの確率質量関数を計算
dhermite(0:8, a = 2, b = 1.4, m = 2)
[1] 0.03337327 0.06674654 0.11346912 0.13794285 0.14839981 0.13660792 0.11478922 0.08744009 0.06203625
# 3頭ちょうど来院する確率のみを取得
dhermite(3, a = 2, b = 1.4, m = 2)
[1] 0.1379428mを変更すると、まとまって発生する件数が変わります。m = 3では3件単位の来院を仮定した分布になります。
# 3頭まとまって来院すると仮定した場合の確率質量関数を計算
round(dhermite(0:8, a = 2, b = 1.4, m = 3), 4)
[1] 0.0334 0.0667 0.0667 0.0912 0.1157 0.1023 0.0980 0.0974 0.0781一般化エルミート分布の累積分布関数を計算：phermiteコマンド
| オプション | 意味 | 初期値 |
|---|---|---|
| q | 累積確率を求める非負整数のベクトル | なし |
| a | エルミート分布の第一パラメーター | なし |
| b | エルミート分布の第二パラメーター | なし |
| m | 一般化エルミート分布の次数、2が標準的なエルミート分布に対応 | 2 |
| lower.tail | TRUEでq以下となる確率、FALSEでqを超える確率を返す | TRUE |
# 来院数が0頭から8頭以下となる累積確率を計算
round(phermite(0:8, a = 2, b = 1.4, m = 2), 4)
[1] 0.0334 0.1001 0.2136 0.3515 0.4999 0.6365 0.7513 0.8388 0.9008
# 8頭を超える確率を計算
phermite(8, a = 2, b = 1.4, m = 2, lower.tail = FALSE)
[1] 0.09919495
＜おすすめのRに関する書籍です＞
www.amazon.co.jp
指定した確率に対応する分位点を計算：qhermiteコマンド
| オプション | 意味 | 初期値 |
|---|---|---|
| p | 分位点を求める確率のベクトル | なし |
| a | エルミート分布の第一パラメーター | なし |
| b | エルミート分布の第二パラメーター | なし |
| m | 一般化エルミート分布の次数、2が標準的なエルミート分布に対応 | 2 |
| lower.tail | TRUEで累積確率がp以上となる最小の整数、FALSEで1-pに対応する値を返す | TRUE |
累積確率が指定した値以上となる最小の整数を返します。準備すべき診察枠の見積もりなどに利用できます。
# 中央値から99パーセント点までの分位点を計算
qhermite(c(0.5, 0.8, 0.9, 0.95, 0.99), a = 2, b = 1.4, m = 2)
[1] 5 7 8 10 12
# 上側5パーセント点を計算
qhermite(0.05, a = 2, b = 1.4, m = 2, lower.tail = FALSE)
[1] 10一般化エルミート分布に従う乱数を生成：rhermiteコマンド
| オプション | 意味 | 初期値 |
|---|---|---|
| n | 生成する観測数 | なし |
| a | エルミート分布の第一パラメーター | なし |
| b | エルミート分布の第二パラメーター | なし |
| m | 一般化エルミート分布の次数、2が標準的なエルミート分布に対応 | 2 |
内部では平均aのポアソン乱数と、平均bのポアソン乱数をm倍したものの和が返されます。
# 乱数の再現性を確保するためにシード値を固定
set.seed(1234)
# 300日分の来院数を生成
rnd <- rhermite(300, a = 2, b = 1.4, m = 2)
# 先頭20件を確認
head(rnd, 20)
[1] 2 10 4 4 12 4 0 11 4 6 9 6 3 6 1 3 1 3 3 7
# 理論値の平均4.8および分散7.6と比較
c(平均 = mean(rnd), 分散 = var(rnd))
平均 分散
5.103333 8.447480エルミート回帰モデルの最尤推定：glm.hermiteコマンド
| オプション | 意味 | 初期値 |
|---|---|---|
| formula | 目的変数と説明変数を記述したモデル式 | なし |
| data | モデルに含まれる変数を格納したデータフレーム | なし |
| link | リンク関数の指定、”log”または”identity” | “log” |
| start | パラメーターの初期値ベクトル | NULL |
| m | 次数mの値、NULLの場合は関数内で推定 | NULL |
戻り値はcoefs（係数）、loglik（対数尤度）、vcov（分散共分散行列）、hess（ヘッセ行列）、fitted.values（当てはめ値）、w（尤度比検定統計量）、pval（そのp値）を含むリストです。wとpvalは、ポアソン分布を帰無仮説とする尤度比検定の結果です。
glm.hermiteコマンド実行時に
警告メッセージ:
model.matrix.default(mt, mf, contrasts) で:
non-list contrasts argument ignored
が表示されますが問題ありません。
# 説明変数を含まないモデルで来院数の分布を推定
mod0 <- glm.hermite(raiin ~ 1, link = "log", m = 2)
# 切片、分散指標、次数の推定値を確認
mod0$coefs
(Intercept) dispersion.index order
1.561707 1.811578 2.000000
# ポアソン分布に対する尤度比検定の結果を確認
c(検定統計量 = mod0$w, p値 = mod0$pval)
検定統計量 p値
2.805052e+01 5.909459e-08 分散指標が1.81と1を上回り、p値も十分に小さいため、ポアソン分布よりも一般化エルミート分布のほうが当てはまりがよいと判断できます。なお、本コマンドは実行のたびに「non-list contrasts argument ignored」という警告を表示しますが、パッケージ内部の処理によるもので推定結果には影響しません。
説明変数を加える場合は、dataオプションにデータフレームを渡します。結果はsummaryコマンドで整形して確認できます。
# 来院数と休診日翌日フラグをデータフレームにまとめる
byoin <- data.frame(raiin = raiin, yokujitsu = yokujitsu)
# 休診日翌日フラグを説明変数としてエルミート回帰を実行
mod1 <- glm.hermite(raiin ~ yokujitsu, data = byoin, link = "log", m = 2)
警告メッセージ:
model.matrix.default(mt, mf, contrasts) で:
non-list contrasts argument ignored
# 推定結果を整形して表示
summary(mod1)
Call:
glm.hermite(formula = raiin ~ yokujitsu, data = byoin, link = "log",
m = 2)
Deviance Residuals:
Min 1Q Median 3Q Max
-8.201381 -8.201381 -2.199341 -1.199341 -1.199341
Coefficients:
Estimate Std. Error z value p-value
(Intercept) 0.7881579 0.1677142 4.6994103 2.609138e-06
yokujitsu 1.4311958 0.1994811 7.1745942 7.252198e-13
dispersion.index 1.1717825 0.4338111 0.1815149 3.350373e-01
order 2.0000000 NA NA NA
(Likelihood ratio test against Poisson is reported by *z value* for *dispersion.index*)
AIC: 130.7737 休診日翌日の係数が1.43と正の値になっており、対数リンクのもとで来院数が約4.2倍に増えると読み取れます。
＜おすすめのRに関する書籍です＞
www.amazon.co.jp
実行例
パッケージの機能を組み合わせた、実践的な使い方の例です。
推定結果からパラメーターを復元して実測度数と比較する
glm.hermiteが返すのは平均の対数と分散指標であり、dhermiteが必要とするaとbではありません。平均をμ、分散指標をdとすると、m = 2のときb = μ(d – 1) / 2、a = μ – 2bで復元できます。復元したパラメーターを使い、実測の度数分布、一般化エルミート分布、ポアソン分布の3者を並べて比較します。
# 推定された切片を指数変換して平均を取得
heikin <- exp(unname(mod0$coefs[1]))
# 推定された分散指標を取得
shisu <- unname(mod0$coefs[2])
# 分散指標と平均から第二パラメーターbを復元
b_hat <- heikin * (shisu - 1) / 2
# 平均とbから第一パラメーターaを復元
a_hat <- heikin - 2 * b_hat
# 復元したパラメーターを小数第4位に丸めて確認
round(c(a = a_hat, b = b_hat), 4)
a b
0.8982 1.9344
# 実測度数と2つの理論分布を並べた比較表を作成
hikaku <- data.frame(
来院数 = 0:10,
実測割合 = as.numeric(table(factor(raiin, levels = 0:10))) / length(raiin),
エルミート = dhermite(0:10, a = a_hat, b = b_hat, m = 2),
ポアソン = dpois(0:10, lambda = heikin)
)
# 小数第4位に丸めて表示
round(hikaku, 4)
来院数 実測割合 エルミート ポアソン
1 0 0.1333 0.0589 0.0085
2 1 0.1000 0.0529 0.0405
3 2 0.1333 0.1376 0.0966
4 3 0.1000 0.1094 0.1536
5 4 0.1000 0.1576 0.1830
6 5 0.0667 0.1130 0.1745
7 6 0.0667 0.1186 0.1386
8 7 0.0667 0.0776 0.0944
9 8 0.0333 0.0661 0.0563
10 9 0.0667 0.0400 0.0298
11 10 0.0333 0.0291 0.01420頭の日と9頭以上の日のいずれについても、ポアソン分布は確率を過小評価しています。一般化エルミート分布は両端の確率を高く見積もっており、実測の裾の重さをよく捉えています。
分位点から必要な診察枠を見積もる
復元したパラメーターをqhermiteに渡すと、指定した割合の日をカバーできる来院数がわかります。診察枠や在庫の準備数を決める場面で利用できます。
# カバーしたい日の割合を設定
wariai <- c(0.5, 0.8, 0.9, 0.95, 0.99)
# 各割合に対応する来院数を計算
waku <- qhermite(wariai, a = a_hat, b = b_hat, m = 2)
# 割合と必要な診察枠を対応させた表を作成
data.frame(カバー率 = wariai, 必要枠 = waku)
カバー率 必要枠
1 0.50 4
2 0.80 7
3 0.90 9
4 0.95 10
5 0.99 13
# 12頭を超える日が発生する確率を計算
round(phermite(12, a = a_hat, b = b_hat, m = 2, lower.tail = FALSE), 4)
[1] 0.012395パーセントの日をカバーするには10頭分、99パーセントでは13頭分の枠が必要という見積もりが得られます。
乱数生成によるシミュレーションで見積もりを検証する
rhermiteで大量の仮想データを生成し、分位点から得た見積もりが妥当かを確認します。
# 乱数の再現性を確保するためにシード値を固定
set.seed(20260804)
# 復元したパラメーターで10000日分の来院数を生成
sim <- rhermite(10000, a = a_hat, b = b_hat, m = 2)
# 生成データの平均と分散を確認
c(平均 = mean(sim), 分散 = var(sim))
平均 分散
4.760300 8.886933
# 12頭を超える日の実際の割合を確認
mean(sim > 12)
[1] 0.0132シミュレーションによる平均4.76は、実測の平均4.77とほぼ一致しています。12頭を超える日の割合は約1.3パーセントで、phermiteから求めた理論値の約1.2パーセントとよく対応しています。
この記事が誰かの役に立ちますように。
