このような場面で便利なのがオフセット項わざ

6.5-6.9 GLM
FM13005 大楠拓也
目次
6.5 交互作用項の入った線形予測子
6.6 割算値の統計モデリングはやめよう
• 6.6.1 割算値いらずのオフセット項わざ
6.7 正規分布とその尤度
6.8 ガンマ分布のGLM
6.9 6章のまとめ
使用するデータがあるURL:
http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
2
6.5 交互作用項の入った線形予測子
• 線形予測子”𝛽1 + 𝛽2 𝑥𝑖 + 𝛽3 𝑓𝑖 ”を複雑化して、交互作用項を追加
• この例題の場合の交互作用とは
 植物の体サイズ”𝑥𝑖 ”と施肥処理の効果"𝑓𝑖 ”の「積」の効果である
• 交互作用項をいれた線形予測子は以下のようになる
logit(𝑞𝑖 )= 𝛽1 + 𝛽2 𝑥𝑖 + 𝛽3 𝑓𝑖 + 𝛽4 𝑥𝑖 𝑓𝑖
参考URL(14~):
http://hosho.ees.hokudai.ac.jp/~kubo/stat/2013/e/kubostat2013e.pdf
3
• 施肥処理の要因は「無処理(C)」「肥料をやる(T)」の２水準
この因子型の説明変数𝑓𝑖 : Cの時は0、Tの時は1となることにする
• 交互作用項𝛽4 𝑥𝑖 𝑓𝑖 は単純に「𝑥𝑖 と 𝑓𝑖 の積に係数𝛽4 の値を掛けたもの」
と考えればよい
例えば、交互作用が大きな影響を持つ場合、図6.8のように平均生存種子数
のサイズ𝑥𝑖 依存性は、施肥処理𝑓𝑖 によって大きく変わる
4
図6.8
• 交互作用項が大きいので、
サイズ依存性が施肥処理
によって大きく変わる場合
の一例
5
生存種子数データ(図6.2)の交互作用を推定
glm(cbind(y,N-y)~x*f,family=binomial,data=d)
• モデル式の右辺”x*f”は”x+f+x:f”を省略する記法
 省略していないモデル式の”x:f”項が交互作用を表している
• パラメータ推定を行うと以下の結果が得られる
Coefficients:
(Intercept)
x
fT
x:fT
-18.52332 1.85251 -0.06376
0.21634
※d <- read.csv("data4a.csv")が必要
※” でエラーが出た場合は自分で ” を入力
Degrees of Freedom: 99 Total (i.e. Null); 96 Residual
Null Deviance: 499.2
Residual Deviance: 122.4
AIC: 273.6
6
• 推定結果と、6.4.4項でAIC最良モデルとして選ばれたx + f モデルfT
の係数の推定値(fTの推定値は2.02,p.123)が異なるように見える
• 施肥処理をしなかった場合(C)
logit(𝑞𝑖 )= −18.5 + 1.85𝑥𝑖
• 施肥処理をした場合(T)
logit(𝑞𝑖 )= −18.5 − 0.0638 + (1.85 + 0.216)𝑥𝑖 = −18.6 + 2.07𝑥𝑖
• しかし図示した図6.9を見ると、モデルの予測としてはほとんど変化し
ていないことが分かる
7
図6.9 交互作用の有無を調べる図示
x + f モデル
x * f モデル
8
• x + f モデルのAICは272だったが、交互作用を入れると274に悪化した
• この例題に関して言えば、交互作用項を追加してもモデルがわかり
にくくなっただけで、予測能力は何も改善されなかったと結論できる
• 交互作用項はその係数だけ見ていても解釈できない場合が多く、
推定結果にもとづく予測などを図示する必要がある
• 交互作用項についての利用で注意すべき点
「むやみに交互作用項をいれない」こと
9
むやみに交互作用項を入れてはいけない
1. 説明変数が多い場合には交互作用項の個数が「組み合わせ論的
爆発」で増加してパラメーター推定が困難になる
2. 交互作用項を入れてみたものの、それが何を表しているのか解釈
できなくなる事例をよく見かける
※交互作用の推定とは、観測データから複雑なパターンの抽出を狙ったもの
実際の観測データに基いてそれを特定するのはなかなか困難な場合がある
10
• 現実のデータにGLMをあてはめた場合、交互作用項を多数含んだ
統計モデルのAICが最良になることがある
• しかしながらこれは交互作用項の効果を過大推定している
言い換えれば、ニセの交互作用でつじつま合わせをしている可能性もある
• 現実のデータでは、説明変数では説明できない「個体差」「場所差」
によるばらつきが発生する
• これらを考慮していないGLMをあてはめた場合、過度に複雑なモデ
ルが選ばれる傾向がある
参考URL(クリスピー(チキン)を例にしている):
http://kogolab.chillout.jp/elearn/hamburger/chap7/sec3.html
11
6.6 割算値の統計モデリングはやめよう
• ２項分布とロジットリンク関数を使ったロジスティック回帰を使う利点
• 「種子が生存している確率」等といった何かの生起確率を推定する
「(観測データ)/(観測データ)」といった割算値を作り出す必要がなくなる
ためである
• 割算値による回帰や相関解析は現在でも見かける手法である
この例題のような”種子の生存確率が何に依存しているのかを知りたい”と
いった時に、観測値どうしの割算をしてしまいがちである
12
• 観測値をこねくりまわして指標を創作してしまう別の例として、
観測値の変数変換もありがちな作法である
観測値の対数をとったり、あるいは複数の観測値をひとつの平均値になお
してしまうような操作
• 変形した観測値を統計モデルの応答変数にするのは不必要である
ばかりではなく、場合によっては間違った結果を導きかねないもので
ある
13
観測データどうしの割算によって生じる問題
1. 情報が失われる
• 例えば野球の打率で、”1000打数300安打”の打者と”10打数3安打”の打者
は、どちらも同じ程度に確からしい「３割打者」と言えるのだろうか？
 このように打数・安打数という２つのデータをひとつにまとめることで、確からしさの情報
