Np-Urのデータ分析教室

オーブンソースデータなどWeb上から入手できるデータを用いて、RとPython両方使って分析した結果を書いていきます

柔軟な確率分布を仮定して分析できる!一般化線形モデル(GLM)とは?

前回までは線形回帰の理論とそれを使った分析の実例について紹介しました。
【理論編】
randpy.hatenablog.com
【実践編】
randpy.hatenablog.com
randpy.hatenablog.com

しかし全てのデータを線形回帰で分析しようとすると、良い結果が得られないことがあります。
そこでよく用いられるのが一般化線形モデル(GLM : generalized linear model)という手法です。

今回は、線形回帰の問題点とそれを回避する一般化線形モデルについて紹介します。

線形回帰だとうまくいかない例

例えば、コンビニの1分間あたりの来店数と商品のバリエーション数の相関を分析したいとします。来店数が目的変数で商品バリエーション数が説明変数です。

つまり、置いてある商品の種類が多いコンビニほど、来店数も多いだろうという仮定を置いている、ということになります。

ここで、例えば以下のようなデータが得られたとします。

商品種類数 来店数
49 15
64 20
41 10
43 8
56 18
・・・ ・・・


横軸に商品種類数・縦族に来店数をとってグラフをプロットしてみるとします。

f:id:Np-Ur:20170625192638p:plain

何となく相関がありそうなのが分かるかと思います。
これを前回習った線形回帰に当てはめてみます。

最小二乗法に従って計算を行うと、最もフィットする直線を引くことができますね。

f:id:Np-Ur:20170625192832p:plain

しかしこの直線には不適切なところがあります。

商品種類数が0に近いところをよく見てみましょう。
来店数がマイナスという結果になっています。

「1分間に-4人来店する……?」という、
現実とはかなりかけ離れた推定となってしまいました。

上記のような直線ではなく、自由な線が引けるとしたら、下のようになるでしょう。
f:id:Np-Ur:20170625224527p:plain

「そんな線なんて推定できるの?」と思ったそこのあなた!!ご安心を!

本記事を最後まで目を通し、一般化線形モデルについての理解を深めれば、RやPythonなどを使って簡単に推定することができちゃいます!

それでは、一緒に勉強していきましょう!

一般化線形モデル

一般化線形モデルは、データのばらつきを正規分布以外(ポアソン分布や二項分布など)にも適用することができ、また目的変数から説明変数への影響の仕方も、この後説明する「線形予測子」と「リンク関数」を用いることで柔軟に変えることができます。

ここから先に進むために、とりあえずポアソン分布と二項分布については簡単に復習しておくと良いでしょう。

ポアソン分布について詳しくはポアソン分布:ある時間帯にかかってくる電話の数がわかる!??から。
ポアソン分布とは簡単に説明すると、ある事象が平均\(\lambda\)回発生すると分かっているときに、その事象が\(y\)回起こる確率を
$$\mathcal{f}(y \mid \lambda) = \frac{\lambda^y exp(-\lambda)}{y!}$$
と表したものです。

二項分布については二項分布とベルヌーイ分布 登校中にヤンキーに遭遇してしまう確率…?から復習しておきましょう!
二項分布も簡単に説明すると、ある事象が起こる確率が\(q\)と分かっているサンプルサイズ\(N\)の標本があったときに、そのうち\(y\)個にその事象が発生する確率を
$$\mathcal{f}(y \mid N, q) = {}_N C _y q^y (1-q)^{(N-y)}$$
と表したものです。

線形予測子とリンク関数

前述したコンビニの1分間あたりの来店数と、その商品種類数との因果関係について再度考えます。ここでは、来店数は非負で離散であり上限が無いため、前章で述べたポアソン分布に従うと仮定するのが良いでしょう。

目的変数\(y_i\)をあるコンビニ\(i\)で観測された1分間あたりの来店数、説明変数\(x_i\)をそのコンビニに置いてある商品の種類数としたとき(\(i=1,...,n\))、(何も考えずに)線形回帰を行うとすれば以下のような式になるでしょう。
\begin{eqnarray*}
y_i &=& a + b x_i + \epsilon_i\\
\end{eqnarray*}

しかし左辺(来店数)は実際には非負の値を取りますが、線形回帰の場合\(x\)の値によっては右辺が負になる場合があり不適切です(先ほど図で示したように)。

そこで、通常ポアソン分布に従うデータを回帰分析したい場合は、
$$\lambda_i = exp(a + b x_i)$$
のように変換する処理を行うことが多いです。観測された\(y_i\)というデータは平均\(\lambda_i\)に従うポアソン分布から生成されたと考え、その\(\lambda_i\)は\(exp(a + b x_i)\)によって決まると考えたモデルとなります。
$$p(y_i \mid \lambda_i) = \frac{\lambda_i^{y_i} exp(-\lambda_i)}{(y_i)!}$$

この処理のメリットとして、\(x\)がどんな値を取ろうとも、\(\lambda_i = exp(a + b x_i)\geq 0\)となっている点があります。

