Stata - Poisson回帰

転倒、入院、感染のように「何回起きたか」を数えたデータを Stata で解析する手順です。R での同じ解析は [[R - Poisson回帰をやってみよう]] にあり、考え方の説明はそちらに書きました。介護施設の入所者を追跡し、睡眠薬の使用が転倒回数と関連するかを調べる設定を使います。 > [!WARNING] 出力について > この記事はコマンドと解釈の説明にとどめ、実行結果の出力を載せていません。手元で流して確認しながら読んでください。数値の読み方は [[R - Poisson回帰をやってみよう]] に実際の出力付きで書いてあります。 ## ステップ1：データを用意して分布を確認する 100人を1年間追跡し、転倒回数を数えたデータを作ります。 ```stata clear all set seed 16 set obs 100 gen id = _n gen sedative = (id > 50) label define sed 0 "睡眠薬なし" 1 "睡眠薬あり" label values sedative sed gen falls = rpoisson(0.4) if sedative == 0 replace falls = rpoisson(1.2) if sedative == 1 ``` `rpoisson()` がポアソン分布に従う乱数を返します。括弧の中が平均です。睡眠薬なしの50人は平均0.4回、ありの50人は平均1.2回で転ぶという設定にしました。転倒回数の散らばりを見ます。 ```stata tabulate falls histogram falls, discrete frequency ``` 0回が最も多く、回数が増えるほど人数が減る非対称な形になります。件数データの典型で、正規分布には見えません。次に平均と分散を比べます。ポアソン分布なら両者は等しくなるはずです（[[ポアソン分布とは]]）。 ```stata summarize falls display "平均 = " r(mean) display "分散 = " r(Var) ``` `summarize` の後は `r(mean)` と `r(Var)` に平均と分散が入っています。両者が近ければ、ポアソン分布を仮定してよいと判断します。観測された人数とポアソン分布から計算した人数を並べると、当てはまりが目で見て確かめられます。 ```stata summarize falls local lambda = r(mean) forvalues k = 0/4 { quietly count if falls == `k' local obs = r(N) local expected = 100 * poissonp(`lambda', `k') display "転倒 `k' 回 : 観測 `obs' 人 / 期待 " %5.1f `expected' " 人" } ``` `poissonp(m, k)` が、平均 m のポアソン分布で k 回起きる確率を返します。100倍すれば期待される人数になります。 ## ステップ2：説明変数を入れずに当てはめるまず説明変数を入れないモデルから始めます。`poisson` の後に従属変数だけを書くと、「全員が同じ λ を持つ」というモデルになります。ステップ1で確認したポアソン分布を、そのまま当てはめる操作にあたります。 ```stata poisson falls ``` 切片だけが推定されます。対数スケールなので、指数変換して戻します。 ```stata display exp(_b[_cons]) ``` この値は、ステップ1で計算した全体の平均と一致します。`poisson falls` は、データに最もよく合う λ を1つ探してくる操作であり、手で平均を計算するのと結果は変わりません。 `_b[_cons]` は直前の推定結果の切片を指します。Stata では推定後に `_b[変数名]` で係数を取り出せます。ここまでは100人全員が同じ転びやすさを持つと考えたことになります。この λ を人によって動かす仕組みが、説明変数です。 ``` 説明変数なし全員が同じ λ ↓ 説明変数あり λ が群によって変わる ``` ## ステップ3：群の違いを入れて発生率比を出す睡眠薬の使用で分けて、1年あたりの平均転倒回数を計算します。 ```stata tabstat falls, by(sedative) statistics(n sum mean variance) ``` 使用群のほうが平均が大きくなっているはずです。この差を評価するため、睡眠薬の使用を説明変数として加えます。 ```stata poisson falls i.sedative ``` `i.` は因子変数の指定で、0を基準として1の効果を推定します。出力される係数は対数スケールです。 `irr` オプションを付けると、係数を指数変換した値が表示されます。 ```stata poisson falls i.sedative, irr ``` R では `exp(coef())` と自分で変換しましたが、Stata はオプション1つで済みます。表示される値が発生率比（incidence rate ratio, IRR）で、信頼区間も指数変換された形で出ます。 IRR が2.5であれば、睡眠薬を使っている人は使っていない人の2.5倍の頻度で転倒しているという意味です。信頼区間が1をまたがなければ統計学的に有意と判断します。推定後に係数だけを確認したいときは次のように書きます。 ```stata display "IRR = " exp(_b[1.sedative]) ``` ## ステップ4：モデルが実際の分布に沿っているか確かめるモデルが出した予測が、実際のデータをどこまで再現できているかを見ます。`predict` で各人の期待件数を出し、群ごとに観測された人数と並べます。 ```stata quietly poisson falls i.sedative predict lambda, n forvalues g = 0/1 { display _newline "=== sedative = `g' ===" forvalues k = 0/4 { quietly count if falls == `k' & sedative == `g' local obs = r(N) quietly summarize lambda if sedative == `g', meanonly local expected = 50 * poissonp(r(mean), `k') display "転倒 `k' 回 : 観測 `obs' 人 / 期待 " %5.1f `expected' " 人" } } ``` `predict lambda, n` がモデルの予測した期待件数を返します。その値を持つポアソン分布で各回数が起きる確率を計算し、群の人数に直したものが期待人数です。観測と期待が近ければ、モデルがデータの形を再現できていると判断します。山の位置が大きく外れていたり、0回の人数だけが極端に多かったりする場合は、ポアソン分布の想定が合っていない可能性を疑います。 > [!NOTE] 全体でまとめて見ないこと > 100人をひとまとめにした度数分布で確認すると、当てはまりを見誤ります。全体の分布は平均さえ合っていればおおむね再現されてしまい、群による違いを入れたことの良し悪しが見えないためです。必ず群ごとに分けて確認します。実際の数値を使った比較は [[R - Poisson回帰をやってみよう]] のステップ4にあります。 ## ステップ5：追跡期間が人によって違うときここまでは全員をちょうど1年追跡していました。しかし実際の研究では、途中で入院した人、転出した人、亡くなった人がいて、追跡期間は人によって違います。追跡年数 `py` を持つデータを作ります。睡眠薬を使っている人ほど状態が悪く、早く追跡が終わる状況を想定します。 ```stata clear all set seed 75 set obs 100 gen id = _n gen sedative = (id > 50) label define sed 0 "睡眠薬なし" 1 "睡眠薬あり" label values sedative sed gen py = runiform(1.0, 2.5) if sedative == 0 replace py = runiform(0.2, 1.2) if sedative == 1 replace py = round(py, 0.1) gen falls = rpoisson(0.48 * py) if sedative == 0 replace falls = rpoisson(1.24 * py) if sedative == 1 ``` 群ごとに転倒の総数と追跡年数の総和を出します。 ```stata tabstat falls py, by(sedative) statistics(sum) bysort sedative: egen tot_falls = total(falls) bysort sedative: egen tot_py = total(py) gen rate = tot_falls / tot_py tabstat rate, by(sedative) ``` 転倒の総数は両群で近い値になりますが、追跡年数は睡眠薬なしの群が大きく上回ります。1年あたりに直すと差がはっきりします。追跡期間を無視して、ステップ3と同じように当てはめます。 ```stata poisson falls i.sedative, irr ``` IRR が1に近い値になり、「睡眠薬と転倒は無関係」という結論が出ます。これは誤りです。転倒の総数が近かったのは、使用群の追跡期間が短かったからであって、転びにくかったからではありません。短い期間しか見ていない群と長く見た群の件数をそのまま比べたために、本当にある差が消えてしまいました。 ## ステップ6：exposure を指定する追跡期間の違いを調整するには、`exposure()` オプションで観察期間の変数を指定します。 ```stata poisson falls i.sedative, exposure(py) irr ``` これで IRR に差が現れます。R で `offset(log(py))` と書いたものに対応します。 `exposure()` は指定した変数の対数を自動で計算してオフセットに入れます。対数を自分で取る必要はありません。自分で取る場合は `offset()` を使います。 ```stata gen logpy = log(py) poisson falls i.sedative, offset(logpy) irr ``` この2つは同じ結果になります。`exposure(py)` と `offset(py)` は別物です。後者は対数を取らずにオフセットへ入れてしまい、誤った推定値を返します。 > [!NOTE] R との対応 > R には `exposure()` にあたる書き方がなく、`offset(log(py))` と自分で対数を取ります。Stata の `exposure(py)` は、この `log()` を代わりにやってくれるオプションです。どちらを使っても推定結果は同じです。 `margins` で、追跡1年あたりの転倒回数を群ごとに求めます。 ```stata poisson falls i.sedative, exposure(py) margins sedative, at(py = 1) predict(n) ``` `predict(n)` が期待件数を指定し、`at(py = 1)` で「1年追跡した場合」に揃えます。ステップ5で手計算した1年あたりの率と一致するはずです。グラフにする場合は次のように続けます。 ```stata marginsplot ``` ## ステップ7：過分散を確認するポアソン回帰は平均と分散が等しいことを前提にしています。この前提が崩れていないかを見ます。 ```stata poisson falls i.sedative, exposure(py) estat gof ``` `estat gof` が逸脱度とピアソン統計量による適合度検定を返します。p値が小さければ、ポアソン分布の想定が合っていないサインです。過分散があるときは負の二項回帰に切り替えます。 ```stata nbreg falls i.sedative, exposure(py) irr ``` `nbreg` の出力の最後に、過分散パラメータ alpha がゼロかどうかの尤度比検定が表示されます。この検定が有意であれば、ポアソン回帰ではなく負の二項回帰を使うべきだと判断できます。標準誤差だけを補正する方法もあります。 ```stata poisson falls i.sedative, exposure(py) irr vce(robust) ``` 過分散があるのに素のポアソン回帰を使うと、標準誤差が小さく見積もられ、有意差が出やすくなります。 ## 補足：exposure と offset の使い分け | 書き方 | 対数 | 用途 | |--------|------|------| | `exposure(py)` | Stata が自動で取る | 観察期間の変数をそのまま渡す | | `offset(logpy)` | 自分で取っておく | 対数済みの変数を渡す、または特殊な調整をする | | `offset(py)` | 取らない | 誤り。率のモデルにならない | オフセット項がなぜ対数になるのか、対数を取ると負の値になる点をどう考えるかは、[[R - Poisson回帰をやってみよう]] の補足節に書きました。式の導出は言語によらず同じです。 ## 関連項目 - [[R - Poisson回帰をやってみよう]] - [[ポアソン分布とは]] - [[Stata - 生存時間解析]] - [[率・割合・比を区別して使う]]