この世は直線関係や正規分布に従う変数ばかりではない

　今まで私たちはいろんな統計モデルを議論してきたかのように思えるが、その実態はほぼすべてが変数の関係が直線関係（線形）であり、誤差が正規分布に従うことを仮定した解析であった（分散分析、線形回帰、重回帰、非線形最小二乗法は非線形であったが）。けれども、この世のデータはそのような関係ばかりではなく、例えばデータが離散的だったり、左右非対称な分布だったり、直線関係とは呼べないような変数間の関係があることを、前回紹介した。次のモチベーションとしては、これらのデータも無理なく、解析したくなってくるだろう。そして、可能であれば、せっかく培ってきた線形回帰の考え方をベースに理解していきたいものである。そこで、編み出されたのが、今回紹介する一般化線形モデルGenelarized Linear ModelあるいはGLMと呼ばれる手法である。この方法は、前述の望みの多くをかなえてくれる。名前の通り線形回帰モデルを拡張した解析手法でありながら、”ある工夫”をすることで正規分布以外のデータや直線関係でないデータを解析できるのである。本稿では、GLMの考え方の概要について紹介する。

一般化線形モデルはリンク関数がミソ

　今一度、線形回帰（線形モデル）がどのようなモデル式になっていたかをおさらいしよう。例えば、1変数の場合なら以下のように表記できる。

もし、もっと多くの変数があるのであれば下記のように表記できる。

あるいはベクトルを使った表記を使い、ベクトルの内積を認めれば、以下のようにシンプルに表記も可能である。

このような説明変数を線形結合（定数倍して足す）して表した式を、線形予測子linear predictorと呼ぶ。さて、もともと線形回帰とは、説明変数が線形結合で表されているゆえに、基本的に直線の変化の関係しか記述できないのであった。だから、例えば下記のような、被説明変数が指数関数のような挙動を示す場合、今までの線形回帰の解析の範疇から飛び出してしまう。

けれども、もし、説明変数の線形結合で被説明変数を無理なく表現できるのであれば、解釈が容易となるため利便性が高いだろう。そこで考えてみると、この式においては指数関数の指数部分には線形予測子がそのまま入っている。だとすれば、両辺の対数をとり、log(y) = Yとすれば、以下のような説明変数の線形結合に帰着できるだろう。

新たな説明変数Y＝線形予測子という形になるので、基本的には今までと同様にしてパラメータβを推定することが可能である。このように、一般の回帰曲線の式を線形の式に変形するような関数、今回であれば対数関数logをリンク関数link functionと呼ぶ。

このリンク関数による変形がGLMの工夫の一つである。逆に、リンク関数の逆関数が元の回帰曲線の式の形を表している。後述するように、元の回帰曲線の性質に合わせた確率分布を指定することで、より妥当性の高い解析を実現できる。

リンク関数の逆関数に合わせた確率分布が指定できる

　また、GLMではリンク関数の逆関数に合わせたデータの確率分布を指定することができ、より妥当な解析を実現可能である。

例えば、上記の対数関数がリンク関数となっているとき、その逆関数は指数関数である。もとの回帰式の形が指数関数であったことと一致する。指数関数の性質として、0以下の値をとらない。つまり、被説明変数yは0より小さな値をとらないような実測値であるはずだし、そのyを生み出す背後のメカニズム（確率分布）も同様に0より小さな値をとらないようになっていてほしい。しかし、今までの線形回帰ではこれは不可能である。なぜならば、線形回帰の仮定では平均値に対して被説明変数のデータのばらつき（残差）は正規分布に従って生成されると考えていた。正規分布は-∞から+∞までの値をとる可能性のある確率分布である。したがって、どのような平均値であっても、正規分布を仮定する以上、0より小さな値が出る確率を0にすることは不可能である。

　このような指数関数と相性のいい確率分布の一つが、ポアソン分布Poisson distributionである。この確率分布は0以上の離散値をとる確率分布である。あるいは、ガンマ分布Gamma distributionも候補として上がるだろう。こちらは0より大きな連続値をとる確率分布である。これらの詳細は実際の解析時に紹介するが、データのばらつきがポワソン分布やガンマ分布に従うと考えられれば、無理なく指数関数の回帰曲線を描くことが可能となるだろう。

