一般化線形モデル（GLM） generalized linear Models 推定結果の確認と検定 Rで学ぶデータサイエンス10 一般化線形モデル 粕谷英一（２０１２）　共立出版

Generalized Linear Models 　response variable ～　intercept + slope * explanatory variable 　 lm(y~ x + f ・・・)，lm(y~x + f -1) (no intercept) Generalized Linear Model 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)を用いると， 　　　　　　　　　　　　　　　　　　　　　　　となる． 実際に赤が8回，白が2回でた場合には，そのことが起こる確率は， 　で，これを母数φの関数と見なしたものを尤度関数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)

尤度の最大化（最尤推定） データがあり，確率分布の種類は決まっているが，パラメータ（母数）値がわからないとき。 得られているデータがもたらされる確率（尤度）が高いパラメータ値だったと考えるのが自然． 尤度が最大になるパラメータ値を推定値として使う． 赤が8回，白が2回でた場合の尤度関数 　これを母数φで微分すると， 最大値はφ=8/10＝0.8で取る． Lik <- function(phi) {dbinom(8,size=10,phi)} optimize(Lik,c(0,1),maximum=TRUE)

尤度の最大化（最尤推定） 赤が8回，白が2回でた場合の尤度関数， 対数尤度関数は， これを母数φで微分すると， 　対数尤度関数は， これを母数φで微分すると， 最大値は最後の分子が0になる， φ=8/10＝0.8で取る． LLik <- function(phi) {log(dbinom(8,size=10,phi))} optimize(LLik,c(0.01,0.99),maximum=TRUE)

スコア関数（尤度の微分） 尤度関数をパラメータで偏微分したものをスコア関数と呼ぶ 最尤推定値は，スコア関数＝0の解． パラメータが複数あるときは，各パラメータに対するスコア関数＝0の連立方程式を解く． #スコア関数の定義 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))

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 glm(y ~ x + f, data = d, family = poisson) #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) #残差 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”) #残差 #目的変数観測値の残差 residuals(fit, type="response") #デビアンス誤差：デビアンスの平方根 residuals(fit, type="deviance") #ピアソン誤差：残差/分散の平方根 residuals(fit, type="pearson")

検定（ｐ26） #Wald検定:summaryで表示される summary(fit) #尤度比検定：anovaを呼び出す anova(fit, test="Chisq") anova(fit, test="LRT") #スコア検定：anovaを呼び出す anova(fit, test="Rao") #AICの値 AIC(fit) extractAIC (fit)