過分散などの対処すべき課題

　以前紹介したように一般化線形モデルでもまだ対処できない課題はある。例えば、カウントデータの過分散はポアソン分布や二項分布だけでは対処できない。この問題以外にもいやらしい問題として、下図のように全体傾向が群内傾向と異なる場合もある。

全体傾向が群内傾向と異なる場合というのは、解析結果次第では全く異なる結論を下すことになる。このようになる良い思考実験(架空データの話なので真実とは限らないが、それっぽい話を用意した)は「体重」と「足の速さ」の関係だ。小学生の体重と足の速さを測ってみると、「体重が増えるほど足が速い」という結果が得られる（だろう）。

でも、なんとなく、これは直感に反している気がする。この問題、実は年齢（学年）の効果を考えていない。つまり、学年が上がるほど、体重も増え、足も速くなるために、「体重が増えるほど足が速い」結果が得られるのだ。同一学年内で比較すると、「体重が増えるほど足が遅くなる」傾向が検出されるようになる（と直感的には思える）。つまりは、一種の擬相関である。ただし、実際のデータ解析では、年齢のようなわかりやすい、説明変数として組み込めるデータとして得られないことも多い。何が具体的に異なるかはわからない、A群、B群……というグループとしてデータが存在することが多いだろう。

　この「過分散」と「全体傾向≠群内傾向」の問題は、一見すれば全く別問題に見える。けれども、それぞれのメカニズムを考えたときに、これらの問題が生じたのは、全体の中にある階層構造の効果を見逃したために起こっている、共通の問題だと考えることができる。また、上記のような問題が一見生じてなくても、階層構造が存在するなら、それを解析に組み込むのは当然の判断である。今後の課題は、こういった階層構造の効果をどうやって取り込んでいくのか、という問題に帰着するわけである。

一般化線形混合モデルGeneralized linear mixed model

　そこで、一般化線形モデルをさらに拡張することを考える。一般化線形モデルは、fをリンク関数とすれば、以下のような線形予測子の形に帰着して、各説明変数の係数を最尤法によって推定する。

ここにこの線形予測子にばらつきを生じさせる項rを導入する。

モデルの式の形としては、これが一般化線形混合モデルの一つの形である。この時、自分が明示的に計測できた説明変数、あるいは効果に興味がある変数を固定効果fixed effect、どんな違いがあるかわからない変数、あるいは効果に興味のない変数をランダム（変量）効果random effectと呼ぶ。ここで、ランダム効果はある確率分布から生成された値である、と考える。一般化線形混合モデルGeneralized linear mixed model (GLMM)とは、これらの固定効果とランダム効果の組み合わせによってモデルを表現することから、そう呼ばれる。また、データに階層構造を仮定していることから、マルチレベルモデルMulti-level modelや階層線形モデルHierarchial linear modelなどの1種として扱われる。

　だがこれだけだとまだ意味が分からないと思う。そこでもっと具体的な例を考えよう。ここでは、まず、ひたすらに分布の混合によって何が起こるかを見ていくことにして、次の項でなぜこんなことをするのかを考える。

・ランダム効果が正規分布に従う対数リンク-ポアソン回帰

　例えば、ランダム効果は平均0、分散σ^2の正規分布に従うとする。そして被説明変数yはポアソン分布から生成されると考える。そうすると、一般化線形混合モデルは以下のように書き直すことができる。ややこしくて恐縮だが、ここでの~（チルダ）は、「y ~ 分布(線形予測子)」とあった時に、「yは線形予測子をパラメータとして持つ分布から生成される」、と読んでほしい。今まで使ってきたチルダの用法とは異なっている。

f^(-1)はリンク関数fの逆関数である。対数リンクを使ったなら指数関数になっているわけである。ランダム効果の文字を見かけ上、消すこともできて、以下のようにできる。

さらに、対数リンクの時、λの対数log(λ)が正規分布に従っている。このように乱数の対数が正規分布に従っているとき、もとの乱数λは対数正規分布log-normal distributionに従う。ゆえに以下のように書き換えることができる。

対数正規分布とは、以下のように正規分布が与えられた時に、その横軸を指数変換した軸の分布である。

ここまでで、一般化線形混合モデルの実態を理解するためのピースがそろった。今までの一般化線形モデルでは、yが平均λのポアソン分布に従ってばらつき、λは線形予測子に一致してばらつかない、と考えていた。一方、一般化線形混合モデルは、λも線形予測子によって決まった値を示すのではなく、対数正規分布に従ってばらつく。つまり、2段階でばらつくわけである。

