About connecting the dots.

data science related trivial things

1からみなおす線形モデル (2) - 一般化線形モデル

はじめに

ということで,前回で駆け足ではありますが一般線形モデルについてざっと説明しました.そして今回は一般線形モデルについてみていきます.

一般線形モデルの限界

一般線形モデルについて,概要は前回説明しましたが,ここで一般線形モデルの定義を明確にしておきたいと思います.といっても土居さんのレジュメを簡単に言い換えただけですが,以下の通りになります.

y の平均がパラメータの一次結合で表わされる
・説明される変数 y が正規分布に従う

パラメータの一次結合というのは,定数倍したパラメータ同士を足し算したもの,という意味です.そのまま,前回にも出てきた,以下のような式のことですね.こういう式を使ったものを一般線形モデルと呼びましょうと,まぁそれだけのことです.

y=\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}

なので,ここで大事なのは後者の,y が正規分布であるという点です.

つまり,全然正規分布ではないような変数については,それを予測するような一般線形モデルをつくってはいけないということです.いやまぁ,モデル自体はRさんが計算してくれるので,作ること自体はできるんですけどね... 作ることができるからといって,作ったものが妥当であるとは限らない,ということです.

具体例として,日本人の平均年収を年齢から予測する,という線形モデルを考えてみます.しかしながらパラメトリックとノンパラメトリック - About connecting the dots.で書いたように,日本人の平均年収は明らかに正規分布ではないわけです.なので,この場合に通常の線形モデルを当てはめるべきではありません.

それ以外にも,もっと極端な例として,勉強時間と睡眠時間から受験の合否を予測する,といったモデルを考えてみます.これなんて2値変数なので,どうがんばっても正規分布とかしようがないわけです.

では,このような場合にはどうすればよいのでしょうか.ということで上記の問題点を克服したのが一般化線形モデルです.

一般化線形モデル

先ほどと同じように,一般化線形モデルの定義をまとめると,以下の通りになります.

・説明される変数 y が正規分布以外の分布でもよい(例えば二項分布やポワソン分布などが当てはまります*1
・y の平均(や確率)をある関数で変化させたら、パラメータの一次結合で表わされる

つまり,正規分布以外の分布に従うyを予測するための線形モデルが,一般化線形モデルということになるわけです.とはいっても,先ほどの例のような2値変数を扱えるようにするためには,実際には新しいワザを使う必要があります.それが「リンク関数」と「線形予測子」です.これについて説明するために,先ほどの受験の合否の例をもう一度持ってきましょう.

線形予測子

受験には,合格と不合格の2値しかないわけですが,仮に数学と英語の偏差値が全く同じ人100人が入試を受けたとして,全員が合格したり,逆に全員が不合格だったり,ということはあまりないでしょう.それなりに偏差値が高い100人の場合には,100人のうち80人が合格するでしょうし,偏差値が低い100人であれば10人しか合格しないかもしれません.つまり,この合格不合格というのは,確率分布として表すことができます.一般的にコインの裏表のような2値のバラツキを表す分布として,2項分布というものがあります.高校の数学の教科書にも載っていあるあれですね.

つまり合否の2値し方取らなかった予測変数yについて,合格確率pというふうにとらえることで,その合格確率が2項分布に従うと考えることができるようになるわけです.そしてそれによって,線形モデルに載せられそうになってきました.実際に,数学と英語の偏差値から合否を予測する,というモデルを考えたときに,「数学と英語の偏差値から」の部分は以下のような線形結合で洗わせるわけです.

z=\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}

このzを線形予測子という名前で呼びます.この線形予測子を使って,合格確率を予測すればいいわけです.ですが,このままだと問題があるのはすぐにわかるかと思います.上記のzの式だと,原理的には−∞〜∞の値を取りうるわけです.これに対して合格確率pは,当然のことながら0〜1の値しか取り得ないわけです.なので,以下のような式はそもそもおかしい,ということになります.

p=\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}

この問題を解決するためのものが,リンク関数というやつです.

リンク関数

今の問題は,先ほどの線形予測子zがpの取りうる範囲を超えてしまうことなわけです.これを解決するためには,zになんらかの変換をかけて,pの取りうる0-1の範囲に収まるようにしてしまえばいいわけです.具体的には,以下のような変換式を使えば,zが−∞〜∞の値まで変化したとしても,変換した結果は0-1の範囲になります.この変換式のことをリンク関数と呼びます

\frac{1}{1+exp(-z)}

