yt coffee

Study hard, play harder.

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

  • 2020-01-18 01:54
  • #julia

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

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

import Random
Random.seed!(0)
nothing

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

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

using DataFrames, Distributions
sizes = rand(Normal(10.0, 0.7), 100)
data = DataFrame()
data.y = [rand(Poisson(3.5 + 0.5size)) for size in sizes]
data.x = sizes
data.f = vcat(["C" for _ in 1:50], ["T" for _ in 1:50])

describe(data)

3 rows × 8 columns

variablemeanminmedianmaxnuniquenmissingeltype
SymbolUnion…AnyUnion…AnyUnion…NothingDataType
1y8.5118.518Int64
2x9.982637.572849.9565611.6197Float64
3fCT2String

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

using Plots
scatter(data.x[data.f .== "C"], data.y[data.f .== "C"], label="C")
scatter!(data.x[data.f .== "T"], data.y[data.f .== "T"], label="T")
xlabel!("data.x")
ylabel!("data.y")
8910112.55.07.510.012.515.017.5data.xdata.yCT

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

箱ひげ図も描画します:

using StatsPlots
boxplot(data.f, data.y, legend=false, yticks=0:2:16)
CT246810121416

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

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

p(yiλi)=λiyiexp(λi)yi!p(y_i|\lambda_i) = \frac{\lambda_i^{y_i}\exp(-\lambda_i)}{y_i!}

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

λi=exp(β1+β2xi)\lambda_i = \exp(\beta_1 + \beta_2 x_i)

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

genlink = (β₁, β₂) -> (x) -> exp(β₁ + β₂ * x)
xs = -4:0.1:4
plot(xs, genlink(-1, 0.4), label="{-1, 0.4}")
xlabel!("x_i")
ylabel!("lambda_i")
plot!(xs, genlink(-2, -0.8), label="{-2, -0.8}")
-4-20240123x_ilambda_i{-1, 0.4}{-2, -0.8}

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

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

logλi=β1+β2xi\log \lambda_i = \beta_1 + \beta_2 x_i

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

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

logL(β1,β2)=ilogλiyiexp(λi)yi!\log L(\beta_1, \beta_2) = \sum_i \log \frac{\lambda_i^{y_i} \exp(-\lambda_i)}{y_i!}

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

using GLM
fit = glm(@formula(y ~ x), data, Poisson(), LogLink())
StatsModels.DataFrameRegressionModel{GeneralizedLinearModel{GlmResp{Array{Float64,1},Poisson{Float64},LogLink},DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Formula: y ~ 1 + x

Coefficients:
────────────────────────────────────────────────────
              Estimate  Std.Error  z value  Pr(>|z|)
────────────────────────────────────────────────────
(Intercept)  1.6692     0.470286   3.54932    0.0004
x            0.0472266  0.0468658  1.0077     0.3136
────────────────────────────────────────────────────

この出力内容は p.50 の summary(fit) を実行したときと対応しています。 β^1\hat{\beta}_1β^2\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() 関数を使います:

loglikelihood(fit)
-248.9474224949767

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

seq₁ = 1.6:0.01:1.7
seq₂ = 0.04:0.001:0.05
logL = (β₁, β₂) -> sum([logpdf(Poisson(exp(β₁ + β₂ * data.x[index])), y) for (index, y) in enumerate(data.y)])
plot(seq₁, seq₂, logL, st=:surface)

β̂₁ = 1.6692
β̂₂ = 0.0472266
maxL = -248.9474224949767
plot!([β̂₁], [β̂₂], [maxL], seriestype=:scatter, legend=false)
1.601.621.651.681.700.04000.04250.04500.04750.0500-256-254-252-250-257-256-255-254-253-252-251-250-249

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

scatter(data.x[data.f .== "C"], data.y[data.f .== "C"], label="C")
scatter!(data.x[data.f .== "T"], data.y[data.f .== "T"], label="T")
xlabel!("data.x")
ylabel!("data.y")
xx = [minimum(data.x), maximum(data.x)]
plot!(xx, ((x) -> exp(1.6692 + 0.0472266x)).(xx), label="prediction")
8910112.55.07.510.012.515.017.5data.xdata.yCTprediction

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

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

今回の fif_i は 2 水準ですが、仮に nn 水準あったとすると:

λi=exp(β1+jβjdij)\lambda_i = \exp\Biggl(\beta_1 + \sum_j \beta_j d_{ij}\Biggr)

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

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

つまり fif_i = C の場合は:

λi=exp(β1)\lambda_i = \exp(\beta_1)

fif_i = T の場合は:

λi=exp(β1+β3)\lambda_i = \exp(\beta_1 + \beta_3)

のようになります。(β2\beta_2xix_i に使っているので、あとの説明との整合性のために 3 にしています)

fit2 = glm(@formula(y ~ f), data, Poisson(), LogLink())
StatsModels.DataFrameRegressionModel{GeneralizedLinearModel{GlmResp{Array{Float64,1},Poisson{Float64},LogLink},DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Formula: y ~ 1 + f

Coefficients:
───────────────────────────────────────────────────────
               Estimate  Std.Error    z value  Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept)   2.15409    0.0481682  44.7201      <1e-99
f: T         -0.0258534  0.0685648  -0.377065    0.7061
───────────────────────────────────────────────────────

β^1\hat{\beta}_1β^3\hat{\beta}_3 はそれぞれ 2.152.150.026-0.026 と推定されています。f: T と書かれていることからここでいう β3\beta_3fif_i が T のときの傾きであると分かります。

loglikelihood(fit2)
-249.38600696156576

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

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

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

fit3 = glm(@formula(y ~ x + f), data, Poisson(), LogLink())
StatsModels.DataFrameRegressionModel{GeneralizedLinearModel{GlmResp{Array{Float64,1},Poisson{Float64},LogLink},DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Formula: y ~ 1 + x + f

Coefficients:
───────────────────────────────────────────────────────
               Estimate  Std.Error    z value  Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept)   1.69225    0.479588    3.52854     0.0004
x             0.0457861  0.0472654   0.968701    0.3327
f: T         -0.0174208  0.0690733  -0.252207    0.8009
───────────────────────────────────────────────────────
loglikelihood(fit3)
-248.9156152929941

まとめると

λi={exp(1.69+0.05xi)(fi=C)exp(1.69+0.05xi0.02)(fi=T)\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\exp(-0.02) = 0.98 なので肥料をやると平均が 0.98 倍になると予測されています。