　逆に取得したデータが特定の値域を持つことが想定される場合は、今までのようにただ線形回帰すればよいとは限らないことを意味している。例えば、種子の数をデータとして取得した場合、その数は0より小さな値になることはないだろう。あるいは重さをデータとして取得した場合、重さは0以下の値をとることはないだろう。それゆえ、確率分布に正規分布を仮定することが適切ではない場面が生じることがある（数や重さであっても正規分布を仮定して解析していい場面もあるが）。このようなデータの解析をするときは、おのずとGLMの必要性に駆られるわけである。

最尤法：手元のデータを最も高い確率で再現できるモデルを探索する

　線形回帰モデルでは最小二乗法という、「予測からのデータのずれ＝残差（平方和）が最小になるようにパラメータを設定する」という方法をとっていた。しかしこの方法は、残差が正規分布に従うことを仮定した方法であり、確率分布を正規分布以外に仮定した場合は必ずしも適していない。また、モデル式はリンク関数の存在や正規分布以外の分布をとるなどで複雑になっており、残差平方和を最小にするというのは容易でない。そこで、モデル式が目指すところから見方を変えてみよう。

　モデル式は手元のデータをうまく予測してほしいわけなのだから、目指すべきモデルは手元のデータを最も高い確率で再現できるものだろう。このようなモデルの探索の仕方を最尤法（さいゆうほう）most likelihood methodと呼ぶ。その見方の別の視点が、最小二乗法のような予測からの誤差が最小になっているモデル、でもあったわけである。さて、手元のデータを最も高い確率で再現できる、とは、数式でどのように表現すればよいだろうか。

　正規分布を例にとって話をしてみよう。説明変数は0個で（線形回帰ページの追記2の状況）、ある正規分布に従うと考えられるサンプルサイズnのデータが存在する。例えば、下図のようにデータ（白丸）が一直線上に並んでいる状態を考えよう。

ここで、このデータから母集団となる正規分布の平均値と分散を最尤推定することを考える。上記のデータに対して、なにか適当な正規分布に対応させてみる。

すると、上図において、ある正規分布に対して各データ点が生じる確率（正確には正規分布の場合は確率密度）を青線のように図示することができる。ここで、青線＝データ点が生じる確率のすべての積を尤度likelihoodと呼ぶ。尤度はまさに手元のデータが生じる確率、というべきものである。下図のように、平均と分散を自由に変更して、つまり尤度をパラメータの関数と考えて、この尤度が最大になった平均と分散を推定値とする方法が、最尤法なのである。このとき、尤度のパラメータに関する関数を尤度関数likelihood functionと呼ぶ。

パラメータの関数として考えるとか、関数の最大・最小を求めるとか、今までもさんざんやってきたから、もはやよくなじんでいるだろう。最小二乗法と最尤法は、モデルの当てはまりの指標が残差平方和か、尤度か、最小値を求めるのか、最大値を求めるのかが違うだけで、その手続きは結構似ている。また、結果も類似していて、線形回帰において最小二乗法で求めた平均や分散、そのほか係数などのパラメータは、誤差を正規分布、リンク関数を恒等関数としたときの最尤法で求めたものと一致する（補足）ことが知られている。ただし、上記の仮定を満たさない条件であっても最尤法は利用可能である点で、最小二乗法の上位互換ともいえる方法である。今後、最尤法や尤度というコンセプトは今後重要になるので、そこだけは押さえて、今後の解説に挑んでほしいと思う。

補足：線形回帰の最小二乗法は誤差を正規分布、リンク関数を恒等関数としたときの最尤法と同じ結果を返す

　最尤法によるパラメータ推定を理解するためには、最小二乗法の時以上に、データがある確率分布から生成される過程を具体的に想像できる必要が出てくる。最尤法の具体的手続きを理解するついでに、「線形回帰の最小二乗法は誤差を正規分布、リンク関数を恒等関数としたときの最尤法と同じ結果を返す」ことを示そう。