　このとき、どんな分布になるかというと、まず、対数正規分布に従って、ある確率p1でλ1が生成され、さらにそのλ1に対してある確率分布のポアソン分布poisson(λ1)が決まる。今度はある確率p2でλ2が生成され……となり、最終的にこれらのポアソン分布の和が求める分布となる。つまり、p1×poisson(λ1)+p2×poisson(λ2)……という感じである。λは連続的に変化するので、実際には、λに関する対数正規分布とポアソン分布の積の積分となる。

当然、そうなったときに出来上がる分布は、下記のようにポアソン分布よりもばらつくだろう。

試しに、自分で上記の設定に合う確率分布を生成してみよう。例えば、下記のように関数を組むことができる。

------------------------------------------------------

library(plotn)

n <- 100000

rlognorm_pois <- function(n, mean, sd){

  sapply(exp(rnorm(n = n, mean = mean, sd = sd)), 

         function(tmp) rpois(1, lambda = tmp))

}

------------------------------------------------------

上記の関数は、まず、正規分布から乱数を生成する。その生成された乱数を指数変換する。指数変換された乱数は、対数をとれば正規分布に従っていたのだから、指数変換された乱数は対数正規分布に従う。この指数変換された乱数を平均とするポアソン分布を生成するのが、上記の関数である。では、いくつかのパラメータを入力して、挙動を見てみよう。

------------------------------------------------------

m <- 0

s <- 0.5

d_mix <- rlognorm_pois(n, m, s)

histn(d_mix, freq = F) #図1の描画

m <- -1

s <- 0.5

d_mix <- rlognorm_pois(n, m, s)

histn(d_mix, freq = F) #図2の描画

m <- 1

s <- 0.5

d_mix <- rlognorm_pois(n, m, s)

histn(d_mix, freq = F) #図3の描画

m <- 1

s <- 1

d_mix <- rlognorm_pois(n, m, s)

histn(d_mix, freq = F) #図4の描画

------------------------------------------------------

確かに生成された乱数の分布はポアソン分布よりも大きくばらつき、狙った通りの効果が得られている。

・ランダム効果がガンマ分布に従う恒等リンク-ポアソン回帰

　ポアソン回帰において重要なのは、λが非負の値しかとらないパラメータであることだ。先ほどまでは、線形予測子が正規分布に従うことと、対数リンクを使うことで、本来-∞から+∞の値をとるはずのものを、非負の値に変換した。この非負の値が従うのが、対数正規分布でありそこからλが生成されると考えていたわけだ。では、はじめから非負の値を生成する確率分布からλが生成されると考えてはどうだろう。例えば、ガンマ分布は大変丁度良い性質がある。形状パラメータkを持つガンマ分布からλが生成されたと考えれば、ガンマ分布は非負の確率分布だからλは特に変換をしなくても（つまり、恒等リンクに相当）、そのままポアソン分布のパラメータとして使える。

　最終的にyが従う分布は、λに関するガンマ分布とポアソン分布の積の積分となる。実は、この分布は負の二項分布negative binomial distributionになることが知られている。ガンマ分布の形状パラメータk、尺度パラメータθとすると、負の二項分布におけるサイズパラメータはk、確率パラメータは1/(1+θ)となる。負の二項分布はポアソン分布よりもばらつきが大きい分布で、ポアソン分布に対して過分散なカウントデータに使える分布である。

　この場合のGLMMは負の二項分布に従うGLMと解釈してもいい場合がある。なお、線形予測子をガンマ分布に取り込む時、ガンマ分布は正の値しかとれないので、線形予測子も-∞から+∞の値ではなく、0より大きな値をとるように変換が必要になる。リンク関数の考え方と同じように、exp(線形予測子)がガンマ分布に従う、と考えればよいだろう。

　では、この確率分布に従う確率分布を生成してみよう。

------------------------------------------------------

rgamma_pois <- function(n, shape, scale){

  sapply(rgamma(n, shape = shape, scale = scale), 

         function(tmp) rpois(1, lambda = tmp))

}

shape <- 1

scale <- 1

d_mix <- rgamma_pois(n, shape, scale)

histn(d_mix, freq = F) #図5の描画

shape <- 3

scale <- 1

d_mix <- rgamma_pois(n, shape, scale)

histn(d_mix, freq = F) #図6の描画

------------------------------------------------------