まとめ

以上から,ようやく2項分布に従う2値変数yを予測する線形モデルを作ることができました.

p=\frac{1}{1+exp(-(\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}))}

ただ,このままだと右辺がややこしくてわかりづらいので,式変形してあげて次のような形に直します.これでようやく見慣れた感じの式になりましたね.

log\frac{p}{1-p}=\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}

ということで,「線形予測子」と「リンク関数」を使うことで,見慣れた一般線形モデルに落とし込むことができるようになったわけです.

Rでの実例

ということで,ようやくですが,Rでのサンプルコードに進みましょう.

> res <- c(0, 1, 0, 0, 1, 0, 1, 0, 1, 0)
> math <- c(36, 71, 37, 49, 64, 36, 69, 48, 51, 43)
> eng <- c(30, 72, 52, 40, 54, 43, 68, 38, 59, 47)
> d <- cbind.data.frame(res, math, eng)
> d

   res math eng
1    0   36  30
2    1   71  72
3    0   37  52
4    0   49  60
5    1   64  54
6    0   36  43
7    1   69  68
8    0   48  38
9    1   51  59
10   0   43  47

> fit <- glm(res~math+eng, data=d, family=binomial(link=logit))
> summary(fit)

Call:
glm(formula = res ~ math + eng, family = binomial(link = logit), 
    data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4614  -0.4044  -0.2253   0.1280   1.6059  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -17.1720    11.6290  -1.477    0.140
math          0.1834     0.1495   1.227    0.220
eng           0.1472     0.1528   0.964    0.335

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 13.4602  on 9  degrees of freedom
Residual deviance:  6.4701  on 7  degrees of freedom
AIC: 12.47

Number of Fisher Scoring iterations: 6

合否と数学と英語の3変数からなるデータセットを用意して,一般化線形モデルのメソッドであるglm()を使います.記法に関してはだいたい一般線形モデルと同じですが,先ほど説明したように,yに当てはめる分布を指定することができます.ここでは2項分布を当てはめているので,binomialになります.またその中でリンク関数を指定できますが,ここではロジスティック回帰を行うため,logitを指定しています*2

その結果が下に示した通りですが,ここではサンプル数が少なすぎるため,有意確率は数学が0.22で英語が0.34とさほど高くはありません.ですが係数はどちらもプラスで,つまりどちらも偏差値が高いほど合格しやすいということを表しています.

なお,重回帰における決定係数のようなものとして,以下の式で表される疑似決定係数(Pseudo R-squared)という指標があります.この疑似決定係数にはいろいろな算出式がありますが,BaylorEdPsychパッケージを入れることで一括で算出することができます(指標の詳細についてはこちらを参照してください).いずれも高い値を示していますが,なにぶんサンプル数が少ないので,そこまで意味のある値になっているとはいえないでしょう.

> install.packages("BaylorEdPsych")
Installing package(s) into '/Library/Frameworks/R.framework/Versions/2.15/Resources/library'
(as 'lib' is unspecified)
trying URL 'http://cran.rstudio.com/bin/macosx/leopard/contrib/2.15/BaylorEdPsych_0.5.tgz'
Content type 'application/x-gzip' length 111959 bytes (109 Kb)
opened URL
==================================================
downloaded 109 Kb

tar: Failed to set default locale

The downloaded binary packages are in
	/var/folders/69/4pf7w3g930l6n3d7tz04yxr00000gn/T//RtmpZsQY9N/downloaded_packages
> library(BaylorEdPsych)
> PseudoR2(fit)
        McFadden     Adj.McFadden        Cox.Snell       Nagelkerke McKelvey.Zavoina           Effron            Count        Adj.Count              AIC    Corrected.AIC 
       1.0000000        0.4056567        0.7397268        1.0000000        0.9999509        1.0000000        1.0000000        1.0000000        6.0000000       10.0000000 

ということで,今回は一般化線形モデルについてざっくりまとめてみました.今回の例では,2項分布にしか触れませんでしたが,これ以外にもカウントデータによく用いられるポワソン分布なんかも,一般化線形モデルとしてはよく使われます(family=poissonで指定できます).

さて,この次は一般化線形混合モデルについてみていきたいと思います.

*1:とはいっても,どんな分布でもいいわけではなくて,厳密には「指数型分布族」に含まれる分布)でなければいけないそうです

*2:分布にbinomialを指定した場合には,何もしなくてもデフォルトでリンク関数にlogitが適用されます