一般化線形モデル 一般化線形モデル

Table of Contents

前回 はチュートリアルの動機と目標を確認し、 線形モデルを作成しました。 モデルの作成には data.frame, formula, model が必要であり、 例として lm でモデルを作成しました。 しかし線形モデルは応答変数$y$が説明変数$x$に対し直線の関係にあることや、 小数や負の値を取れることを仮定しています。 今回は応答変数が 説明変数と直線の関係にないケース や 様々な分布に従うケース を「一般化線形モデル」で扱います。

この一般化線形モデルを理解するため、 まずは実データに触れながら (i) 線形モデルを使う際の問題 を確認します。 そして (ii) 問題を解決する「リンク関数」 を抑え、 最後に (iii) glm を用いた一般化線形モデル を実際に作成します。

直線の当てはめで起きそうな問題

まずは直線の問題点を見るために xyplot() を使いましょう。 前回のデータ sleepstudy に「回答の正解/不正解(1/0)」の列 Correct を追加しました。 今回は寝不足Days がCorrect を説明するかを検証するので、 復習がてら Days と Correct の関係を formula として f に与えて可視化しましょう。 可視化のしやすさのため、 sleepstudy は 350番の被験者データに限定しています1。

library(lattice); library(lme4) set.seed(43) # data rt = sleepstudy$Reaction logistic = function(x) {1 / (1 + exp(-x))} scale = function(x) (x-min(x))/(max(x)-min(x)) ps = logistic(-5*scale(rt)+2) xs = as.integer(rbinom(ps, n=length(ps), size=1)) sleepstudy$Correct = xs # f = Reaction ~ Days # 前回の範囲 f = xyplot(f, sleepstudy, subset=Subject==350, # 350番に限定 type=c("p", "r") # plot, regression ) # 余裕があれば `lm` を用いた回帰分析も確認 # lm(f, sleepstudy, subset=Subject==350) f = Correct ~ Days xyplot(f, sleepstudy, subset=Subject==350, # 350を指定 type=c("p", "r") # plot, regression ) lm(f, sleepstudy, subset=Subject==350) test_function("xyplot") success_msg("Great job!")
前回は Reaction ~ Days を見ました。今回は Correct を見ましょう。

当てはめられた回帰直線が右肩下がりなので、 Days が増すほど正解が減ることは分かります。 しかし4日目あたりの正答数0.8とはどういうことでしょうか。 また、Days によっては 0–1 の範囲を超えてしまいそうです。 lm を用いた回帰分析をした結果からも 切片は1を超えており Days が0 の時におかしなことになります。 したがって $Y = ax + b$ の当てはめでは 以下の問題が発生しうることを確認しましょう。

  1. $Y$ が小数点を取らないとき(正解数0.8は意味不明)
  2. $Y$ に範囲の制限があるとき(直線だと飛び出す)

GLMの基本要素

整数を取ってほしいのに小数点の値になってしまう、 範囲指定があるときに直線の当てはめをすると範囲外に飛びだす、 などの問題に対応するため、 (1) 確率分布と(2) 線形予測子、(3) リンク関数 の概念を抑えていきます。 どれも最初は抽象的でわかりづらく思えますが、 これが分かればGLMの考え方や結果の見方をマスターできますし、 おそらくRでもSPSSでも問題ないでしょう。 ただ、ここからは話が抽象的になるので、 手を動かしていきましょう。

確率分布

線形モデルの $ax+b$ が直線である以上、 傾きが0でない限りいつかは指定の範囲を飛び出しますし 小数点も取ります。 そこでまず、 問題一つを与えた時の正答数 $Y$ が「二項分布」に従う場合を考えます。 この二項分布は確率 $p$ で現れる事象が $size$回のうちに起こる回数(x)を分布(y=確率)として表現します。 話が抽象的なので手を動かしていきましょう。

問題を正答できる確率 $p$ が0.8だったとして、 $size$ が1の二項分布を考えてみましょう。 これは「1回の試行で成功確率が0.8の時、X回成功する確率はYだ」という分布です。 そもそも1回の試行なので、成功する回数は0か1です。 試しに100回二項分布に0と1を出させてみす。 rbinom(p=0.8, size=1, n=100) の結果を y に格納して、 hist(y) としてみましょう。

y = # edit here hist() # edit here y = rbinom(p=0.8, size=1, n=100) hist(y) test_function("hist") success_msg("Great job!")
use rbinom and hist functions

0が少なく1が多い結果になりました2。 つまり、確率0.8の現象が1回のうちに起こる回数(0/1)は 「1の場合のほうが起こりやすそうだ」という分布になりそうです。 このように $Y$ の分布が確率 $p$ とサイズ $size$ の二項分布に従うことを $Y \sim \textrm{Binomial}(p, size)$ と表現します。 この時点で$Y$が離散かつ制限のある場合を表現できているのですが、 まだ $ax+b$ との関係が不明です。