確かに生成された乱数の分布はポアソン分布よりも大きくばらつき、狙った通りの効果が得られている。

・ランダム効果が正規分布に従うロジットリンク-二項回帰

　続いては、二項分布との混合分布を考えよう。ランダム効果は平均0、分散σ^2の正規分布に従うとする。そして被説明変数yは二項分布から生成されると考える。そうすると、一般化線形混合モデルは以下のように書くことができる。

ランダム効果の文字を見かけ上、消すこともできて、以下のようにできる。

さらに、ロジットリンクの時、pのロジット変換logit(p)が正規分布に従っている。このように乱数のロジット変換が正規分布に従っているとき、もとの乱数pはロジット-正規分布logit-normal distributionに従う。ゆえに以下のように書き換えることができる。

ロジット-正規分布とは、以下のように正規分布が与えられた時に、その横軸をロジスティック変換した軸の分布である。

pは連続的に変化するので、pに関するロジット-正規分布と二項分布の積の積分が実際にデータが従う確率分布ということになる。

　では、この確率分布に従う確率分布を生成してみよう。

------------------------------------------------------

rlogitnorm_binom <- function(n, mean, sd, size){

  sapply(1/(1+exp(-rnorm(n = n, mean = mean, sd = sd))), 

         function(tmp) rbinom(1, size = size, prob = tmp))

}

size <- 40

m <- 0

s <- 1

d_mix <- rlogitnorm_binom(n, m, s, size)

histn(d_mix, freq = F) #図7の描画

m <- -1

s <- 1

d_mix <- rlogitnorm_binom(n, m, s, size)

histn(d_mix, freq = F) #図8の描画

m <- 1

s <- 1

d_mix <- rlogitnorm_binom(n, m, s, size)

histn(d_mix, freq = F) #図9の描画

m <- 0

s <- 3

d_mix <- rlogitnorm_binom(n, m, s, size)

histn(d_mix, freq = F) #図10の描画

------------------------------------------------------

確かに生成された乱数の分布は二項分布よりも大きくばらつき、狙った通りの効果が得られている。さらにとても興味深いのが、図10でランダム効果の従う正規分布の標準偏差が大きくなると、混合分布からは中間的な値よりも0や最大値の40付近の極端な値が生じやすくなる。これは二項分布にはなかった混合分布における新しい性質で、このような確率分布が記述できれば、極端な値ばかりが存在する状況、例えば、結実率が0% or 100%の場合が多い、みたいな状況でも解析を行うことができる。

・ランダム効果がベータ分布に従う恒等リンク-二項回帰

　二項回帰において重要なのは、確率pが0~1の値しかとらないパラメータであることだ。先ほどまでは、線形予測子が正規分布に従うことと、ロジットリンクを使うことで、本来-∞から+∞の値をとるはずのものを、0~1の値に変換した。では、はじめから0~1の値を生成する確率分布からpが生成されると考えてはどうだろう。例えば、ベータ分布は大変丁度良い性質がある。dispersionパラメータφを持つベータ分布からpが生成されたと考えれば、ベータ分布は0~1で定義される確率分布だからpは特に変換をしなくても（つまり、恒等リンクに相当）、そのまま二項分布のパラメータとして使える。

　最終的にyが従う分布は、pに関するベータ分布と二項分布の積の積分となる。実は、この分布はベータ二項分布beta binomial distributionになることが知られている。ベータ二項分布は二項分布よりもばらつきが大きい分布で、二項分布に対して過分散なカウントデータに使える分布である。

　この場合のGLMMはベータ二項分布に従うGLMと解釈してもいい場合がある。なお、線形予測子をベータ分布に取り込む時、ベータ分布は0~1の値しかとれないので、線形予測子も-∞から+∞の値ではなく、0~1の値をとるように変換が必要になる。リンク関数の考え方と同じように、logistic(線形予測子)がベータ分布に従う、と考えればよいだろう。

　では、この確率分布に従う確率分布を生成してみよう。

------------------------------------------------------

rbeta_binom <- function(n, mu, phi, size){

  sapply(rbeta(n, shape1 = mu*phi, shape2 = (1-mu)*phi), 

         function(tmp) rbinom(1, size = size, prob = tmp))

}

size <- 40

mu <- 0.5

phi <- 5

d_mix <- rbeta_binom(n, mu, phi, size)

histn(d_mix, freq = F) #図11の描画

mu <- 0.7

phi <- 5

