Transcript glm2

一般化線形モデル（GLM）
generalized linear Models
推定結果の確認と検定
Rで学ぶデータサイエンス10
一般化線形モデル
粕谷英一（２０１２） 共立出版
Generalized Linear Models
• Linear Model
 yi  1   2 xi  3 f i  
response variable ～ intercept + slope * explanatory variable
 lm(y~ x + f ・・・)，lm(y~x + f -1) (no intercept)
Generalized Linear Model
 f ( yi )  1   2 xi  3 f i  
Model &Link function ～ intercept + slope * explanatory variable
 glm(y ~ x, data = d, family = poisson)
尤度（ｐ12）
• ある確率分布でパラメータの値θが決まれば，デー
タXの値ｘについてその値が得られる確率（確率密度
）が計算できる．
– f(x|θ）
– R上では ｄ確率分布名（ｘ，θ）の形
#一様分布 (unif)の例
#確率密度関数のグラフ
curve(dunif(x,min=0,max=2),xlim=c(-0.5,3),ylim=c(0,1),xlab="y",ylab="probability density")
#ある値に対する確率密度の値はdunif関数
dunif(0.2, min=0,max=2.0)
#分布関数，累積分布関数：変数がある値以下を取る確率：puniｆ関数
curve(punif(x,min=0,max=2),xlim=c(-0.5,3),ylim=c(0,1),xlab="y",ylab="probability")
#分位数(quantile）その値以下を取る確率がｐであるような点の値，分布関数の逆関数
qunif(0.75,min=0,max=2.0)
#乱数の発生：runif関数 ，乱数の個数とパラメータを与える
runif(3,min=0, max=2.0)
尤度（ｐ12）
• ある確率分布でパラメータの値θが決まれば，データX
の値ｘについてその値が得られる確率（確率密度）が
計算できる．
– f(x|θ）
– R上では ｄ確率分布名（ｘ，θ）の形
• 逆に，データX=xが与えられたとき，パラメータ
の値θに対して，その値ｘが得られる確率を尤
度：ゆうど(likelihood）という．
二項分布の例と尤度関数
• つぼのなかに赤球ｒ個，白球ｗ個あり，1つ取り
出して色を記録して戻すことをｎ回繰り返す。
• 赤が出る回数Yがｙを取る確率は，一つの母数
φ＝r/(r+w)を用いると，
n-y
y
P(Y = y | f ) = n Cyf (1- f ) となる．
• 実際に赤が8回，白が2回でた場合には，その
2
8
ことが起こる確率は， 10 C8 f (1- f )
で，これを母数φの関数と見なしたものを尤度
関数L(φ)と呼ぶ．
二項分布の例と尤度関数
#二項分布の関数形：Rではdbinom
barplot(dbinom(0:10,size=10,prob=0.6),ylab="probability
",space=0, names=as.character(0:10), col="white")
#赤が8回，白が2回でた場合の尤度関数L(φ)
Lik <- function(phi) {dbinom(8,size=10,phi)}
curve(Lik(x), 0, 1)
#尤度関数の対数値を対数尤度関数(LogLikelihood)
LLik <- function(phi) {log(dbinom(8,size=10,phi))}
curve(LLik(x), 0.05, 0.95)
尤度の最大化（最尤推定）
• データがあり，確率分布の種類は決まっているが，パ
ラメータ（母数）値がわからないとき。
0.15
0.10
0.05
0.00
Lik <- function(phi) {dbinom(8,size=10,phi)}
optimize(Lik,c(0,1),maximum=TRUE)
Lik(x)
= 10 C8 f 7 (1- f ) (8 -10f )
最大値はφ=8/10＝0.8で取る．
0.20
0.25
0.30
• 得られているデータがもたらされる確率（尤度）が高い
パラメータ値だったと考えるのが自然．
• 尤度が最大になるパラメータ値を推定値として使う．
• 赤が8回，白が2回でた場合の尤度関数
2
8
C
f
1f
(
)
10
8
これを母数φで微分すると，
7
C
f
(1- f ) {8(1- f ) - 2f }
10 8
0.0
0.2
0.4
0.6
x
0.8
1.0
尤度の最大化（最尤推定）
• 赤が8回，白が2回でた場合の尤度関数，10 C8 f (1- f )
対数尤度関数は，
8
log( 10 C8 )+8log f + 2log (1- f )
これを母数φで微分すると，
8
2
8(1- f ) - 2f 8 -10f
=
=
f (1- f )
f (1- f )
f (1- f )
-10
-15
LLik <- function(phi) {log(dbinom(8,size=10,phi))}
optimize(LLik,c(0.01,0.99),maximum=TRUE)
-20
LLik(x)
-5
最大値は最後の分子が0になる，
φ=8/10＝0.8で取る．
0.2
0.4
0.6
0.8
2
スコア関数（尤度の微分）
• 尤度関数をパラメータで偏微分したものをス
コア関数と呼ぶ
• 最尤推定値は，スコア関数＝0の解．
• パラメータが複数あるときは，各パラメータに
対するスコア関数＝0の連立方程式を解く．
0.00
-0.02
-0.04
Scor(x)
#スコア関数の定義
Scor <- function(phi){phi^7*(1-phi)*(8-10*phi)}
#スコア関数のグラフ
curve(Scor(x),0,1)
abline(h=0.0)
#スコア関数の＝0の数値解
uniroot(Scor, c(0.05,0.95))
0.02
7
C
f
(1- f ) (8-10f )
10 8
0.0
0.2
0.4
0.6
0.8
1.0
3つの検定方法（ｐ16）
• 2つのモデルを比較する
– 帰無仮説のモデル：パラメータが0（θo）
– 対立仮説のモデル：パラメータは最尤推定値
• Wald検定：帰無仮説モデルが正しければ，最尤推
定量はθoを中心とする正規分布に従うことを用いて
，最尤推定値が得られる確率を計算する
• スコア検定：帰無仮説のモデルが正しければ，ス
コアの絶対値が大きくなる可能性が小さいことを利
用
• 尤度比検定：帰無仮説のモデルと対立仮説のモ
デルの対数尤度の差が，カイ2乗分布に従うことを
利用
3つの検定方法（ｐ16）
• Wald検定：帰無仮説モデルが正しければ，最尤推定量は
θoを中心とする正規分布に従うことを利用
• スコア検定：帰無仮説のモデルが正しければ，スコアの絶
対値が大きくなる可能性が小さいことを利用
• 尤度比検定：帰無仮説のモデルと対立仮説のモデルの対
数尤度の差が，カイ2乗分布に従うことを利用
スコア検定
対数尤度
θoでの接線の傾き
尤度比検定
対数尤度の差
デビアンス
尤離度
2つのモデルの対数
尤度の差の2倍
Wald検定
θの最尤推定値の離れ
帰無仮説
モデルでの値θo
最尤推定値
θML
説明変数が目的変
数に対して「偶然」を
超える効果を与えて
いるかを検討する．
その説明変数がない
モデル（帰無モデル）
でも十分起こりうる結
果か？を確認．
パラメータθ
予測の最適化とAIC（ｐ18）
• 得られているデータに対する尤度について，説明変
数を追加すると尤度は大きくなるか変わらない
– 全ての説明変数を使うと，得られているデータにはよく当
てはまるものの，それに引きずられ，次の別の説明変数
の値が新たに得られた時に目的変数の値を予測するに
は適切でない可能性がある。
• 赤池情報量基準（AIC)の小さいモデルを選ぶ
AIC=-2×（そのモデルでの最大対数尤度)+2×（パラメータ数）
– パラメータの数を多く入れすぎない
Poisson Model (p49)
（counting data of occurrence)
Poisson Model for number of seeds of a plant,
regressed on plant size and nutrification (p49)
Maximize log-likelihood y
log L( 1 ,  2 ,  3)   log
i exp( i )
i
yi !
log( i )  1   2 xi  3 f i
glm(y ~ x + f, data = d, family = poisson)
i
#page 19
x1 <- c(6.5,3.8,3.4,2.4,3.0,5.5,2.4,6.6)
x2 <- c(3.7,4.9,1.0,1.8,4.6,4.8,3.8,2.7)
y1 <- c(8, 5, 2, 0, 1, 11, 4, 9)
fit <- glm(y1~x1+x2, family=poisson)
summary(fit)
#説明変数の係数と切片
coef(fit)
#目的変数値と，線形予測子の値
predict(fit, type="response")
predict(fit, type="link")
#説明変数値が異なる値の場合の予測値
predict(fit,newdata=data_new1,type="link")
推定結果の利用（ｐ23）
#page 19
#残差
x1 <- c(6.5,3.8,3.4,2.4,3.0,5.5,2.4,6.6)
#目的変数観測値の残差
x2 <- c(3.7,4.9,1.0,1.8,4.6,4.8,3.8,2.7)
residuals(fit, type="response")
y1 <- c(8, 5, 2, 0, 1, 11, 4, 9)
#デビアンス誤差：デビアンスの平方根
fit <- glm(y1~x1+x2, family=poisson)
residuals(fit, type="deviance")
summary(fit)
#ピアソン誤差：残差/分散の平方根
#説明変数の係数と切片
residuals(fit, type="pearson")
coef(fit)
#目的変数値と，線形予測子の値
predict(fit, type=“response”)
predict(fit, type=“link”)
#説明変数値が異なる値の場合の予測値
predict(fit,newdata=data_new1,type=“link”)
検定（ｐ26）
#Wald検定:summaryで表示される
summary(fit)
#尤度比検定：anovaを呼び出す
anova(fit, test="Chisq")
anova(fit, test="LRT")
#スコア検定：anovaを呼び出す
anova(fit, test="Rao")
#AICの値
AIC(fit)
extractAIC (fit)