さて、上の Histogram of y から判断すると $p$ は0.5と0.8、どちらが尤もらしいでしょうか。 おそらく0.8の方が尤もらしいと感じると思います。 つまり「データに対して最も尤もらしい(最尤の)$p$」は 同じデータで複数の$p$の比較し、最尤の $p$ が答えになります。 なお、このような最尤パラメータの推定を「最尤推定」と呼びます。 また推定ごとに $p$ の値はぶれるので、 そのブレはそのまま「$p$ の値のばらつき」と考えられます。

ここでようやく $p=ax+b$ が成り立つと仮定しましょう。 そうすると、データに対して適当に $a$ や $b$ を変えて 比較すれば $a$ や $b$ を最尤推定できます。 そして $a$ がわかればモデリングの目的の一つである 「要因の大きさを知ること」が達成できます。 ただ $p$ は確率なので0–1の値を取るのに対し、 $ax+b$ は直線で範囲を超えてしまうので $p=ax+b$ とは表現できません。 この先送りにされた問題は次の節で考えます。

ここまでをまとめると以下になります。

  1. $Y$の範囲が制限されている、値が離散であるべき、という問題は「確率分布」が解決
  2. 手元のデータから確率分布の$p$は推定可能で $p=ax+b$ が成り立つなら $a$ や $b$ も推定可能
  3. 実際には成り立たないので、推定のためにはもう一つ作業が必要

確率分布と線形予測子

確かに $p=ax+b$ とは表現できないのは事実です。 しかし「$\textrm{logistic}$関数」は $ax+b$ の結果(実数)を必ず 0–1 の間に収めます。 先ほど $Y$ を見て適当に $p$ を動かして 最尤の $p$ を決めるという話をしたのですが、 $p=\textrm{logistic}(ax+b)$が成り立つなら、 適当に $a$ や $b$ を動かしても 範囲を飛び出さずに $a$ や $b$ を決められるはずです3。 そして $a$ がわかればモデリングの目的の一つである 「要因の大きさを知ること」が達成できます。

ここで一旦 $ax+b$ を lp (linear predicor, 線形予測子) として、 lp の結果が -4 から 4 のように 0–1の外になるケースを見てみましょう。 そして logistic が-4から4の値 lp をどのように 0 から1 に変換するかを見ます。 R では lp = -4:4 とすれば lp に -4から4 の整数を格納します。 この lp を logistic に与えた結果を p として plot しましょう4。

library(lattice); library(lme4) logistic = function(x) {1 / (1 + exp(-x))} lp = # -4から4までの整数を作成 (ax+bの結果) p = # logistic()にlpを与えた結果をpに格納 plot(x=lp, y=p, type="l") # "l" for line lp = -4:4 p = logistic(lp) plot(x=lp, y=p, type="l") test_function("plot") success_msg("Great job!")

$\textrm{logistic}$ が実数を 0–1 の間に収めそうなことを 確認できましたか5。 まだ気になる人は logistic(-Inf) や logistic(Inf) も試してみて、 リミットを見てみましょう。 ともかく、これで $ax+b$ の値を 0–1 に収められ確率の$p$ の条件を満たしたので $p=\textrm{logistic}(ax+b)$ と表現でき、 $\textrm{logistic}$ 経由で$ax+b$ のパラメータを最尤推定できます。

ここまでを振り返ると、データ $Y$ から「確率分布の$p$」, それと $p=\textrm{logistic}(ax+b)$ から $a$ と $b$ も最尤推定できそうです。 ただ、そこで推定された値の解釈は困難です。 例えば $b=0$ の時、 「$x$ (寝不足)が1日増えると正答確率$p$ は $\textrm{logistic}(a)$ だけが増える」 と解釈できるのですが、 解釈のたび $\textrm{logistic}(ax + b)$ を計算するのは面倒です。 この問題を解決するのが「リンク関数」です。

リンク関数と解釈

解釈のたびに $\textrm{logistic}(ax + b)$ を計算するのが面倒なら、 その計算をなかったことにする6、 「$\textrm{logit}$関数」を考えてみましょう。 つまり、解釈の難しかった $p = \textrm{logistic}(ax+b)$ の両辺に $\textrm{logit}$関数を当てると、 $\textrm{logit}(p) = ax+b$ と書き換えられます。 そうすると、$x$ (寝不足)が1日増えると $a$分だけ $\textrm{logit}(p)$ が増すと言い換えられます。 もし $\textrm{logit}(p)$ 自体が解釈しやすければ $\textrm{logit}(p)$ を計算せず、そのまま解釈できます。