d_mix <- rbeta_binom(n, mu, phi, size)

histn(d_mix, freq = F) #図12の描画

mu <- 0.5

phi <- 10

d_mix <- rbeta_binom(n, mu, phi, size)

histn(d_mix, freq = F) #図13の描画

mu <- 0.5

phi <- 1

d_mix <- rbeta_binom(n, mu, phi, size)

histn(d_mix, freq = F) #図14の描画

------------------------------------------------------

確かに生成された乱数の分布は二項分布よりも大きくばらつき、狙った通りの効果が得られている。さらにとても興味深いのが、図14でランダム効果の従うベータ分布のdispersion parameterが小さいと、混合分布からは中間的な値よりも0や最大値の40付近の極端な値が生じやすくなる。これは二項分布にはなかった混合分布における新しい性質で、このような確率分布が記述できれば、極端な値ばかりが存在する状況、例えば、結実率が0% or 100%の場合が多い、みたいな状況でも解析を行うことができる。

　最終的に出力される分布がわかるのであれば、尤度の関数も決められる。あとはGLMと同じように混合分布を支配するパラメータを最尤推定すれば、求めたい値が求まるという寸法である。いずれの分布の混合を見ていても気づくと思うが、分布の混合によって、より自由度高い確率分布の記述が可能となっている。特に、ポアソン分布や二項分布では記述しきれない過分散は、GLMMによって容易に克服が可能となる。

　ただし、やはりというべきか、分布の混合で必要になる積分計算はやや重たい。後述するような、大量の説明変数があったり、さらにそれらに対して「傾き」にもランダム効果を入れる場合のように複雑なモデルを組むほど、尤度の計算は困難を極め、極端に重たくなっていく。

なぜ、固定効果ではなくランダム効果にするのか？

　なぜ、ランダム効果を固定効果として考えずに解析をしているのだろう。それは自由度の節約、という考えにある。今一度、GLMMのモデル式に立ち返って考える。

例えば、ここでランダム効果rがyごとに自由に決まるパラメータだったとする。ランダム効果の意味合いとして、自由に値が決まってもいいだろう。しかし、そのモデルはすぐに破綻することがわかる。もし、データ点がn個、係数を決めるパラメータの数がm個だとしたら、このモデルで推定するべきパラメータはn + m個になる。データ点n個に対して、n個のパラメータを持つモデルがフルモデルと呼ばれていたのだから、モデルが持てるパラメータの数は最大n個である。つまり、このn + m個のパラメータを持つモデルは少なくともGLMの範疇で推定することが不可能である（GLMでは無理というだけで、世の中的にはパラメータのほうが多くても、モデル自身がパラメータの縮減を行いながら推定を行うことは、機械学習の分野で普通に行われている）。

　じゃあ、考慮するべき区画ごとにrが決まっている場合、つまり区画を固定効果のように扱うなら、過剰なパラメータの問題も避けられる。しかし、この状況でも区画をランダム効果として考えたほうがお得な場面がある。例えば、施肥量の異なる植木鉢に1個体の植物を植え付け、その植物の種子生産を測ることを考える。かなりの数を測ることにしたので、1か所にまとめて置くことはできなかった。そこで、いくつかの場所に分けて、置いておくことにしたとしよう。まとめておいた場所をA区、B区……と名付けておく。当然、特定の区画に特定の施肥量の植木鉢が固まらないように、各区画に均等に配置したとする。

では、いざ、種子生産を測ってみたとすると、もし区画を区別せずにひとまとめのデータにしたらおそらく過分散になっている。というのも、A区、B区……と分けておいていたのであれば、区画内の管理は非常に似ているけれども、区画間の管理は微妙に違う、ということが起こりえる。区画内は比較的厳密にポアソン分布に従っていたとしても、区画間の違いは間違いなく種子生産に影響を与えて、ばらつきを過剰にさせるだろう。けれども、もはやどんな違いが具体的に存在したのかはわからないわけだ。

そこで区画の効果を説明変数に加えることを考える。シンプルに考えれば、区画の効果をカテゴリカル変数にして組み込むことを考えるだろう。

ここで区画ごとにrB、rC……というパラメータを推定するので、パラメータ数は抑えられている。もちろん、こう考えても絶対に間違いだとは言い切れない。けれども、この解析において、一番興味があるのは施肥量と種子生産の関係だ。区画分けはやむを得ず行ったわけで、区画ごとの違いにはあまり興味もない、そもそも区画間の違いが何に起因するのかわからない。こんな状態では、区画の効果を再現することもできないだろう。そういう意味で、固定効果のように操作した効果というよりも、再現不能なばらつきともいえる。

