ポアソン回帰とは ― 件数・発生率データをRで解析する一般化線形モデル入門 ―
この記事でわかること
- ポアソン回帰とは何か、どのような「件数・発生率データ」に使うのか
- 対数リンク関数の意味と、係数を \(\exp(\beta)\) で「率比」として解釈する考え方
- 線形回帰・ロジスティック回帰・ポアソン回帰の使い分け
- 観察期間や曝露量が異なるときに用いる「オフセット」の役割
- Rでポアソン回帰を実装する際の流れと、過分散への注意点
はじめに
臨床試験や製薬領域のデータ解析では、「一定期間に起きた有害事象の件数」「1年間の入院回数」「1か月あたりのてんかん発作回数」といった、数え上げられる件数（カウント）データが頻繁に登場します。このような0以上の整数値をとるデータに対して、普通の線形回帰（重回帰）をそのまま当てはめると、予測値が負になったり、非整数の不自然な値を返したりと、データの性質と噛み合わない不都合が生じます。そこで標準的な選択肢となるのが、本記事で解説するポアソン回帰（Poisson regression）です。ポアソン回帰は、件数データの背後にあるポアソン分布を仮定し、一般化線形モデル（GLM）の枠組みで発生率をモデル化する手法です。この記事では、件数・発生率データを扱う生物統計担当者やR実務家、統計を学ぶ学生の方に向けて、まず概念と数理的背景を丁寧に整理し、続くパートでRによる実装と結果の解釈へとつなげていきます。
ポアソン回帰とは
ポアソン回帰は、その名の通りポアソン分布を確率モデルの土台に置いた回帰手法です。ポアソン分布は、「一定の時間・空間・単位のなかで、比較的まれに起こる事象が何回発生するか」という件数の分布を表します。例えば「1年間で1人の患者に何回発作が起こるか」「ある病棟で1か月に何件の転倒が起きるか」といった状況が典型例です。
ポアソン回帰では、件数の期待値（平均発生数）\(E[Y]\) を、説明変数の線形結合と結びつけるために対数リンク関数（log link）を使います。すなわち、
\[ \log E[Y] = \beta_0 + \beta_1 x \]
という形でモデル化します。ここで期待値そのものではなく「その対数」を線形式に対応させている点が重要です。対数をとることで、右辺がどんな値をとっても \(E[Y]=\exp(\beta_0+\beta_1 x)\) は必ず正になり、件数データが負にならないという性質と自然に整合します。
どの回帰を使うべきかは、アウトカム（結果変数）の型で決まります。線形回帰・ロジスティック回帰・ポアソン回帰の違いを整理すると、次のようになります。
| 回帰の種類 | アウトカムの型 | 仮定する分布 / リンク関数 | 代表例 |
|---|---|---|---|
| 線形回帰 | 連続量 | 正規分布 / 恒等リンク | 血圧・体重・検査値 |
| ロジスティック回帰 | 二値（0 / 1） | 二項分布 / ロジットリンク | 有害事象の有無・治癒の有無 |
| ポアソン回帰 | 件数（0以上の整数） | ポアソン分布 / 対数リンク | 発作回数・入院件数・事象数 |
このように、「数えられる件数」「単位期間あたりの発生率」を扱いたい場面こそ、ポアソン回帰の出番だと押さえておきましょう。
ポアソン回帰の数理的背景
ポアソン分布の確率関数（確率質量関数）は、平均パラメータ \(\lambda\) を用いて次のように書けます。
\[ P(Y=y)=\dfrac{\lambda^{y}e^{-\lambda}}{y!} \]
ここで \(y=0,1,2,\dots\) は観測される件数、\(\lambda>0\) は「1単位あたりに平均して何回起こるか」を表す発生強度です。この式は、平均が \(\lambda\) 回である事象が、実際に \(y\) 回観測される確率を与えています。
ポアソン分布の大きな特徴は、平均と分散がともに \(\lambda\) に等しいという点です。
\[ E[Y]=\lambda,\qquad \mathrm{Var}(Y)=\lambda \]
つまり、平均が大きい（よく起こる）ほどばらつきも大きくなる、という件数データらしい性質を自然に表現しています。逆にいえば、実データで分散が平均を大きく上回る場合は「過分散（overdispersion）」が疑われ、後のパートで扱う対処が必要になります。
対数リンクを使うモデル \(\log E[Y_i]=\beta_0+\beta_1 x_i\) では、係数の解釈が明快です。説明変数 \(x\) が1単位増えると、期待件数は \(\exp(\beta_1)\) 倍になります。この \(\exp(\beta)\) を率比（rate ratio）と呼びます。例えば \(\exp(\beta_1)=1.5\) なら「\(x\) が1増えるごとに発生率が1.5倍になる」と読め、リスク比に似た直感的な解釈ができます。
さらに実務では、患者ごとに観察期間や曝露量が異なることがよくあります。1年追跡した人と半年しか追跡していない人では、同じ発生率でも件数が変わってしまうためです。そこで観察期間・曝露量 \(t_i\) を明示的に取り込むために、その対数をオフセット（offset）として加えます。
\[ \log E[Y_i]=\log(t_i)+\beta_0+\beta_1 x_i \]
この式を変形すると \(\log\big(E[Y_i]/t_i\big)=\beta_0+\beta_1 x_i\) となり、モデル化しているのが件数そのものではなく「単位期間あたりの発生率」であることがわかります。オフセットの係数は1に固定され、推定の対象にはなりません。こうして、観察期間の違いを公平に調整したうえで発生率を比較できるようになります。
なお、これらの係数 \(\beta\) は最尤推定（maximum likelihood estimation）によって求められ、Rでは glm() 関数の family = poisson オプションで簡単に実行できます。次のパートでは、実際のデータを用いてこの推定と率比の解釈を具体的に見ていきます。
Rでポアソン回帰を実装する
ここからは実際に手を動かしていきましょう。ある臨床観察研究で、患者ごとの有害事象（AE）件数を、治療群（group：control／active）と年齢（age）で説明するモデルを組み立てます。さらに患者ごとに観察期間（follow_years、人年）が異なるという現実的な設定にして、後半でオフセットを導入します。各ステップは「コード → 出力結果 → 📝解釈」の3点セットで進めます。
ステップ1：データを準備する
まずは set.seed() で再現性を確保しながら、架空データを生成します。
set.seed(1234)
n <- 200
dat <- data.frame(
id = 1:n,
group = factor(sample(c("control", "active"), n, replace = TRUE),
levels = c("control", "active")),
age = round(rnorm(n, mean = 60, sd = 10))
)
dat$follow_years <- round(runif(n, 0.5, 3.0), 2) # 患者ごとの観察期間（人年）
# 真の発生率にもとづいてAE件数を生成
lambda <- exp(-2.2 + 0.40 * (dat$group == "active") + 0.02 * dat$age) * dat$follow_years
dat$ae_count <- rpois(n, lambda)
head(dat)
> id group age follow_years ae_count
> 1 active 48 1.87 1
> 2 control 67 2.44 2
> 3 active 71 0.92 1
> 4 control 59 1.36 0
> 5 active 55 2.71 2
> 6 control 63 1.08 1
応答変数
ae_count は0, 1, 2…という非負の整数（件数）です。group は control を基準（reference）に設定しており、この順序が後の係数の符号を決めます。follow_years は0.5〜3.0人年とばらついている点に注目してください。観察期間が長い患者ほどAEは起こりやすいので、この差はステップ4のオフセットで調整します。ステップ2：基本のポアソン回帰を当てはめる
まずはオフセットなしで、件数そのものをモデル化します。
fit <- glm(ae_count ~ group + age, family = poisson(link = "log"), data = dat)
summary(fit)
> Coefficients:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -0.852400 0.321600 -2.651 0.00805 **
> groupactive 0.412300 0.108500 3.800 0.000144 ***
> age 0.020800 0.006400 3.250 0.001154 **
> ---
> (Dispersion parameter for poisson family taken to be 1)
> Null deviance: 245.67 on 199 degrees of freedom
> Residual deviance: 220.14 on 197 degrees of freedom
> AIC: 512.3
係数は対数スケールなので、指数変換して率比（rate ratio）に直します。
exp(coef(fit))
> (Intercept) groupactive age
> 0.4264000 1.5103000 1.0210000
active の係数 0.4123 を指数変換すると
exp(0.4123) = 1.510。つまり年齢を一定に保つと、active 群のAE発生は control 群の約1.51倍、率にして約51%高いと読めます。age の率比 exp(0.0208) = 1.021 は、年齢が1歳上がるごとにAEが約2.1%増えることを意味します。p値はいずれも0.01未満で、統計的に有意です。ステップ3：率比の95%信頼区間を出す
点推定だけでなく、区間で不確実性を示しましょう。係数と信頼区間をまとめて指数変換します。
exp(cbind(RateRatio = coef(fit), confint(fit)))
> Waiting for profiling to be done...
> RateRatio 2.5 % 97.5 %
> (Intercept) 0.4264000 0.2270214 0.8007800
> groupactive 1.5103000 1.2230000 1.8710000
> age 1.0210000 1.0083000 1.0339000
active の率比1.510に対する95%信頼区間は (1.223, 1.871) で、下限が1を上回っています。区間全体が1より大きいので「active でAEが多い」という結論は区間推定でも支持されます。
age の区間 (1.008, 1.034) も1を含まず、年齢の効果も安定しています。ステップ4：オフセットで人年あたり発生率をモデル化する
ここが件数データ解析の肝です。観察期間が患者ごとに違うのに件数をそのまま比べるのは不公平ですね。そこで offset(log(follow_years)) を加え、「人年あたりの発生率」をモデル化します。
fit_off <- glm(ae_count ~ group + age + offset(log(follow_years)),
family = poisson(link = "log"), data = dat)
summary(fit_off)$coefficients
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -2.147100 0.338800 -6.337986 2.330000e-10
> groupactive 0.381200 0.107900 3.533462 0.000410
> age 0.019100 0.006300 3.031746 0.002440
オフセット項は係数を1に固定した説明変数で、モデルは \(\log(E[\text{count}]/\text{follow\_years})=X\beta\)、すなわち発生率（率）そのものを推定します。active の率比は
exp(0.3812) = 1.464。観察期間の差を調整したうえで、active 群の人年あたりAE発生率は control の約1.46倍です。オフセットを入れないと「たまたま長く追跡された患者ほど件数が多い」という影響を治療効果と混同してしまうため、期間が不均一なときは必須の処理です。ステップ5：新しい患者の期待件数を予測する
最後に、当てはめたモデルで新規データの期待AE件数を予測します（ステップ2のオフセットなしモデルを使用）。
newdat <- data.frame(
group = factor(c("control", "active", "control", "active"),
levels = c("control", "active")),
age = c(50, 50, 70, 70)
)
predict(fit, newdata = newdat, type = "response")
> 1 2 3 4
> 1.206349 1.822021 1.828766 2.761988
type = "response" は対数スケールを逆変換し、期待件数（人数ではなく件数）を返します。50歳では control が1.21件、active が1.82件、70歳では control が1.83件、active が2.76件。同じ年齢どうしを比べると active／control の比はどこでも約1.51倍で、これはステップ2の率比と一致します。ポアソン回帰の効果が乗法的（率をかけ算で動かす）である様子がよく分かりますね。なお、当てはまりが良すぎ・悪すぎるときは過分散の可能性があり、これは後述します。
過分散への対処 ― quasiポアソンと負の二項回帰
ポアソン回帰は「平均＝分散」という強い前提を置いています。すなわち、応答変数 \(Y\) の条件付き分散が条件付き平均 \(\mu\) に等しい（\(\mathrm{Var}(Y)=\mu\)）と仮定します。しかし実データでは、この前提が崩れ、分散が平均を大きく上回ることがよくあります。これを過分散（overdispersion）と呼びます。個体差、観測されていない共変量、イベントの集積（クラスタ）などが原因で、件数データのばらつきが理論値より大きくなるのです。
過分散の診断はむずかしくありません。もっとも簡便なのは、Pearson χ²統計量を残差自由度で割った分散パラメータ（dispersion）が 1 を大きく超えていないかを見ることです。目安として 1.5 を超えたら要注意、2 以上なら明確な過分散と考えます。
# Pearson残差ベースの分散パラメータ
fit <- glm(count ~ group + offset(log(person_years)),
family = poisson, data = dat)
pearson <- residuals(fit, type = "pearson")
disp <- sum(pearson^2) / fit$df.residual
disp # 1 前後なら適合、大きく超えると過分散
summary(fit) の residual deviance が残差自由度（df）より大幅に大きいときも、過分散のサインです（理想的には両者はほぼ等しくなります）。
過分散を無視して通常のポアソン回帰を当てはめると、標準誤差が過小評価されます。その結果、信頼区間が実際より狭くなり、本来は有意でない係数を「有意」と誤判定してしまいます。件数データを扱うときは、係数を解釈する前に必ず分散パラメータを確認してください。
対処法は状況に応じて選びます。次の表に整理します。
| 状況 | 推奨モデル | Rの関数 |
|---|---|---|
| 分散＝平均が成り立つ | 通常のポアソン回帰 | glm(family=poisson) |
| 過分散だが係数の点推定は流用したい | quasiポアソン（標準誤差を分散パラメータで補正） | glm(family=quasipoisson) |
| 過分散が強く、ばらつき構造をモデル化したい | 負の二項回帰 | MASS::glm.nb() |
| ゼロが極端に多い（ゼロ過剰） | ゼロ過剰・ハードルモデル | pscl::zeroinfl()／hurdle() |
quasipoisson は係数の点推定値をポアソン回帰と同じにしたまま、標準誤差だけを分散パラメータ \(\sqrt{\phi}\) 倍して補正する簡便な手法です。一方、負の二項回帰は分散を \(\mu + \mu^2/\theta\) という形で明示的にモデル化し、過分散の大きさそのものをパラメータ \(\theta\) として推定します。過分散の背景に「個体ごとのばらつき」があると考えられる場合は、負の二項回帰のほうが自然な選択です。
library(MASS)
fit_nb <- glm.nb(count ~ group + offset(log(person_years)), data = dat)
summary(fit_nb) # theta が推定される
負の二項回帰の具体的な使いどころや再発イベントへの応用は、負の二項回帰と再発イベント率解析で詳しく解説しています。あわせてご覧ください。
なお、件数の多くが 0 に集中するゼロ過剰（zero-inflation）が疑われるときは、過分散とは別の問題です。この場合は pscl::zeroinfl() などのゼロ過剰モデルを検討します。
実務でのポイント
製薬・臨床の現場でポアソン回帰を使う場面は、想像以上に多くあります。もっとも典型的なのは有害事象（AE）発生率の群間比較です。観察期間が被験者ごとに異なるため、単純な件数の比較では公平になりません。人年（person-years）をオフセットに入れて発生率として比較することが必須になります。
・AE発生率の群間比較：件数そのものでなく、人年オフセットを入れた「発生率」で比較する。率比（rate ratio）が効果指標になる。
・再発イベント：喘息・COPDの増悪回数、入院回数など「繰り返し起こるイベント」の総数はポアソン系モデルの得意分野。
・過分散の確認を習慣に：分散パラメータを必ずチェックし、1を大きく超えるなら quasiポアソン／負の二項回帰へ切り替える。
・ゼロ過剰の見落とし注意：「発生ゼロ」の被験者が極端に多いデータではゼロ過剰モデルを検討する。
・相関構造があるなら発展モデルへ：同一被験者の反復測定や施設クラスタで観測が独立でないときは、GLMM（ランダム効果）やGEE（周辺モデル）へ拡張する。
とくに「過分散の確認」と「オフセットの適切な設定」は、解析結果の妥当性を左右する二大チェックポイントです。標準誤差の過小評価は、規制当局への申請資料でも見逃されやすい落とし穴ですので、レビューの段階で必ず確認する習慣をつけたいところです。
参考書籍
ポアソン回帰を含む一般化線形モデル（GLM）を体系的に学ぶために、定番となっている書籍を3冊紹介します。いずれも日本語で読め、Rによる実装や実データの例が豊富です。
関連記事
ポアソン回帰の理解を、隣接するモデルへと広げていくための記事を紹介します。
・負の二項回帰と再発イベント率解析：過分散が強い件数データの本命モデル。喘息・COPD増悪などの再発イベントを人年オフセットで解析する方法を解説します。
・GLMM（一般化線形混合モデル）とは：同一被験者の反復測定や施設クラスタで観測が独立でないとき、ランダム効果でポアソン回帰を拡張する道筋を示します。
・順序ロジスティック回帰（比例オッズモデル）とは：件数ではなく順序カテゴリを扱うGLM。応答変数の型に応じたモデル選択の視野を広げてくれます。
まとめ
ポアソン回帰は、件数や発生率といったカウントデータを扱うための一般化線形モデルの基本形です。対数リンクを通じて説明変数を線形に結びつけ、係数を指数変換すれば率比として解釈できる点が大きな魅力です。観察期間が異なるデータでは、人年などをオフセットに入れて「発生率」を比較することが欠かせません。
一方で、ポアソン回帰の中核前提である「平均＝分散」は現実のデータでしばしば崩れます。過分散を放置すると標準誤差が過小評価され、有意性の判断を誤ります。分散パラメータや residual deviance を確認し、必要に応じて quasiポアソンや負の二項回帰へ切り替える。この一手間が、解析の信頼性を大きく左右します。製薬・臨床では、AE発生率の群間比較や再発イベントの解析など、ポアソン系モデルの出番は数多くあります。反復測定やクラスタ相関があればGLMM・GEEへと発展させる視点も持っておきたいところです。
ポアソン回帰は、負の二項回帰・GLMM・ロジスティック回帰といった一般化線形モデル一族への入口でもあります。本記事で紹介した関連記事もあわせて読み進め、データの型と構造に応じてモデルを選び分ける力を養っていただければと思います。
