一般化線形モデル - ポアソン回帰

「データ解析のための統計モデリング入門」の第3章です。前回の確率分布と統計モデルの最尤推定の続きです。

まずは準備として、再現性のために乱数の種を設定します：

第3章では、架空植物 100 個体を調査して得られた、個体ごとの種子数のデータの分析をします。植物個体 $i$ の種子数は $y_i$ 個であり、また個体の属性の一つである体サイズ $x_i$ が観測されています。さらに、全体のうち 50 個体（$i \in {1, 2, \dots, 50}$）は何も処理をしていないけれど、残り 50 個体（$i \in {51, 52, \dots, 100 }$）には肥料を加える処理をほどこしています。

今回もダミーのデータを作ります。第3章のデータを作り方の説明が本文中にないのですが、ポアソン回帰を題材としているので、応答変数 $y_i$ はポアソン分布を使っていると想定し、説明変数 $x_i$ は p.44 のデータの概要からそれっぽい正規分布を使って作ってみました：

得られたデータを可視化してみます：

（データの作り方から当然ではありますが） $x_i$ が大きくなるに連れて $y_i$ が大きくなる傾向が見て取れます。それと比較すると $f_i$ は影響があるのか無いのか微妙な感じがします。

箱ひげ図も描画します：

ポアソン回帰の統計モデル

前回の記事でカウントデータはポアソン分布を使って表現できることを確認しました。今回のデータも種子数を数えたカウントデータ、ということになっているので、ポアソン分布で表現できそうです。ある個体 $i$ において種子数が $y_i$ である確率 $p(y_i|\lambda_i)$ はポアソン分布に従っていると仮定すると：

$$ p(y_i|\lambda_i) = \frac{\lambda_i^{y_i}\exp(-\lambda_i)}{y_i!}$$

散布図を眺めた感じから施肥効果 $f_i$ はあまり種子数に効果がなさそうなので、体サイズ $x_i$ だけに依存する統計モデルについて考えてみます。この場合平均種子数 $\lambda_i$ を説明変数 $x_i$ の関数として定義すればいいので、詳しい種明かしはあとでするとしてひとまず以下のようにおいてみます：

$$ \lambda_i = \exp(\beta_1 + \beta_2 x_i) $$

いくつか $\beta_1$ と $\beta_2$ を変えて、この関係を図示すると:

$\lambda_i$ が非負の値になっていることが分かります。 ポアソン分布の平均 $\lambda_i$ は必ず非負でならなければなりませんが、これであれば説明変数 $x_i$ がどのような値でも自動的にこの条件が守られて都合がよさそうです。

ところで、この式は以下のように変形できます：

$$ \log \lambda_i = \beta_1 + \beta_2 x_i $$

このときの右辺 $\beta_1 + \beta_2 x_i$ は ${\beta_1, \beta_2}$ の線形結合になっていることから、 線形予測子（linear predictor）と呼ばれます。またこの式のように（応答変数の期待値）=（線形予測子）となっているものを リンク関数（link function）と呼びます。今回の場合は対数関数なので 対数リンク関数（log link function）です。リンク関数はモデルごとによく使われるものがある程度決まっていて、前述の理由からポアソン回帰ではたいていは対数リンク関数が使われます。

ポアソン回帰はポアソン分布を使った統計モデルの対数尤度 $\log L$ が最大になるパラメータ $\hat{\beta}_1$ と $\hat{\beta}_2$ の推定値を求めます。この対数尤度は：

$$ \log L(\beta_1, \beta_2) = \sum_i \log \frac{\lambda_i^{y_i} \exp(-\lambda_i)}{y_i!}$$

パラメータが 2 つあるので解析的に答えを求めるのは簡単ではありませんが、 Julia の GLM パッケージを使うと切片 $\beta_1$ と傾き $\beta_2$ の最尤推定値を簡単に求めることができます。

この出力内容は p.50 の summary(fit) を実行したときと対応しています。 $\hat{\beta}_1$ と $\hat{\beta}_2$ はそれぞれ 1.67 と 0.05 と推定されています。 Std.Error は 標準誤差（standard error, SE）、 z value は最尤推定値を SE で除した z 値、あるいは Wald 統計量（Wald statistics）です。詳しくは書籍を参照してください。

@formula は Julia の StatsModels package にあるマクロです。本に書かれた R コードと記述の方法が同じなので、おそらく R から移植されたものなのだと思いますが、y とか x という変数は存在しないのに、いきなり第一引数で使うのがなかなか違和感があります：

fit <- glm(y ~ x, data = d, familly = poisson)

話を戻して、得られたモデルの最大対数尤度を評価するには StatsBase パッケージの loglikelihood() 関数を使います：

最尤推定値とその周辺の様子とプロットしてみると、たしかにこの点で最大値をとっていることが分かります：

得られたポアソン回帰の推定結果を使って、さまざまな $x_i$ における平均種子数 $\lambda$ の予測をしてみると、若干の右肩上がりのグラフが得られました。

説明変数が因子型の統計モデル

ここまで使ってきた説明変数 $x_i$ は数値でした。次は因子型データである $f_i$ です。この場合も使い方は同じで、そのまま @formula 関数で指定します。

今回の $f_i$ は 2 水準ですが、仮に $n$ 水準あったとすると：

$$ \lambda_i = \exp\Biggl(\beta_1 + \sum_j \beta_j d_{ij}\Biggr) $$

このとき $d_{ij}$ は本の表現を借りるならダミー変数で $f_i$ が j 番目の因子だったときだけ 1 になります。たとえば今回は C と T の 2 水準なので：

$$ d_{ij} = \begin{cases}0 & (f_i = \text{C}) \ 1 & (f_i = \text{T}) \end{cases}$$

つまり $f_i$ = C の場合は:

$$ \lambda_i = \exp(\beta_1) $$

$f_i$ = T の場合は：

$$ \lambda_i = \exp(\beta_1 + \beta_3) $$

のようになります。（$\beta_2$ は $x_i$ に使っているので、あとの説明との整合性のために 3 にしています）

$\hat{\beta}_1$ と $\hat{\beta}_3$ はそれぞれ $2.15$、 $-0.026$ と推定されています。f: T と書かれていることからここでいう $\beta_3$ は $f_i$ が T のときの傾きであると分かります。

このモデルでの最大対数尤度は先程のモデルより小さくなりました。

説明変数が数量型＋因子型の統計モデル

$x_i$ と $f_i$ 両方を組み合わせた統計モデルも同様に求めることができます。

まとめると

$$ \lambda_i = \begin{cases}\exp(1.69 + 0.05 x_i) & (f_i = \text{C}) \ \exp(1.69 + 0.05 x_i - 0.02) & (f_i = \text{T}) \end{cases}$$

この式は （平均）=（定数）×（サイズの効果）×（施肥処理の効果）と分解できます。 $\exp(-0.02) = 0.98$ なので肥料をやると平均が 0.98 倍になると予測されています。