　では、こう考えてはどうだろう。区画ごとに違いを生み出すパラメータを個別に推定するのではなく、区画ごとの違いは何か確率分布に従って生成されるという縛り（仮定）を設けるのである。この考え方は至極妥当だ。いくら区画ごとに種子生産数が違っていて、何か区画ごとに違っていたとしても、大筋は似ているはずだからだ。全然似ていないのであれば、明らかに考慮するべき事項を見逃している。そして、こう考えることで、自由度の過剰な消費を抑えることができる。

　以前、オフセット項の時も紹介したが、自由度の大きさはパラメータの推定精度に大きな影響を与える。自由度が大きいほど、推定精度は良いものとなる。モデルの自由度はサンプルサイズを推定するパラメータ分だけ差し引く。つまり、過剰にパラメータを推定するほど、自由度は無駄に低下することになるわけである。

　さて、区画がs個あり、これを固定効果と考えるのだったら、s-1個のパラメータを推定する必要がある。これは個別に推定していると考えていることになっている。サンプルサイズnだとして、施肥量と区画を固定効果で考えたときの自由度はn - 2 - (s-1) = n - s - 3である。対して、もし、s個の区画が互いに似ていて、正規分布に従っていると考えられれば、正規分布のばらつきを決めるパラメータσ、この一つのパラメータを推定するだけでいい。区画をランダム効果と考えるなら自由度はn - 2 - 1 = n - 3となって、区画の個数まるまま節約される。

　このように、固定効果をランダム効果をうまいこと切り分けることで、一番興味のある現象を精度よく推定することが可能となるわけだ。だからこそ、初めにランダム効果のことを「興味のない変数」と表現したのである。

「傾き」にもランダム効果を入れられる

　実は今までの説明は、ランダム効果を切片項だけに入れた場合のみを考察していた。また、GLMMのモデル式に立ち返る。

ランダム効果を示す値の位置を変えれば、即座に切片項に影響を与えている変数であることがわかるだろう。一方で、このモデル式では「傾き」はランダム効果として考えている効果、例えば区画間で同じであることを仮定している。

　モデル式中の「傾き」にもランダム効果を入れることが可能である。それは以下のようにする。

発想は切片項に対するランダム効果と同じで、「傾き」は互いに似ているが、ある確率分布に従ってばらつく、と考えているだけに過ぎない。上記の施肥と区画の効果を例にとれば、区画ごとに施肥の効果が違っている状況に相当する。これは、性質としては、施肥と区画の交互作用を考えているような感じに近いのだが、パラメータ数としては区画を固定効果としたときに比べてずいぶん少なくてよく、自由度の節約に貢献している。

GLMが良いのか、GLMMが良いのか？

　途中で、混合分布の形が自明なものがあることを紹介した。例えば、ガンマ分布とポアソン分布の混合分布の負の二項分布である。こういう時、GLMの確率分布を負の二項分布に指定して解析することも可能である。では、過分散データがあるとき、負の二項分布を仮定したGLMと、ガンマ分布（あるいは正規分布）をランダム効果としたポアソン分布を仮定したGLMM、どちらで解析するべきなのだろうか。これには決まった解答はない。けれども、指針はある。

　もし、ランダム効果に相当する要因、例えば、植木鉢の置き場所のようなものの記録が一切ないにもかかわらず過分散なデータに出会ってしまえば、それはもうランダム効果を入れようがない。それゆえに、「負の二項分布のGLMを行わざるを得ない」ということができるだろう。もし、ランダム効果に相当する要因が解析に組み込めるなら順当にGLMMを行うべきだろう。

　GLMにおいて、確率分布を正規分布やガンマ分布にしたとき、実は過分散という概念は存在しない。なぜなら、これらの確率分布は平均とは別に分散の値を決めることができるからだ。この観点から行くと、正規分布やガンマ分布に従うデータは必ずしもGLMMが必要がないようにも感じられる。しかし、一番初めに問題提起したように、解決するべきデータの性質はなにも過分散だけではない。過分散はそれらの問題の一つの現れ方なのだ。GLMにするか、GLMMにするか、という問題は、データの構造に強く依存している問題なのである。

固定効果とランダム効果を識別する

　久保先生のページが役立つ。後日、自分なりに考えをまとめてみよう。