よって\((a + b x_i)\)は\(-\infty \sim+\infty\)の範囲を動きますが、expの中に入れることで\(0 \leq \lambda_i \leq \infty\)である\(\lambda_i\)を推定できるようになります。

また両辺に対数を取って式変更を行うと、
$$\log \lambda_i = a + b x_i$$
と書けます。

この右辺の線形結合された式を「線形予測子」、線形予測子を導くような左辺の関数を「リンク関数」と呼びます!

中でも今回は対数をとっているので、「対数リンク関数」といいます。

ポアソン分布を仮定したデータを使った一般化線形モデルの場合、リンク関数に対数リンク関数を用いることが一般的です。
ちなみに一般線形モデル(線形回帰)は、データのばらつきに正規分布・リンク関数をそのまま(変形なし)と考えると、一般化線形モデルの1つとみなすことができます。

なので線形回帰と一般化線形モデルは別物と捉えずに、「線形回帰は一般化線形モデルの特殊系である」と考えるのが良いでしょう!

最尤法

前回線形回帰について説明した際は、\(a\)や\(b\)のパラメータを推定するときは最小二乗法(OLS)を使うと書きました。

この最小二乗法は線形回帰でしか適用できず、代わりに一般化線形モデルでは「最尤法(さいゆうほう)」という手法を使ってパラメータを推定します。
なお、最尤法は線形回帰でも使うことができます。

モデル 一般線形モデル(線形回帰) ポアソン回帰
確率分布 正規分布 ポアソン分布
リンク関数 そのまま 対数
解法 最尤法 or 最小二乗法 最尤法

最尤法については別の機会に改めてしっかり記事を書きたいと思いますが、ざっくり言うと、どんなパラメータを与えたときに、最もその手持ちのデータを再現できそうかを求める手法です。

最尤法を直観的に理解する例として、コインを2回投げてみて表が1回、裏が1回出たときを考えます。
このとき表が出る確率(パラメータ)を最尤法を使って推定してみます。

例えば表が出る確率が\(1/4\)だと仮定してみましょう。このとき表が1回、裏が1回出る確率は$$1/4 * 3/4 = 3/16$$となりますね。
この\(3/16\)という値を尤度と言います。尤度はパラメータの値を定めたときに、そのデータが起こりうる(表が1回、裏が1回出る)確率を意味しています。

また、表が出る確率を\(1/5\)としたとき、尤度は$$1/5 * 4/5 = 4/25 < 3/16$$となります。表が出る確率が\(1/5\)のときよりも\(1/4\)のときの方が尤度が大きくなりました。

「表が出る確率を\(1/5\)とするよりも\(1/4\)とした方が尤度が大きい」ということは、「表が出る確率を\(1/5\)とするよりも\(1/4\)とした方が、手持ちのデータ(表が1回、裏が1回出たという事象)を再現しやすい」……ということは、「表が出る確率は\(1/5\)よりも\(1/4\)とした方が良さそうじゃない!?」
と考えられますね。

このようにして、度が大きいパラメータを推定する方法が最尤法です!

ちなみにこの例では表が出る確率を\(1/2\)としたとき尤度は\(1/4\)で最大となります。
当然といえば当然ですね…。

ポアソン回帰のイメージ

ポアソン回帰(ポアソン分布における一般化線形モデル)を最尤法を使って実際に解く過程を確認します。

先ほどのコンビニの来店数\(y\)と商品の種類数\(x\)について、
$$\lambda_i = exp(a + b x_i)\\
p(y_i \mid \lambda_i) = \frac{\lambda_i^{y_i} exp(-\lambda_i)}{(y_i)!}$$
これらの式を思い出してください。

求めたいのはパラメータ\(a,b\)です。(たまに「あれ?何求めるんだっけ?」と陥ることがあるので確認しておきましょー)

ここで、パラメータ\(a,b\)がある値に定まったと考えます。このとき例えばコンビニA(商品の種類数\(x_A\)と来店数\(y_A\))があるとすると、
$$\lambda_A = exp(a + b x_A)$$
とまず\(\lambda_A\)が求まります(代入するだけ!)。

この\(\lambda_A\)によって来店数\(y_A\)のポアソン分布の形が決定されるので、このとき\(y_A\)という事象が起こりうる確率は