下では先程と同様、$ax+b$ を lp として -4 から 4 を与え、 lp を logitic に入れた結果を p としました。 上に述べたように logit(logistic(lp)) == lp は成り立つのでしょうか。 また、 logit(p) == lp も正しいのでしょうか。 そこで計算上の誤差を 0.000001 としたとき、 上のそれぞれが成り立つなら両辺の差は計算上の誤差より小さくなる (左辺-右辺 < 0.000001) はずです。 以下は実行するだけなのですが、確認しましょう。

library(lattice); library(lme4) logistic = function(x) {1 / (1 + exp(-x))} logit = function(p) log({p / (1-p)}) lp = -4:4 p = logistic(lp) # lpを与えた結果をpに格納 # **以下は実行してTRUEを確認するだけ** logit(logistic(lp)) - lp < 0.000001 logit(p) - lp < 0.000001 lp = -4:4 tol = 0.000001 # 計算上の誤差 p = logistic(lp) # にlpを与えた結果をpに格納 logit(logistic(lp)) - lp < tol logit(p) - lp < tol success_msg("Great job!")

$\textrm{logistic}$ の結果を $\textrm{logit}$ が打ち消せること、 また $\textrm{logit}$ を $p = \textrm{logistic}(ax+b)$ の両辺に当てると、 等号の関係を変えないまま $\textrm{logit}(p) = ax+b$ となることを確認できましたか。 これで $x$(寝不足)が1日増えると$a$分だけ $\textrm{logit}(p)$ が増すと言えます。 $\textrm{logit}(p)$ の実体さえ分かれば推定したパラメータの意味が 解釈しやすくなります。

正答確率を $p$ としたとき、$\textrm{logit}(p)$ は $\textrm{log}\frac{p}{1-p}$ と表現できます。 みためはややこしいですが、 $\frac{p}{1-p}$ の部分は「オッズ」と呼ばれており、 要は「正答確率 $p$ が誤答確率 $1-p$ に比べてどのくらい大きいか」を表現しています。 例えば $p=1-p$ なら$1$, $p$ が倍なら$2$になります。 つまり、ロジット関数はこのオッズの対数を取った値になります。

対数を取ると分かりづらいので、 両辺に exp を与えて対数を落とします7。 つまり exp(logit(p)) = exp(a*x + b) とすると p/(1-p) = exp(a*x) * exp(b) になります。 したがって、もし最尤推定の結果が $a=2$, $b=0$ だったら $x$ が1増えるとexp(2*1) * exp(0)を計算し、結果は7なので オッズ(正答確率/誤答確率)が7倍になると解釈できます。

ここまでの流れをまとめると以下になります。

  1. 確率分布とは観測したデータ $Y$ が従う分布で、 $Y$の範囲が制限されている、値が離散であるべき、という問題を解決する (例: 正答数や5段階評価なども二項分布)
  2. 手元のデータから確率分布のパラメータ (例: 成功確率$p$) は最尤推定でき、 また $p = \textrm{logistic}(ax+b)$ と変換を挟めば $a$ や $b$ も最尤推定可能
  3. $\textrm{logistic}(ax+b)$の計算はめんどうだが $\textrm{logit}(p) = ax+b$ とリンク関数で書き換えることで 最尤推定した $a$ の効果を容易に解釈可能

一つ前の「線形モデル」が $Y = ax+b$ だけで表現できていたのですが、 それと比べると新たに確率分布とリンク関数が加わっています。 ただ、実際にモデリングするときに気をつけるのは確率分布だけで大丈夫で、 リンク関数は「知っていると解釈のときに便利」程度です。 今回はこの確率分布に「二項分布」を仮定しましたが、 これのどこが「一般化」なのでしょうか。

実は一つ前の「線形モデル」は確率分布に正規分布を仮定していました。 GLMは確率分布に正規分布だけでなく今回のように二項分布も仮定できます。 他にも0以上の整数の分布(ポワソン分布)や、 0以上で偏っている分布(対数正規分布)など、 「指数型分布族(exponential family)」に属す様々な分布を扱えます。 このように $Y$ が従う様々な確率分布を一般化線形モデルは扱えるのです。 この扱える $Y$ の拡張こそが「一般化」になります。

それでは最後に glm を使ってモデリングして 一般化線形モデルの説明を終えます。

glm を用いた一般化線形モデル

GLM でもlmと同様に data.frame, formula, model が必要で、 formula にはデータの列の関係を記述します。 glm が lm と違う点は $Y$ が従う確率分布を glm内で family=binomial のように指定する点です。 なお、リンク関数は family=binomial とした時点で 自動的に logit が指定されるので定義は不要です。 以下に Correct と Days の関係を formula に定義し、 family の指定をしましょう。

