R - Poisson回帰をやってみよう

転倒、入院、感染のように「何回起きたか」を数えたデータを解析するのが Poisson 回帰です。追跡期間が人によって違うときは、オフセット項を加えて調整します。ここでは介護施設の入所者を追跡し、睡眠薬の使用が転倒回数と関連するかを調べる設定を使います。 ## ステップ1：データを用意して分布を確認する 100人を1年間追跡し、その間の転倒回数を数えたデータを作ります。乱数で作るので、`set.seed()` を同じにすれば手元でも同じ結果になります。 ```r set.seed(16) fall <- data.frame( id = 1:100, sedative = factor(rep(c("no", "yes"), each = 50), levels = c("no", "yes")), falls = c(rpois(50, 0.4), rpois(50, 1.2)) ) head(fall) ``` `sedative` が睡眠薬の使用の有無、`falls` が1年間の転倒回数です。100人全員をちょうど1年追跡した状況を考えます。`sedative` はステップ3まで使わず、しばらくは `falls` だけを見ていきます。転倒回数がどう散らばっているかを見ます。 ```r table(fall$falls) ``` ``` 0 1 2 3 4 41 40 13 4 2 ``` 0回の人が最も多く、回数が増えるほど人数が減る非対称な形です。件数データの典型的な形で、正規分布には見えません。次に平均と分散を比べます。ポアソン分布であれば両者は等しくなるはずです（[[ポアソン分布とは]]）。 ```r mean(fall$falls) var(fall$falls) ``` ``` [1] 0.86 [1] 0.8690909 ``` 平均0.86に対して分散0.869で、ほぼ一致しています。この先はポアソン分布を仮定して進めます。観測された人数とポアソン分布から計算した人数を並べます。 ```r lambda <- mean(fall$falls) data.frame( 転倒回数 = 0:4, 観測人数 = as.numeric(table(factor(fall$falls, levels = 0:4))), 期待人数 = round(dpois(0:4, lambda) * 100, 1) ) ``` ``` 転倒回数観測人数期待人数 1 0 41 42.3 2 1 40 36.4 3 2 13 15.6 4 3 4 4.5 5 4 2 1.0 ``` ![[poisson_dist_20260925.png|600]] 棒が実際の人数、赤い線がポアソン分布から計算した人数です。よく重なっています。 ## ステップ2：説明変数を入れずに当てはめるまず説明変数を入れないモデルから始めます。`~` の右に `1` とだけ書くと、「全員が同じ λ を持つ」というモデルになります。これはステップ1で確認したポアソン分布を、そのまま当てはめる操作にあたります。 ```r m0 <- glm(falls ~ 1, data = fall, family = poisson(link = "log")) summary(m0) ``` ``` Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.1508 0.1078 -1.399 0.162 ``` 切片だけが推定されました。対数スケールなので `exp()` で戻します。 ```r exp(coef(m0)) ``` ``` (Intercept) 0.86 ``` 0.86 という値は、ステップ1で計算した全体の平均と同じです。 ```r mean(fall$falls) ``` ``` [1] 0.86 ``` つまり `glm(falls ~ 1, family = poisson)` は、データに最もよく合う λ を1つ探してくる操作です。手で平均を計算するのと結果は変わりません。推定された λ から期待人数を計算すると、ステップ1の表の右列と一致します。 ```r round(dpois(0:4, exp(coef(m0))) * 100, 1) ``` ``` [1] 42.3 36.4 15.6 4.5 1.0 ``` ここまでは、100人全員が同じ転びやすさを持つと考えたことになります。しかし実際には、睡眠薬を使っている人とそうでない人で転びやすさが違うかもしれません。この λ を人によって動かす仕組みが、説明変数です。 ``` 説明変数なし全員が同じ λ λ = 0.86 ↓ 説明変数あり λ が人によって変わる λ = 0.48（睡眠薬なし） λ = 1.24（睡眠薬あり） ``` ## ステップ3：群の違いを入れて発生率比を出す睡眠薬の使用で分けて、1年あたりの平均転倒回数を計算します。 ```r aggregate(falls ~ sedative, data = fall, FUN = mean) ``` ``` sedative falls 1 no 0.48 2 yes 1.24 ``` 使用群は平均1.24回、非使用群は0.48回です。この差を評価するため、睡眠薬の使用を説明変数として加えます。書き方はステップ2と同じで、`~` の右を `1` から `sedative` に変えるだけです。 ```r m1 <- glm(falls ~ sedative, data = fall, family = poisson(link = "log")) summary(m1) ``` ``` Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7340 0.2041 -3.596 0.000323 *** sedativeyes 0.9491 0.2404 3.948 7.89e-05 *** ``` 係数は対数スケールの値なので、このままでは読めません。`exp()` で戻すと解釈できる値になります。 ```r exp(coef(m1)) ``` ``` (Intercept) sedativeyes 0.480000 2.583333 ``` この2つの数字は、先ほど計算した群ごとの平均から導けます。 - `exp(切片)` = 0.48 は、非使用群の1年あたり平均転倒回数そのものです - `exp(切片 + 係数)` = exp(-0.7340 + 0.9491) = 1.24 で、使用群の平均と一致します - `exp(係数)` = 2.58 は、1.24 ÷ 0.48 にあたる比ですこの比を発生率比（incidence rate ratio, IRR）と呼びます。睡眠薬を使っている人は、使っていない人の2.58倍の頻度で転倒しているという意味です。信頼区間も指数変換します。 ```r round(cbind(IRR = exp(coef(m1)), exp(confint.default(m1))), 3) ``` ``` IRR 2.5 % 97.5 % (Intercept) 0.480 0.322 0.716 sedativeyes 2.583 1.613 4.138 ``` IRR 2.58（95%CI 1.61-4.14）と報告します。区間が1をまたがないので、統計学的に有意です。 ## ステップ4：モデルが実際の分布に沿っているか確かめるモデルが出した予測が、実際のデータをどこまで再現できているかを見ます。群ごとに、観測された人数とモデルが予測した人数を並べます。 ```r k <- 0:4 for (g in c("no", "yes")) { sub <- fall[fall$sedative == g, ] lam <- unique(fitted(m1)[fall$sedative == g]) # その群の λ print(data.frame( 群 = g, 転倒回数 = k, 観測 = as.numeric(table(factor(sub$falls, levels = k))), 期待 = round(dpois(k, lam) * 50, 1) )) } ``` `fitted(m1)` がモデルの予測した λ で、睡眠薬なしの群は0.48、ありの群は1.24でした。その λ を持つポアソン分布で各回数が起きる確率を計算し、50人分に直したものが期待人数です。 ``` 群転倒回数観測期待 1 no 0 29 30.9 2 no 1 18 14.9 3 no 2 3 3.6 4 no 3 0 0.6 5 no 4 0 0.1 群転倒回数観測期待 1 yes 0 12 14.5 2 yes 1 22 17.9 3 yes 2 10 11.1 4 yes 3 4 4.6 5 yes 4 2 1.4 ``` ![[poisson_fit_20260925.png|700]] 棒が実際の人数、赤い点と線がモデルの予測です。睡眠薬なしの群は0回に、ありの群は1回に山があり、モデルもそのとおりに予測しています。裾の落ち方も再現できています。 1回の人数が両群とも予測より多めですが、50人ずつの標本ではこの程度のずれは通常の範囲です。山の位置が大きく外れていたり、0回の人数だけが極端に多かったりする場合は、ポアソン分布の想定が合っていない可能性を疑います。 > [!NOTE] 全体でまとめて見ないこと > 100人をひとまとめにした度数分布で確認すると、当てはまりを見誤ります。全体の分布は平均さえ合っていればおおむね再現されてしまい、群による違いを入れたことの良し悪しが見えないためです。必ず群ごとに分けて確認します。 ## ステップ5：追跡期間が人によって違うときここまでは全員をちょうど1年追跡していました。しかし実際の研究では、途中で入院した人、転出した人、亡くなった人がいて、追跡期間は人によって違います。追跡年数 `py` を持つデータを作ります。睡眠薬を使っている人ほど状態が悪く、早く入院や転院になって追跡が短く終わる状況を想定します。 ```r set.seed(75) py <- c(round(runif(50, 1.0, 2.5), 1), round(runif(50, 0.2, 1.2), 1)) fall2 <- data.frame( id = 1:100, sedative = factor(rep(c("no", "yes"), each = 50), levels = c("no", "yes")), py = py, falls = c(rpois(50, 0.48 * py[1:50]), rpois(50, 1.24 * py[51:100])) ) head(fall2) ``` ``` id sedative py falls 1 1 no 1.3 0 2 2 no 1.0 0 3 3 no 2.0 0 4 4 no 1.5 3 5 5 no 1.1 1 6 6 no 1.2 1 ``` 群ごとに転倒の総数と追跡年数の総和を出します。 ```r agg <- aggregate(cbind(falls, py) ~ sedative, data = fall2, FUN = sum) agg$rate <- round(agg$falls / agg$py, 3) agg ``` ``` sedative falls py rate 1 no 42 86.7 0.484 2 yes 42 33.9 1.239 ``` 転倒の総数はどちらの群も42回で、まったく同じです。しかし追跡年数は非使用群が86.7年、使用群が33.9年と、2.5倍以上の開きがあります。1年あたりに直すと0.484回と1.239回で、大きく違います。追跡期間を無視して、ステップ3と同じように当てはめます。 ```r m_wrong <- glm(falls ~ sedative, data = fall2, family = poisson) round(summary(m_wrong)$coefficients, 4) ``` ``` Estimate Std. Error z value Pr(>|z|) (Intercept) -0.1744 0.1543 -1.1299 0.2585 sedativeyes 0.0000 0.2182 0.0000 1.0000 ``` 係数がちょうど0、p値が1.000です。IRR にすると1.00で、「睡眠薬と転倒は無関係」という結論になります。これは誤りです。転倒の総数が同じだったのは、使用群の追跡期間が短かったからであって、転びにくかったからではありません。短い期間しか見ていない群と長く見た群の件数をそのまま比べたために、本当にある差が消えてしまいました。 ## ステップ6：オフセット項を入れる追跡期間の違いを調整するには、`offset(log(追跡期間))` をモデルに加えます。 ```r m2 <- glm(falls ~ sedative + offset(log(py)), data = fall2, family = poisson) summary(m2) ``` ``` Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7248 0.1543 -4.697 2.64e-06 *** sedativeyes 0.9390 0.2182 4.303 1.68e-05 *** ``` ```r round(cbind(IRR = exp(coef(m2)), exp(confint.default(m2))), 3) ``` ``` IRR 2.5 % 97.5 % (Intercept) 0.484 0.358 0.656 sedativeyes 2.558 1.668 3.923 ``` IRR 2.56（95%CI 1.67-3.92）で、有意な関連が現れました。ステップ3で全員を1年追跡していたときの2.58とほぼ同じ値です。 ![[poisson_offset_20260925.png|600]] 横軸が追跡年数、縦軸が転倒回数です。赤い点（睡眠薬あり）は左側、つまり追跡の短い側に固まっています。直線の傾きが1年あたりの転倒回数にあたり、赤のほうが急です。件数を縦軸だけで比べると同じに見えても、傾きで比べると差があることが読み取れます。オフセットを入れたモデルでは、`py = 1` を与えると「1人を1年追跡したときの転倒回数」が得られます。 ```r nd <- data.frame( sedative = factor(c("no", "yes"), levels = c("no", "yes")), py = 1 ) nd$pred <- predict(m2, newdata = nd, type = "response") nd ``` ``` sedative py pred 1 no 1 0.4844291 2 yes 1 1.2389381 ``` ステップ5で手計算した1年あたりの率と一致します。 ## ステップ7：過分散を確認するポアソン回帰は平均と分散が等しいことを前提にしています。この前提が崩れていないかを見ます。 ```r m2$deviance / m2$df.residual sum(residuals(m2, type = "pearson")^2) / m2$df.residual ``` ``` [1] 1.10224 [1] 0.9584089 ``` どちらも1に近いので問題ありません。目安として2を超えるようなら過分散を疑い、負の二項回帰に切り替えるか、標準誤差を補正します。 ```r # 標準誤差だけを補正する方法 m_qp <- glm(falls ~ sedative + offset(log(py)), data = fall2, family = quasipoisson) # 負の二項回帰に切り替える方法 library(MASS) m_nb <- glm.nb(falls ~ sedative + offset(log(py)), data = fall2) ``` 過分散があるのに素のポアソン回帰を使うと、標準誤差が小さく見積もられ、有意差が出やすくなります。 ## 補足：オフセット項の中身率のモデルを立てて式を整理すると、`offset(log(py))` という形が出てきます。使うだけなら読み飛ばして構いません。 ### なぜ対数をとるのか Poisson 回帰では、期待件数 μ を直接モデル化せず、その対数をモデル化します。 $ \log \mu_i = \beta_0 + \beta_1 x_i $ 件数は負にならないので、右辺がどんな値をとっても μ = exp(右辺) は正に保たれます。係数が差ではなく比として解釈できるようにもなります。使用群の期待件数を μ₁、非使用群を μ₀ として差をとると、次のようになります。 $ \log \mu_1 - \log \mu_0 = \log \frac{\mu_1}{\mu_0} = \beta_1 $ 両辺を指数変換すれば、exp(β₁) が件数の比そのものになります。ステップ3で `exp(係数)` が群平均の比と一致したのは、この式のためです。 ### オフセット項が出てくる筋道追跡期間が違うとき、比べたいのは件数そのものではなく、単位時間あたりの率です。観察期間を t とすると、率は μ / t になります。モデル化したいのはこちらです。 $ \log \frac{\mu_i}{t_i} = \beta_0 + \beta_1 x_i $ 左辺の対数を分解します。 $ \log \mu_i - \log t_i = \beta_0 + \beta_1 x_i $ log t を右辺へ移します。 $ \log \mu_i = \beta_0 + \beta_1 x_i + \log t_i $ 最後の `log t` が、係数が1に固定された説明変数の形で現れました。これがオフセット項です。ふつうの説明変数は係数をデータから推定しますが、オフセットは推定せず1に固定します。式の導出から必然的にそうなるためです。 `offset()` を使わず、次のように書いても同じ意味になりません。 ```r # これは間違い。log(py) の係数を推定してしまう m_notoffset <- glm(falls ~ sedative + log(py), data = fall2, family = poisson) round(coef(m_notoffset), 4) ``` ``` (Intercept) sedativeyes log(py) -0.7124 0.9200 0.9787 ``` log(py) に0.9787という係数が推定されています。オフセットとして入れれば、ここは推定されず1に固定されます。今回は偶然1に近い値になったため IRR も2.51とあまり変わりませんが、これはデータがそうだっただけです。追跡期間と件数の関係が比例からずれているデータでは係数が1から離れ、率のモデルとしての解釈が失われます。 ### 対数をとると負になってしまう追跡年数が1年未満の人では、log をとると負の値になります。今回のデータでも100人中40人が1年未満でした。 ```r sum(fall2$py < 1) ``` ``` [1] 40 ``` ```r log(0.2) ``` ``` [1] -1.609438 ``` 負の値が出ても問題はありません。オフセットは線形予測子に足されるだけで、最終的な期待件数は次のように計算されます。 $ \mu_i = \exp(\beta_0 + \beta_1 x_i + \log t_i) = t_i \times \exp(\beta_0 + \beta_1 x_i) $ 指数関数の外へ出すと、率 exp(β₀ + β₁x) に観察期間 t を掛けた形に戻ります。t が1未満なら、期待件数はそれだけ小さくなります。対数が負であることと、期待件数が負になることは別の話です。期待件数が負にならないのは、exp() が常に正の値を返すためです。 ### 時間の単位を変えたらどうなるか追跡年数を月に直して当てはめてみます。こうすると log は正の値になります。 ```r fall2$pm <- fall2$py * 12 m_month <- glm(falls ~ sedative + offset(log(pm)), data = fall2, family = poisson) c(年単位 = coef(m2)[1], 月単位 = coef(m_month)[1]) c(年単位 = exp(coef(m2))[2], 月単位 = exp(coef(m_month))[2]) ``` ``` 年単位.(Intercept) 月単位.(Intercept) -0.7247843 -3.2096909 年単位.sedativeyes 月単位.sedativeyes 2.557522 2.557522 ``` IRR は2.5575のまま変わりません。変わるのは切片だけで、その差は -0.7248 -(-3.2097) = 2.4849 で、これは log(12) と一致します。切片は「1単位の期間あたりの率」を表すので、単位を月に変えれば12分の1になります。一方 IRR は比なので、単位の取り方に影響されません。どの単位を使っても結論は同じですが、切片や予測値を報告するときは、何単位あたりの率なのかを明記する必要があります。 ### log を付け忘れるとオフセットに log を付け忘れるのはよくある間違いです。 ```r m_bad <- glm(falls ~ sedative + offset(py), data = fall2, family = poisson) exp(coef(m_bad))[2] ``` ``` sedativeyes 3.043515 ``` エラーも警告も出ず、3.04という一見もっともらしい値が返ってきます。正しい値2.56とは違いますが、この数字だけを見て誤りに気づくのは困難です。オフセットには必ず log をとった値を入れます。 > [!NOTE] 観察期間が同じなら不要か > 全員の観察期間が同じなら、オフセットは定数になり切片に吸収されるので、入れても入れなくても IRR は変わりません。ステップ3でオフセットを使わなかったのはこのためです。ただし入れておけば切片が「1年あたりの率」として読めるようになる利点があります。 ## 関連項目 - [[ポアソン分布とは]] - [[Stata - Poisson回帰]] - [[率・割合・比を区別して使う]] - [[R - 生存時間解析の基礎]]