[mathjax]
はじめに
前回の記事では、線形回帰に続き、ロジスティック回帰を見ていきました。
目的変数である(y)は、線形回帰においては任意の実数値、ロジスティック回帰においては2値をとりました。では、目的変数(y)が、例えば植物の種の個数など非負整数をとる場合はどのように回帰モデルを設計したら良いでしょうか。結論としては、(y)がポアソン分布に従うとしてモデルを組むことができます。このように、(y)が正規分布に限らず様々な確率分布に従う回帰モデルを統一的に扱うための枠組みとして、一般化線形モデル（generalized linear model, GLM）が存在します。今回は、GLMについて見ていきます。
一般化線形モデル（GLM）を構成する3要素
1) 指数型分布族（exponential family）に従う目的変数 (p(y; theta, phi))
2) 線形予測子（linear predictor）(eta = beta^mathrm{T}x)
3) リンク関数（link function）(g(mu)=eta=beta^mathrm{T}x). ただし、(mu=mathbb{E}[y|x])
なお、これまでのブログでは線形モデルのパラメタとして(theta)を用いてきましたが、GLMの一般的な変数表記に合わせて、指数型分布族におけるパラメタを(theta)とし、代わりにモデルパラメタを(beta)と表記します。
一般化線形モデルでは、目的変数(y)は全て指数型分布族に属する確率分布によりモデリングされます。回帰のゴールは、(y)の期待値である(mu=mathbb{E}[y|x])を求めることです。上の3つの要素について、一つずつ確認していきましょう。
補足
前回までの記事では、(y)の期待値は(h_theta(x))という表記をしています。今回はモデルパラメタが(beta)ですから、(y)の期待値は(h_beta(x))であるとも言えます。
1) 指数型分布族 exponential family
指数型分布族は、次のような数式で表される確率分布のグループのことです。
$$
large
p(y; theta, phi) = exp bigr(frac{ytheta – b(theta)}{a(phi)}+c(y, phi) bigl).
$$
(theta, phi)はパラメタ、(a, b, c)は関数です。
指数型分布族は、一般化の度合いや文脈により表式にはいくつかの種類があるため、注意してください。(theta)は上記のような表式の場合には、natural parameter（自然パラメタ）またはcanonical parameter（正準パラメタ）と呼ばれます。(phi)はdispersion parameter（分散パラメタ）と呼び、確率分布の分散に対応しますが、GLMにおいては分散は最適化に関わらないため省略することがあります。また、関数(c)の部分を(exp)の外に出して表記することがあります。
さて、いくつかの主要な分布（正規分布、ベルヌーイ分布、ポアソン分布）が指数型分布族に属することを確認していきます。
正規分布
正規分布の確率密度関数は次のように与えられます。
$$
large
begin{align}
p(y) &= frac{1}{sqrt{2pi}sigma}exp bigl(- frac{(y-mu)^2}{2sigma^2} bigr) \
&= exp bigl(log frac{1}{sqrt{2pi}sigma} – frac{(y-mu)^2}{2sigma^2} bigr) \
&= exp bigl( frac{ymu – mu^2/2}{sigma^2} – y^2/(2sigma^2) – logsigma – log sqrt{2pi} bigr)
end{align}
$$
natural parameterは(theta=mu)、dispersion parameterは(phi=sigma)です。
(a(phi) = sigma^2, b(theta) = mu^2/2=theta^2/2, c(y, phi) = – y^2/(2sigma^2) – logsigma – log sqrt{2pi})です。
ベルヌーイ分布
ベルヌーイ分布は、成功確率(p)を用いて次のように表されます。
$$
large
begin{align}
p(y) &= p^y (1-p)^{1-y} \
&= exp (y log p + (1-y) log(1-p)) \
&= exp (y log p + log (1-p) – y log(1-p)) \
&= exp (y log (frac{p}{1-p}) + log(1-p))
end{align}
$$
(y)の期待値(mathbb{E}[y]=mu)は(p)であることに注意します。natural parameterは(theta = log (frac{p}{1-p}) = log (frac{mu}{1-mu}) )、dispersion parameterはありません。natural parameterと(theta)の関係から、(p=e^theta/(1+e^theta))であることが分かります。(a(phi)=1, b(theta) = log(1-p)= – log(1+e^theta), c(y, phi) = 0)です。
ポアソン分布
ポアソン分布は確率変数が非負整数をとる、すなわち(yin {0,1,2,3, ldots })である確率分布です。期待値を(mu)として、次のようになります。
$$
large
begin{align}
p(y) &= frac{1}{y!} e^{-mu}mu^y \
&= exp bigl(ylog mu – mu – log (y!) bigr)
end{align}
$$
natural parameterは(theta=log mu)、dispersion parameterはありません。(a(phi)=1, b(theta)=mu=e^theta, c(y, phi) = – log(y!))です。
2) 線形予測子 linear predictor
線形予測子(eta)は、モデルパラメタ(beta=(beta_0, beta_1, beta_2, ldots, beta_n)^mathrm{T})を係数とする、説明変数のベクトル(x=(x_0, x_1, x_2, ldots, x_n)^mathrm{T})の線形結合(beta^mathrm{T}x)です。(eta)をnatural parameterと表記する資料もあるため、注意してください（それはnatural parameterと線形予測子が等しくなるケースであり、次のセクション 3)で見ていきます）。今回の表式では、(eta)は説明変数の線形結合を「ひとまとまり」に扱う変数に過ぎず、それ以上の意味は持ちません。
セクション 1)では、正規分布・ベルヌーイ分布・ポアソン分布が指数型分布族に属し、natural parameterである(theta)によって分布を決定できることを確認しました。ということは、natural parameterを説明変数(x)で表すことができれば、モデルを構築することができます。一般化線形モデルでは、natural parameterを線形予測子 (eta = beta^mathrm{T}x)によりモデリングします（これが一般化線形モデルの「線形」の所以です）。
3) リンク関数 link function
回帰のゴールは、説明変数(y)の期待値(mu=mathbb{E}[y|x])を求めることです。リンク関数は、(mu)と線形予測子(eta)を「接続 (link)」する次のような関数のことです。
$$
large
g(mu) = eta
$$
リンク関数は、一般化線形モデルのデザインにより様々な選択肢があります。このうち、線形回帰やロジスティック回帰で用いてきたリンク関数は、正準リンク関数（canonical link function）と呼ばれるものです。正準リンク関数は、natural parameterを線形予測子を用いて次のようにモデリングすることで、分布に対して一意に決まります。
$$
large
theta = eta
$$
線形回帰（正規分布）、ロジスティック回帰（ベルヌーイ分布）、ポアソン回帰（ポアソン分布）について正準リンク関数を求めていきましょう。
正規分布
$$
large
begin{align}
g(mu) &= eta \
&= theta \
&= mu \
end{align}
$$
1行目はリンク関数の定義、2行目は正準リンク関数の条件、3行目はセクション 1)で導出した正規分布のnatural parameterの等式です。よって、
$$
large
g(mu) = mu
$$
すなわち、正規分布の正準リンク関数は恒等関数であることが分かります。
ベルヌーイ分布
同様に、正準リンク関数はnatural parameterを(mu)の関数で表したものに対応するため、セクション 1)の結果から
$$
large
g(mu) = log(frac{mu}{1-mu})
$$
これは、ロジット関数と呼ばれているものです。ロジット関数の逆関数はシグモイド関数であり、すなわちシグモイド関数は線形予測子(eta)から期待値(mu)を求める関数であったことが分かります。
ポアソン分布
セクション 1)の結果から
$$
large
g(mu) = log(mu)
$$
正準リンク関数を用いない一般化線形モデルも存在します。ベルヌーイ分布では、ロジスティック回帰の他に、リンク関数に標準正規分布の累積密度関数、すなわち(g(mu)=boldsymbol{Phi}^{-1}(mu))を用いるプロビット回帰があります。
最尤法によるモデルパラメタの決定
では、尤度関数を導出し最尤推定によりモデルパラメタ(beta)および(y)の推定値(hat{mu})を決定します。
データを( { x^{(i)}, y^{(i)}}_{i=1}^N )とすると、各データサンプルに対する確率密度関数は次のようになります。
$$
large
begin{align}
p(y^{(i)}; theta^{(i)}, phi) = exp bigl( frac{y^{(i)}theta^{(i)}-b(theta^{(i)})}{a(phi)} + c(y^{(i)}, phi) bigr)
end{align}
$$
各サンプルごとに、dispersion parameterである(phi)は等しいが、natural parameterである(theta^{(i)})は異なることに注意してください。例えば、線形回帰においては、(y^{(i)}|x^{(i)} sim mathbb{N}(mu^{(i)}, sigma^2))です。
尤度関数は、次のようになります。ただし、尤度関数の引数を、各サンプルのnatural parameterを結合したベクトル(boldsymbol{theta})として表記しています。
$$
large
L(boldsymbol{theta}) = prod_{i=1}^N p(y^{(i)}; theta^{(i)}, phi)
$$
よって対数尤度関数は、
$$
large
l(boldsymbol{theta}) = sum_{i=1}^N frac{y^{(i)}theta^{(i)}-b(theta^{(i)})}{a(phi)} + sum_{i=1}^N c(y^{(i)}, phi).
$$
最尤法により、natural parameterは次のように決定します。
$$
large
begin{align}
hat{boldsymbol{theta}} &= argmax_theta l(boldsymbol{theta}) \
&= argmax_theta sum_{i=1}^N frac{y^{(i)}theta^{(i)}-b(theta^{(i)})}{a(phi)} + sum_{i=1}^N c(y^{(i)}, phi) \
&= argmax_theta sum_{i=1}^N y^{(i)}theta^{(i)}-b(theta^{(i)})
end{align}
$$
natural parameter (theta)とモデルパラメタ(beta)の関係はモデルによって異なりますが、ここではcanonical link function（正準リンク関数）を採用する場合を考えます。正準リンク関数を採用すると、natural parameter (theta^{(i)})と線形予測子 (eta^{(i)}=beta^mathrm{T} x^{(i)})は等しくなるため、次が成立します。
$$
large
begin{align}
hat{beta} &= argmax_beta l(beta) \
&= argmax_beta sum_{i=1}^N y^{(i)}beta^mathrm{T} x^{(i)}-b(beta^mathrm{T} x^{(i)})
end{align}
$$
あとは、セクション 1)で求めた(b(theta))の部分を確率分布に合わせて書き下せば、線形回帰やロジスティック回帰の尤度関数がこれまでと同じ形式になることを確認できます。
さらに、回帰による(y^{(i)})の予測値(hat{mu}^{(i)})は次のように求まります。
$$
large
begin{align}
hat{mu}^{(i)} = g^{-1}(hat{beta}^mathrm{T} x^{(i)})
end{align}
$$
例えば、ロジスティック回帰におけるリンク関数(g)であるロジット関数の逆関数(g^{-1})は、前回導入したシグモイド関数です。
まとめ
今回は、一般化線形モデルGLMの概要から、線形回帰やロジスティック回帰、ポアソン回帰など主要な回帰手法が、GLMに含まれることを見ていきました。
参考文献
[1] Generalized Linear Models – CMU statistics
[2] Dobson, Annette J., and Adrian G. Barnett. An introduction to generalized linear models. CRC press, 2018.
[3] Stanford cs229 lecture note 1