が失われる
2. 変換された値の分布はどうなる？
• 分子・分母にそれぞれ誤差が入った数量どうしを割算して作られた割算値
はどんな確率分布に従うのでしょうか？
• カウントデータに１を加えて対数変換すれば、それは正規分布になるので
しょうか？
 確率分布が導出出来る場合もあるが、一般的には難しい問題である。
分子と分母が独立ではない場合、さらに複雑になる。
14
6.6 割算値いらずのオフセット項わざ
• ロジスティック回帰の統計モデルのように、「N個のうちy個で事象が
生じる確率」を明示的に扱う二項分布を使うことで、割算値の使用は
回避できる
• ここでは人口密度のように「単位面積あたりの個体数」といった概念
に興味があって、観測データを観測データで割ったような数量なん
かが使われがちな場面で、どのように統計モデルを作ればよいのか
考えてみよう
• このような場面で便利なのがオフセット項わざである
詳細URL:
http://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/49477/5/kubostat2008d.pdf (14~)
http://hosho.ees.hokudai.ac.jp/~kubo/ce/2009/kubo2009glm.pdf (38~)
15
16
• ポアソン回帰におけるオフセット項わざを、図6.10(A)のような架空
データを使う。この調査の概要は以下のようになっている。
1.
2.
3.
4.
5.
森林のあちこちに調査地100箇所を設定した(i∈ {1,2, … , 100})
調査地”i”毎にその面積𝐴𝑖 が異なる
調査地”i”の「明るさ」𝑥𝑖 を測っている
調査地”i”における植物個体数𝑦𝑖 を記録した
(解析の目的)調査地”i”における植物個体の「人口密度」が「明るさ」𝑥𝑖 に
どう影響されているかを知りたい
※人口密度とは、個体数／面積といった単位面積あたりの個体数といった概
念であるが、観測値𝑦𝑖 と𝐴𝑖 で割算値をこしらえる必要はない
⇒GLMのオフセット項を利用して解決できる
17
• 面積が𝐴𝑖 である調査地i における人口密度
平均個体数λ𝑖
𝐴𝑖
= 人口密度
• 人口密度は正の量なので、指数関数と明るさ𝑥𝑖 依存性を組合せて
λ𝑖 = 𝐴𝑖 ×人口密度 = 𝐴𝑖 exp(𝛽1 + 𝛽2 𝑥𝑖 )
• 上の式のようにモデル化し、変形する
λ𝑖 = exp(𝛽1 + 𝛽2 𝑥𝑖 + log 𝐴𝑖 )
上の式のようになり、𝑧𝑖 = exp(𝛽1 + 𝛽2 𝑥𝑖 + log 𝐴𝑖 )を線形予測子とする対数リンク関
数・ポアソン分布のGLMになる
⇒線形予測子の中でパラメーターがつかないlog 𝐴𝑖 のような項を”オフセット項”と呼ぶ
⇒線形予測子にlog 𝐴𝑖 という「げた」をはかせている、と考えればよい
18
• Rのglm()ではオフセット項を以下のように指定できる
glm(y~x,offset=log(A),family=poisson,data=d)
Coefficients:
(Intercept)
x
0.9731
1.0383
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 261.5
Residual Deviance: 81.61
AIC: 650.3
※ここではデータフレイムdにデータが格納されていて、y,x,A列がそれぞれ
調査地内での個体数・明るさ・面積を表している
※d <- read.csv("data4b.csv")が必要
※” でエラーが出た場合は自分で ” を入力し直す
19
• オフセット項を使うと、個体数平均は調査地の面積𝐴𝑖 に比例する、と
いった仮定を反映させつつ明るさ𝑥𝑖 の効果を推定可能
• 個体数を調査地面積で割って密度にする、といった観測値どうしの
割算はまったく不要
• 単位面積・単位時間当たりのカウントデータだけでなく、概念として
は[(連続値)/(連続値)]となるような比率・密度なども、オフセット項を
使った統計モデリングが可能
20
6.7 正規分布とその尤度
• 正規分布は連続値のデータを統計モデルで扱うための確率分布
であり、ガウス分布とも呼ばれる
• ポアソン分布と同じく平均値のパラメーター”μ”をもち、これは±∞の
範囲で自由に変化できる
• ポアソン分布とは異なり、値のばらつきをあらわす標準偏差を表現
するパラメーター”σ”でデータのばらつきも指定できる
2
標準偏差 = 分散
詳細URL(p.15~):
http://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/49477/3/kubostat2008b.pdf
21
• 平均パラメーターμと標準偏差パラメーターσの正規分布の数式表現
𝑝 𝑦 μ, σ =
(𝑦 − μ)2
exp{−
}
2
2
2σ
2πσ
1
• これは正規分布の確率密度関数である
• Rで正規分布の密度関数を図示(図6.11(A))
y <- seq(-5,5,0.1)
plot(y,dnorm(y,mean=0,sd=1),type ="l")
※d <- read.csv("data4b.csv")が必要
※” でエラーが出た場合は自分で ” を入力し直す
22
• 連続値の確率分布では、確率密度を積分した量が確率と定義され
ている
• 図6.11(A)-(C)では、確率変数yが1.2≤y ≤ 1.8となる確率の大小が、
グレイの領域の面積としてあらわされている
• 正規分布の確率密度関数をp(y|μ,σ)とすると、
この確率p(1.2≤y ≤ 1.8|μ,σ)は “
1.8
p(y|μ,σ)dy”
1.2
と書ける
23
• Rでこの確率を計算したい時はpnorm(x,mu,sd)関数を使用する
• 例えば、μ=0,σ=1の場合にp(1.2≤y ≤ 1.8|μ,σ)を評価(※標準正規分布)
pnorm(1.8,0,1) - pnorm(1.2,0,1)
[1] 0.07913935
確率が0.079ぐらいとわかる
• 図6.11のグレイ領域の面積計算は、y=1.2と1.8の中間の値である
y=1.5における確率密度p(y=1.5|0,1)を「高さ」、1.8-1.2=0.6を幅Δyと
する長方形であると近似してみる
dnorm(1.5,0,1) * 0.6
[1] 0.07771056
※幅Δyが小さいほど良い近似になる
確率が0.078ぐらいとわかる
24
25
正規分布の最尤推定
• N人からなる人間の集団の身長データをY={𝑦𝑖 }、個体iの身長が𝑦𝑖
• ある𝑦𝑖 が「𝑦𝑖 -0.5Δy ≤y ≤ 𝑦𝑖 +0.5Δy」である確率は確率密度関数と
区間幅Δyの積であると近似できるので、正規分布を使った統計
モデルの尤度関数は以下のようになる
L(μ, σ) =
=
𝑖𝑝
𝑦𝑖 μ, σ Δ𝑦
(𝑦𝑖 −μ)2
𝑖 2πσ2 exp{− 2σ2 }Δ𝑦
1
26
• 対数尤度関数
log L(μ, σ) = −0.5𝑁𝑙𝑜𝑔
2πσ2
−
1
2σ2
2
(
𝑦
−
μ)
+ Nlog(Δy)
𝑖 𝑖
• 区間幅Δyは定数だから、パラメーター{μ,σ}の最尤推定値に影響を与
えないため、尤度関数や対数尤度関数の表記では、 Δyやlog(Δy)を
無視して省略表記できる
log L(μ, σ) = −0.5𝑁𝑙𝑜𝑔
2πσ2
−
1
2σ2
2
(
𝑦
−
μ)
𝑖 𝑖
この式を使って
最尤推定を行う
27
• 尤度が確率密度の積である場合には、対数尤度は負の値になると
は限らない
• 連続値の確率分布を使った統計モデルでは対数尤度が正の値になったり、
AICや逸脱度が負の値になることもある
• 標準偏差が一定である正規分布のパラメーターの最尤推定が、
最小二乗法による推定と等しくなる
• 正規分布を使った統計モデルのあてはめとしてよく使われている
直線回帰については既に述べてある(3.7)
• 正規分布を部品とするGLMであり、数量的な説明変数である𝑥𝑖 を使って、
線形予測子を𝑧𝑖 = 𝛽1 +𝛽2 𝑥𝑖 、恒等リンク関数を使って平均を𝜇𝑖 = 𝑧𝑖 と指定
• GLMの最尤推定法によるパラメータ推定と、いわゆる「最小二乗法による直
線のあてはめ」は同等とみなすことができる
28
6.8 ガンマ分布のGLM
• ガンマ分布は確率変数のとりうる範囲が0以上の連続確率分布
• 確率密度関数は以下のように定義される
𝑟 𝑠 𝑠−1
𝑝 𝑦 𝑠, 𝑟 =
𝑦
exp(−𝑟𝑦)
Γ(s)
s:shapeパラメータ,r:rateパラメータ, Γ(s):ガンマ関数
• 平均は𝑠 𝑟, 分散は𝑠 𝑟 2 となり、分散 = 平均 𝑟 となる
• s=1のときは”指数分布”になる
• Rのdgamma(y,shape,rate)関数を使うとガンマ分布の確率密度を
評価可能
29
30
例題(図6.13)
• 架空植物50個体の葉の重量と花の重量の関係を調べる
• 個体iの葉重量を𝑥𝑖 、花重量を𝑦𝑖 とするが、𝑥𝑖 が大きくなるにつれ𝑦𝑖 も
大きくなっている
• 応答変数を𝑦𝑖 、説明変数を𝑥𝑖 として両者の関係をGLMで記述する
• ある個体の花の重量𝑦𝑖 が平均𝜇𝑖 のガンマ分布にしたがっている
• 応答変数𝑦𝑖 は連続値だが、重量なので正の値しかとらないため、
ばらつきは正規分布ではなく、ガンマ分布で説明する
31
32
• 平均花重量𝜇𝑖 が葉重量𝑥𝑖 の単調増加関数であり、さらに何らかの
生物学的な理由がある場合、以下のように仮定する
𝜇𝑖 = A𝑥𝑖𝑏
• この右辺でA=exp(a)とおいてから、全体を指数関数でまとめる
𝜇𝑖 = exp(a)𝑥𝑖𝑏 = exp(a+blog𝑥𝑖 )
• この両辺の対数をとると
𝑙𝑜𝑔𝜇𝑖 = a+blog𝑥𝑖
• 線形予測子a+blog𝑥𝑖 と対数リンク関数を使って平均𝜇𝑖 が与えられた
• この線形予測子では説明変数は𝑥𝑖 ではなくlog𝑥𝑖 となり、推定すべき
パラメーターは切片a と傾き bである
33
• 線形予測子を設定したGLMのパラメーターを、Rのglm()関数で推定
• glm()による推定では、平均𝜇𝑖 を決める線形予測子とリンク関数だけ
を指定すれば良いという便利さがある
• 平均・分散をshape/rateパラメーターとどう対応づけるかといったこと
を気にする必要はない
• 他のGLMの推定値と同じように指定する
• 推定された結果から得られた予測を図6.13に示す
• ここでは花重量の予測だけでなく、ガンマ分布を使って評価された
50%と90%区間の予測も示している
34
※サイトのd.Rdataが必要
glm(y ~ log(x),family = Gamma(link = "log"),data=d)
Coefficients:
(Intercept)
log(x)
-1.0403
0.6833
Degrees of Freedom: 49 Total (i.e. Null); 48
Residual
Null Deviance: 35.37
Residual Deviance: 17.25
AIC: -110.9
35
6.9 6章のまとめ
ーGLMの様々な拡張と応用ー
• GLMでは応答変数のばらつきを表現する確率分布は正規分布だけ
でなく、ポアソン分布・二項分布・ガンマ分布などが選択可能(6.1)
• 「N個の観察対象のうちk個で反応が見られた」というタイプのデータ
に見られるばらつきを表すために二項分布が使える(6.3)
• 生起確率と線形予測子を結びつけるロジットリンク関数を使ったGLM
のあてはめは、ロジスティック回帰と呼ばれる(6.4)
• 線形予測子の構成要素として、複数の説明変数の積の効果を交互
作用項が使える(6.5)
36
• データ解析でしばしば見られる観測値どうしの割算値作成や、応答
変数の変数変換の問題点を上げ、ロジスティック回帰やオフセット項
の工夫をすれば、情報消失の原因となる「データ加工」は不要(6.6)
• 連続値の確率変数のばらつきを表現する確率分布として正規分布・
ガンマ分布などがあり、これらを統計モデルの部品として使うときに
は、離散値の確率分布との違いに注意しなければならない(6.7&6.8)
37

Download Report