library(lattice); library(lme4) set.seed(43) # data rt = sleepstudy$Reaction logistic = function(x) {1 / (1 + exp(-x))} scale = function(x) (x-min(x))/(max(x)-min(x)) ps = logistic(-5*scale(rt)+2) xs = as.integer(rbinom(ps, n=length(ps), size=1)) sleepstudy$Correct = xs f = m = glm(f, sleepstudy, family=) summary(m) f = Correct ~ Days m = glm(f, sleepstudy, family=binomial) summary(m) test_function("summary") success_msg("Great job!")
前回は Reaction ~ Days を見ました。今回は Correct を見ましょう。

結果を格納した m のサマリーの Coefficients を見ると、 以下のようになっています。 切片(Intercept) に対する Estimate は 0.83, 要因Days に対する Estimate は -0.16となっています。 線形モデルの時と同様、要因($Days$) と同じ行にある Pr(>|t|) の値を見ると、 0.05 より小さいので 「その傾きは統計的に意味のあるものなんだよ」と安全に主張できます。

> Coefficients:
>              Estimate Std. Error z value Pr(>|z|)  
> (Intercept)  0.83782    0.29216   2.868  0.00413 **  
> Days        -0.16017    0.05442  -2.943  0.00325 **

また、上のテーブルは リンク関数を使うと logit(p) = -0.16*Days + 0.83 と表現でき、 つまりは p/(1-p) = exp(-0.16*Days) * exp(0.83) です。 ここで$x$ が1の時の $a$ の効果を知りたいときは $b=0$ として、 exp(-0.16*1) * exp(0) を考えます。

a = -0.16 b = 0 odds = exp() * exp() # p/(1-p)つまりoddsを計算 print(odds) a = -0.16 b = 0 odds = exp(a*1) * exp(b) print(odds) test_function("print") success_msg("Great job!")
前回は Reaction ~ Days を見ました。今回は Correct を見ましょう。

結果は 0.85 となっているので、 Days が1増えると正解するオッズは0.85倍になる、 つまりオッズが下がるという意味になります。 言い換えると、Days が1増えると 正解する確率が0.85倍になる、ということです。

まとめ

一般化線形モデルでは確率分布を使って $Y$ の制限を実現し、 リンク関数を使って解釈を容易にしていました。 実際にモデルを使う時の変化は lm ではなく glm を使い、 確率分布を family に指定するだけです。 その過程には $\textrm{logistic}$関数や確率の話が出てきましたが、 これを抜かすといわゆる「ブラックポックス」な統計になってしまいます。

このページをTokyoRというSlackのコミュニティに投稿したところ、 「GLMとロジスティック回帰の違いがわかりづらい」という声がありました。 今回の話は「ロジスティック回帰」と呼ばれるもので、 これは GLM の一種です。 GLMにはロジスティック回帰や「ポアソン回帰」などの種類があるのですが、 基本の考え方は今回と同じです。 応答変数の性質、従う確率分布、そして応答変数と説明変数の関係を考え、 ベストなGLM (確率分布とリンク関数) を選択しましょう。

応答変数の性質 応答変数が従う確率分布 リンク関数の候補
0以上の整数 二項分布(binomial) logit
0以上の整数 ポアソン分布(poisson) log
0以上の実数 ガンマ分布(gamma) log
実数 正規分布(gaussian) identity

正直、Reaction も0以下にはならないので 「実数」というと言い過ぎなのですが、 そこらへんはモデルの厳密さと複雑さとの相談になります。 ともかく、次はチュートリアルの最終目標である 構造的なノイズ の扱いについてです。

脚注

  1. formula を使う部分では大抵 subset オプションがある。 ↩

  2. 一度 R Console に y と打ってみると rbinom(p=0.8, size=1, n=100) の結果を見れます。 ↩

  3. 厳密にはf(x)が短調増加の関数でないといけません。 また、aやbを動かすイメージがつかない場合は一旦aを0だと考えると、 ただの p=f(b) になります。bが増えるときは必ずpが増えるので、 それはpを動かして調整しているのと同じことになります。 ↩

  4. ちなみに logistic は以下のように定義できます。

    logistic = function(x) {1 / (1 + exp(-x))}

    今回の目標から考えると関数の作成は蛇足ですが、 Don’t Repeat Yourself (DRY) という概念も 本当はお伝えしたい気持ちです。 ↩

  5. ちなみに sigmoid の意味は「sっぽい」です。 アンドロイドとか ALife のボイド(boid: bird-oid) も同じ語源です。 ↩

  6. つまり $\textrm{logit}(\textrm{logistic}(ax+b)) = ax+b$ となる。 ↩

  7. logとexpの関係は 上の R Console で exp(log(2)) とすれば logisticとlogitの関係と同じだとわかります。 ↩