\begin{eqnarray*}
p(y_A \mid \lambda_A) &=& \frac{\lambda_A^{y_A} exp(-\lambda_A)}{(y_A)!}\\
&=& \frac{(exp(a + b x_A))^{y_A} exp(-(exp(a + b x_A))}{y_A!}
\end{eqnarray*}
となります(これまた代入するだけ!)。

この数式自体が難しいと思う方もいるかもしれませんが、「対数リンク関数を使ってうまいこと線形予測子を変換して\(\lambda_A\)を求めて、それを使ってポアソン分布に従う\(y_A\)の発生確率を考えている!」という認識で一旦は大丈夫です!

これにて、パラメータ\(a,b\)がある値に定まったときのコンビニA(商品の種類数\(x_A\)と来店数\(y_A\))というデータが発生し得る確率が求まりました。
じゃあコンビニB(商品の種類数\(x_B\)と来店数\(y_B\))は?コンビニCは?…と考えていけば全て同様に求めることができます。

これらを全てのコンビニで掛け合わせることで、以下のようになります。


\begin{eqnarray*}
L(a, b \mid Y) &=& \prod_{i} \frac{(\lambda_i)^{y_i} exp(-\lambda_i)}{y_i!}\\
&=& \prod{i} \frac{(exp(a + b x_i))^{y_i} exp(-(exp(a + b x_i))}{y_i!}
\end{eqnarray*}
この計算はつまり、あるパラメータ\(a,b\)が与えられたときに手持ちのデータ(各コンビニの来店数と商品種類数)が再現できる確率となります。

ちなみに\(\prod\)はあまり馴染みない方もいるかもしれませんが、すべて掛け合わせるという意味の記号です。(詳しくは総乗をご覧ください。)


なお、これを尤度関数と呼びます。
この尤度関数を最大にする(手持ちのデータを最も再現できる)ようなパラメータ\(a\)と\(b\)を、最尤法によって求める、というわけです。

実際にRなどを使って分析する場合は、確率分布とリンク関数を指定しあげれば、コード1行で推定することができます。
なので上記のことは忘れてしまっても使えることは使えるのですが、頭の片隅で何をやっているかイメージしてもらえればと思います!

一般化線形モデルの具体例 - ロジスティック回帰

目的変数\(y_i\)が\(0\)または\(1\)をとるようなとき、適切な説明変数を用意することで\(y_i\)が\(0\)なのか\(1\)なのかの確率を推定する問題を考えてみます。例えば、

  • 通意表の成績やセンター試験での点数から大学入試の合否(合格なら\(y_i=1\)、不合格なら\(y_i=0\)として)の確率を推定したい
  • 健康診断の数値からある病気にかかっているかどうか(病気なら\(y_i=1\)、健康なら\(y_i=0\)として)の確率を推定したい

ということがあるかと思います。これは二項分布とベルヌーイの分布 登校中にヤンキーに遭遇してしまう確率…?で紹介したベルヌーイ分布ですね!

今回簡単のため、大学入試の合否を目的変数\(y_i\)に、センター試験の合計点数を説明変数\(x_i\)としてこのような問題を考えます(実際には説明変数にはもっと多くの変数を用意することが通常です)。

このとき、線形回帰を適用して
\begin{eqnarray*}
y_i &=& a + b x_i + \epsilon_i\\
\end{eqnarray*}
などするのは当然おかしい結果となります。なぜなら、右辺は\(-\infty \sim+\infty\)の範囲を動けるのに対し、左辺は\(0\)または\(1\)しか動けないからです。

線形予測子とリンク関数

そこで、前章で行った時と同様に、線形予測子とリンク関数を用意します!目的変数が2通り(1 or 0)の場合では、次のように変換します。
\begin{eqnarray*}
\log \frac{q_i}{1-q_i} = a + b x_i
\end{eqnarray*}

\(q_i\)とは大学入試の合格する確率、つまり\(y_i=1\)となる確率を意味しています。よって不合格となる確率は\(1 - q_i\)です。

式展開をすると、
$$q_i = \frac{1}{1 + exp(-(a + b x_i))}$$

このようにも表せますね。
前者の関数をロジット関数、後者の関数をロジスティック関数と呼び、このロジスティック関数を用いてパラメータを推定する方法をロジスティック回帰といいます。

ロジスティック関数を使うメリットとして、\(x\)がどんな値を取ろうとも、\(0 \leq q_i \leq 1\)となっている点です!
このように線形予測子がどんな動きを取ろうとも、このリンク関数と組み合わせることで、確率を表現できるというわけです。
これはすごい発明ですね!(皆さんテンション上がっていますか?)

なお、推定にはポアソン回帰同様最尤法が通常用いられます。

パラメータ\(a,b\)がある値に定まったと考えます。このとき例えばセンター試験で\(x_A\)と取ったA君がいたとすると、
$$q_A = \frac{1}{1 + exp(-(a + b x_A))}$$
と大学合格する確率\(q_A\)が求まります。

このときもしA君が合格していたら、尤度はそのまま\(q_A = \frac{1}{1 + exp(-(a + b x_A))}\)となり、不合格していた場合の尤度は\(1-q_A = 1- \frac{1}{1 + exp(-(a + b x_A))}\)となります。

このような尤度を全ての受験生について計算して、すべてかけると尤度関数が求められます。(ポアソン回帰のときと一緒です!)
そしてこの尤度関数を最大にするようなパラメータを求めれば終わりです。

まとめ

今回は前回前々回と学習した正規分布以外に確率分布を用いて、分析を行う「一般化線形モデル」について学習しました。

正直何も考えず正規分布だけ仮定して線形回帰している論文があまりに多いので、これを意識して柔軟に分析できるだけで、かなりレベルアップできると思います!

次回は実践編ということで、あるデータを一般化線形モデルを使って分析してみます!
こうご期待です!!


また一般化線形モデルについて、下の文献がかなり分かりやすく解説しているので、気になる方は読んでみてください!

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)