データ解析
第四回「線形回帰分析の拡張：一般化線形モデル」
鈴木 大慈
理学部情報科学科
西八号館 W707 号室
s-taiji@is.titech.ac.jp

1 / 29

休講情報

7/7 は休講

2 / 29

今日の講義内容
一般化線形モデル (glm)

3 / 29ガウス・マルコフモデルの限界
ガウス・マルコフモデル:

y = β ⊤ x + ϵ (ϵ ∼ N(0, σ 2 )).
y は各 x でガウス分布．
y の期待値は説明変数 x に対して線形.
しかし，従属変数は離散値 (整数) を取ったり，非負値の制約があったり，説明変数と非線形な関係であったりする．
→ 「線形・正規分布」以上のことをしたい．
それを可能にするのが一般化線形モデル．
線形判別分析 ともつながる．

4 / 29ガウスマルコフモデルのイメージ

5 / 29このようなデータでは．
．？

6 / 29線形回帰してみる

あまり良くなさそう．

7 / 29一般化線形回帰こちらの方が良さそう．
→ 一般化線形モデル

8 / 29一般化線形モデルの基本形 (重要)
二つの構成要素あるパラメータ θ によって特徴づけられる分布 Pθ (Y ).

u = β ⊤ x に対して Y の分布 Pθ のパラメータ θ を定める関数 θ = g (u):g −1 のことをリンク関数と呼ぶ．


一般化線形モデル：

yi ∼ Pθ=g (β ⊤ xi )(i = 1, . . . , n).

例：ガウスマルコフモデルでは Pθ = N(θ, σ 2 ), g (u) = u.

Pθ=g (β ⊤ x) の密度関数を p(y |g (β ⊤ x)) と書くと，対数尤度はn∑log(p(yi |g (β ⊤ xi ))),i=1となる．
9 / 29一般化線形モデルの例 ：二値判別ベルヌーイ分布とロジットリンク関数:
Y = 0 または 1, Berθ (Y ) = θY (1 − θ)1−Y (θ ∈ (0, 1), Y = 0, 1).

1
1 + exp(−β T x)
1
1P(Y = 0|X ) = 1 −=.
T
1 + exp(−β x)
1 + exp(β T x)
P(Y = 1|X ) =すなわち，1
θ = g (β ⊤ x) = 1+exp(−β,⊤( x))⊤−1θ
β x = g (θ) = log 1−θ : ロジット関数.

(yi ∼ Ber θ =1
1 + exp(−β ⊤ xi )).

ロジスティック回帰と呼ばれる．
例えば，顔認識では x が画像で，y = 1 のときにその画像が顔画像，y = 0 のときに顔画像以外というように用いる．
10 / 29

0.6
0.4
0.00.2logistic(x)0.81.0

一般化線形モデルの例 ：二値判別（つづき）−4−2024u

Figure: ロジスティック関数 g (u).

u が正に大きければ Y = 0 になりやすく，u が負に大きければ Y = 1 になりやすい．
前のページの記法で言うと，β ⊤ x の大小で Y の出方が決まる．
11 / 29