　今まで正規分布の形だけを示しており、ほとんど式の形については述べてこなかった。ここで、ある確率変数xが平均μ、標準偏差σである正規分布が関数としてどのように表されるかを示す。

式を微分すると、x < μではf'(x) > 0、x > μではf'(x) < 0となることから、この正規分布の関数はx = μで頂点y = 1/√(2πσ^2)となる一山形の曲線になることがわかる。また、x → ±∞でf(x) → +0なので、xがμから離れると0に正のほうから漸近することもわかる。このように連続的な確率変数に対して、その確率密度を表した関数f(x)を確率密度関数（Probability density function, PDF）と呼ぶ（離散的な確率変数の場合は、確率質量関数（Probability mass function, PMF））。確率変数があるx0という値だとしたら、その点での確率密度はf(x0)である、ということである。

　さて、通常、最小二乗法をおこなってパラメータ推定するような例を使って、最尤推定を行ってみる。説明変数を一つ持ち（変数xとする）、被説明変数yに関して、パラメータβ0、β1を推定する。以下のような仮定を置く。

1. y ~ β0 + β1×xという直線の関係がある

2. 誤差 = 実測 - 予測（β0 + β1×x）は正規分布する

本項目の例ではリンク関数に対数関数を紹介したが、リンク関数を恒等関数、つまり端的に言えば変換しない場合でも一般化線形モデルの手法は使える。1の主張は、このような無変換の形で、説明変数と被説明変数が結ばれているということである。2の主張は線形回帰の時同様である。最尤推定における線形回帰の仮定は以上である。

　次に、尤度を計算する。尤度は実測データの確率密度の積だった。今、説明変数x = {x1, x2, ……, xn}に対応する被説明変数がy = {y1, y2, ……, yn}とする。簡単な状況は以下のような図となる。また、その時のあるデータペア（xi, yi）において、xiが得られた時にyiが得られる確率密度（条件付き確率p(yi | xi)とあらわす）を示す。

この時の尤度はつまり、条件付き確率p(yi | xi)をi = 1からnまで掛け合わしたものである。したがって下記のような式になる。

積の式は取り扱いが面倒なので、両辺の対数をとる。すると、和の式に帰着可能である。このような形にした尤度を対数尤度log likelihoodと呼ぶ。

今回、あたかも式の見通しを良くするために対数変換したが、実際のコンピュータの計算でも対数をとったうえで計算を行う。というのも、確率密度は1より小さな値であるから、それらの何十もの積である尤度は極めて桁数の小さな値となる。コンピュータでは、計算できる桁数が決まっている（大体、10の-320乗程度が限界のようだ）ので、このような問題を避けるために実用上でも対数をとっているのである。

　さて、ここでの目的は、尤度を最大にするβ0、β1（とσ）を探すことである。対数尤度は以下のような多変数関数とみなして、尤度を最大にするβ0、β1、σを探そう。

ところでよく考えてみると、σを定数としたとき、対数尤度の前半はβ0、β1によらない定数だから、後半の式を最小にすれば対数尤度は最大値となることがわかる。特に、後半の式の分子を最小にすればよい。これは最小二乗法と同じ問題を解いていることになる。ゆえに解くまでもなく、過去の記録を参照すると答えは以下のとおりである。なお、ある関数f(x)を最大化するようなパラメータxをarg max f(x)、最小化するようなパラメータxをarg min f(x)と表記する。argとは関数の引数argumentのことである。

確かに、最小二乗法で求めた係数と一致する。また、正規分布の分散は以下のように求められる。

つまり、標本分散は残差平方和をサンプルサイズで割ったものとなる。これも線形回帰において暗黙の了解として示した不偏分散と類似した形となっている。そう、実は最尤推定で求めた分散は不偏分散ではない。不偏分散を求めるときは、サンプルサイズnを、自由度分差し引いた値で除す必要がある。つまり、説明変数が一つあるならn-2で除し、説明変数が一つもないならn-1で除す。

　以上のように、確かに誤差を正規分布、リンク関数を恒等関数としたときの最尤法は、線形回帰の最小二乗法と同じ結果を返すのである。