一般化線形モデルの例 ：ロジスティック回帰の尤度関数{1(yi = 0)と書き換える．
−1 (yi = 1)5

簡単のため，ỹi =i=1=n∑

11 + exp(−ỹi β ⊤ xi )()− log 1 + exp(−ỹi β ⊤ xi ) .

i=13)2(loglog(1+exp(−x))=1i=1n∑4log(p(yi |g (β ⊤ xi )))0n∑−4⇒凹関数．最適化しやすい．

−2024xlog(1 + exp(−x)) のグラフ

12 / 29

一般化線形モデルの例 ：二項分布0.30

二項分布とロジットリンク関数:( )Binθ|N (Y ) = YN θY (1 − θ)N−Y (θ ∈ (0, 1), Y = 0, 1, . . . , N).
()
1yi ∼ Bin θ =|N
1 + exp(−β ⊤ x)()1
θ⊤−1g (β ⊤ x) = 1+exp(−β(θ) = log 1−θ: ロジット関数.
⊤ x) , β x = g

0.15
0.10
0.050.00

特に N = 1 の場合は二値判別.

Bin(y)

0.200.25

theta=0.5
theta=0.2051015y13 / 29一般化線形モデルの例 (離散つづき)
ポアソン分布と対数リンク関数:Y −θPoθ (Y ) = θ Ye ! (θ > 0, Y = 0, 1, 2, . . . ).

yi ∼ Po(θ = exp(β ⊤ x))
g (β ⊤ x) = exp(β ⊤ x), g −1 (θ) = log(θ) = β ⊤ x: 対数リンク関数.
β ⊤ x が負の値をとっても大丈夫！
⊤(exp(β ⊤ x))Y e − exp(β x)P(Y |β x) =,Y!
log(P(Y |β ⊤ x)) = Y (β ⊤ x) − exp(β ⊤ x) − log(Y !).
⊤対数尤度関数:n∑i=1log(P(yi |β ⊤ xi )) =n∑[yi (β ⊤ xi ) − exp(β ⊤ xi ) − log(yi !)].
i=1

⇒ β に関して凹関数．最適化しやすい．指数型分布族はすべて凹関数になる．
14 / 29

一般化線形モデルの例 (離散つづき)
負の二項分布とロジットリンク関数:( +k−1) kNbθ|k (Y ) = Y k−1
θ (1 − θ)Y (θ ∈ (0, 1), Y = 0, 1, 2, . . . ).
()1
yi ∼ Nb θ =|k
1 + exp(−β ⊤ x)()1
θ⊤−1g (β ⊤ x) = 1+exp(−β,βx=g(θ)=log⊤ x)1−θ : ロジット関数.

β ⊤ x が −∞ から ∞ までの値をとることで，θ は (0, 1) 区間を動く．
)k ()Y()(Y +k −1
11P(Y |β ⊤ x) =,1 + exp(−β ⊤ x)
1 + exp(β ⊤ x)
k −1( +k−1)log(P(Y |β ⊤ x)) = log Y k−1− k log(1 + exp(−β ⊤ x)) − Y log(1 + exp(β ⊤ x)).
対数尤度関数:n∑i=1log(P(yi |β ⊤ xi )) =n∑[−k log(1 + exp(−β ⊤ xi )) − yi log(1 + exp(β ⊤ xi ))i=1( +k−1)+ log Y k−1].
15 / 29

0.00

0.00

0.05

0.05

0.10

0.10

02468100.15

Nb(y,k=3)

0.15Po(y)

0.20

0.20

0.25

0.25

theta=2
theta=6
theta=0.3
theta=0.65

y024y6810

ポアソン分布と負の二項分布

16 / 290.300.30

一般化線形モデルの例 (連続)ガンマ分布と逆数リンク関数:
α−1 −Y /θeΓθ|α (Y ) = Y Γ(α)θ
(θ > 0, Y ≥ 0).
αyi ∼ Gamma(θ =1|α)β⊤xg (β ⊤ x) = β ⊤1 x , β ⊤ x = g −1 (θ) = θ1 : 逆数.
R では形状パラメータ α も同時に推定される．
正規分布と恒等リンク関数:()−θ)2N(Y |θ = β ⊤ x, σ 2 ) = √ 1 2 exp − (Y2σ(θ Y ∈ R).
22πσyi ∼ N(θ = β ⊤ x, σ 2 )g (β ⊤ x) = β ⊤ x, β ⊤ x = g −1 (θ) = θ: 恒等写像.
これはガウスマルコフモデル．
※ここで挙げた分布とリンク関数の組は，よく使われる例として挙げただけで，この組み合わせでなくてはいけないということはない．
17 / 29一般化線形モデルの推定
前のスライドに書いたように β の対数尤度はn∑log(p(yi |g (β ⊤ xi ))),i=1となる．ここで，ℓ(yi , β ⊤ xi ) := log(p(yi |g (β ⊤ xi ))) と書くと，β の最尤推定量は

β̂ = arg max
βn∑ℓ(yi , β ⊤ xi ),i=1となる．ℓ をロス関数と言ったりもする．最適化は汎用最適化ソルバーなどを用いる (最急降下法，準ニュートン法など)．

pθ や g を介さずに直接ロス関数 ℓ を設計することもある．

18 / 29

漸近正規性

ここで，もしモデルが正しかったら，√d⊤−1n(β̂ − β ∗ ) → N(0, Eβ ∗ [∇β ℓ(Y , β ⊤ X )∇⊤β ℓ(Y , β X )|β=β ∗ ] ).
(漸近正規性)

なお，Fisher 情報行列は⊤Eβ ∗ [∇β ℓ(Y , β ⊤ X )∇⊤β ℓ(Y , β X )|β=β ∗ ]n1∑⊤≃
EY |β̂,xi [∇β ℓ(Y , β ⊤ xi )∇⊤β ℓ(Y , β xi )|β=β̂ ]ni=1で近似可能．
→ 信頼区間の構築や検定が可能に．

19 / 29

今回のデモで必要なライブラリMASS
faraway
install.packages(..)
library(..)
でインストール可能．

20 / 29

一般化線形モデルを R で実際に使ってみる
data(gala,package="faraway")
ガラパゴス諸島の 30 の島と亀の種類との関連
7 変数 30 サンプル
Species：その島の亀の種類の数 (従属変数)
Endemics：亀固有種の数 (説明変数)
Area：島の面積 (km2) (説明変数)
Elevation：島の標高 (m) (説明変数)
Nearest：最近隣の島との距離（km）(説明変数)
Scruz：Santa Cruz 島との距離 (km) (説明変数)
Adjacent：近隣の島のエリア（km2）(説明変数)

21 / 29

一般化線形モデル (glm) 基本形
一般化線形モデルの最尤推定:

gala.pm1<-glm(Species~ ., data = gala,family = poisson(link="log"))
family で分布族を決定．link でリンク関数を決定．
ポアソン回帰の場合はリンク関数のデフォルトが log なので link=”log”は省略
可能．
gala.pm1<-glm(Species~ ., data = gala, family = poisson)

22 / 29

family と link 関数の例
family とデフォルトの link 関数の組:
binomial(link = "logit")
gaussian(link ="identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")
negative.binomial(link = "log")
negative.binomial は library(MASS) が必要．

23 / 29

結果の要約要約の表示: summary(gala.pm1).
Coefficients:

Estimate Std. Error z value Pr(>|z|)(Intercept) 2.828e+00 5.958e-02 47.471 < 2e-16 ***Endemics
3.388e-02 1.741e-03 19.459 < 2e-16 ***Area-1.067e-04 3.741e-05 -2.853 0.00433 **Elevation
2.638e-04 1.934e-04
1.364 0.17264
Nearest1.048e-02 1.611e-03
6.502 7.91e-11 ***Scruz-6.835e-04 5.802e-04 -1.178 0.23877
Adjacent
4.539e-05 4.800e-05
0.946 0.34437

係数 β の推定値，標準偏差，Wald 統計量 (z), p-値．
ここで，係数の標準偏差の導出には最尤推定量の漸近正規性を用い，正規分布に近似して求めている (正規分布で近似してしまえば，最小二乗法と同じ理屈).
Wald 統計量は最尤推定量を (推定された) 標準偏差で割ったもの．βi∗ = 0 なる
帰無仮説のもと，漸近的に正規分布に従う．
つまり，Pr(>|z|) が十分小さければ βi∗ = 0 なる仮説は棄却される．
24 / 29

結果の要約 2

Null deviance: 3510.73 on 29 degrees of freedom
Residual deviance: 313.36 on 23 degrees of freedomdeviance (逸脱度) とは，あてはまりの「悪さ」の指標で，D = −2 log(L̂)
で与えられる．すなわち，最尤推定量の対数尤度に −2 をかけたものである．

Residual deviance は D から「サンプル数分のパラメータを使ってあてはめたモデル」の逸脱度を引いたもの．
値が小さければ小さいほど手元にあるサンプルへの当てはまりが良い．モデルが正しければ漸近的に「degree of freedom」と同じ自由度の χ2 分布に従う.
Null deviance は切片のみのモデルの Residual deviance である．
この例の場合，Residual deviance は大きく，ポアソンモデルがそこまで良いとは言えない．
→ 負の二項分布など別の分布を試してみるとどうなるか? (レポート課題)
25 / 29

モデル選択Residual deviance が小さいモデルは手元にあるデータへの当てはまりは良いが，必ずしも予測力が高いわけではない．
予測誤差が小さくなるようなモデルを選ぶ場合，AIC を用いてモデル選択を行えばよい．

glm オブジェクトに対しても step(.) や AIC(.) が使える．

26 / 29

第二回レポート
自分でデータを用意して線形回帰分析を行ってみよ．
※データは何でもよい．ネットに落ちているものでよい．ただし，データの出自を明記せよ．

gala データに，ガウスマルコフモデルおよび負の二項分布モデルを当てはめ，講義で行ったポアソンモデルとあてはまりの良さや予測誤差 (AIC) 等を比較せよ．
solder データでも同様にガウスマルコフモデルおよびポアソンモデルと比較せよ．
（optional）余力があれば他の分布やリンク関数を当てはめてみよ．

27 / 29

レポートの提出方法
私宛にメールにて提出．
件名に 必ず「データ解析第 n 回レポート」と明記し，R のソースコードと結果をまとめたレポートを送付のこと．
氏名と学籍番号も忘れず明記すること．
レポートは本文に載せても良いが，pdf などの電子ファイルにレポートを出力して添付ファイルとして送付することが望ましい (これを期に tex の使い方を覚えることを推奨します)．
出力結果をただ提示するだけでなく，必ず考察を入れること．
提出期限は講義最終回まで．
※相談はしても良いですが，コピペは厳禁です．

講義情報ページ
http://www.ocw.titech.ac.jp/index.php?module=General&
action=T0300&GakubuCD=100&GakkaCD=15&KougiCD=5522&
Nendo=2015&Gakki=1&lang=JA&vid=03
http://www.is.titech.ac.jp/~s-taiji/lecture/2015/dataanalysis/
dataanalysis.html
28 / 29一般化線形モデル参考文献[1] 久保拓弥: データ解析のための統計モデリング入門‐一般化線形モデル・階層
ベイズモデル・MCMC．岩波書店，2012.
[2] J.F. Faraway: Extending the linear model with R. Chapman and HallCRC,
2005.

29 / 29