一般化線型モデル generalized linear model; GLM

Presentation on theme: "一般化線型モデル generalized linear model; GLM"— Presentation transcript:

1 一般化線型モデル generalized linear model; GLM

2 一般化線型モデルとは 一般線型モデル（general linear model; GLM）という別ものがある
重回帰、分散分析、共分散分析、ロジスティック回帰、ポアソン回帰などが包含される統計モデル 応答変数がしたがう確率分布：正規分布、ガンマ分布、ポアソン分布、2項分布、逆ガウス分布（これらは、指数型分布族といわれる） リンク関数：説明変数の線形結合と応答変数の期待値の関係 モデルのあてはめ（パラメータの推定）：最尤法 実際には、反復再重み付け最小二乗法iteratively reweighted least squares (IRWLS; 最尤法と同じ結果になる）が使われることが多い 一般線型モデル（general linear model; GLM）という別ものがある 重回帰、分散分析、共分散分析などが包含される統計モデル 応答変数がしたがう確率分布：正規分布 説明変数の線形結合がそのまま応答変数の期待値となる モデルのあてはめ（パラメータの推定）：最小自乗法

3 確率関数 離散変数の場合：確率関数 この関数で確率が計算できる（下グラフの縦軸の値が確率になっている） 例： 2項分布、ポアソン分布など
2つのパラメータ（n, p）を持つ n: ベルヌーイ試行の回数、p: 成功確率 取り得る値：0からnまでの整数 例1：玉がたくさん入っている袋（白玉の割合p）からn個の玉を取り出すとき白玉の個数 例2：最初にn本ある樹木の内生き残る本数 ポアソン分布 パラメータλを1つもつ、λは平均かつ分散 所与の空間・時間の範囲内で事象が起こる回数 取り得る値：0から∞までの整数 例1：店の中の客の数 例1：1年間にプロットに発生する実生の数

4 確率密度関数 連続変数の場合：確率密度関数 下グラフの縦軸の値は確率ではない。 例：正規分布、ガンマ分布など 正規分布
以上4つはよく使うものなので、覚えておいてください。 正規分布 中央値、分散（広がり方を決めるパラメータ）をもつ -∞<x<∞ ガンマ分布 分布の形を決める2つのパラメータをもつ 様々な形 x>0で値を持つ

5 応答変数がしたがう確率分布 何を選ぶか？考えるべきこと 連続か離散か？ マイナスの数字をとるか？ 値に上限があるか（離散値）？
連続：胸高直径、樹高 →正規分布、ガンマ分布 離散：実生の数、ネズミの数 →ポアソン分布、2項分布 マイナスの数字をとるか？ 負の数字にもなりうる、予測値が0に近づかない →正規分布 値に上限があるか（離散値）？ 上限なし →ポアソン分布 上限あり （0,1の内どれかなど） →二項分布分布

6 リンク関数 説明変数の線形結合ηと応答変数の期待値μとの関係 x1, x2, … xp: 説明変数 η: 説明変数の線形結合
μ: 応答変数の分布の平均 限定された範囲の数字（正値、非負値、[0,1]）しかとらない応答変数への対応 - ηは[-∞, ∞]の数字なので、ηがそのままμにすると問題が生じる 非直線的な反応、応答変数の変数変換がリンク関数で処理できる

7 例1 ポアソン回帰 応答変数がしたがう 確率分布：ポアソン分布 母シュート=60cmの時の確率 応答変数が正規分布にしたがうと仮定すると…
例1　ポアソン回帰 応答変数が正規分布にしたがうと仮定すると… あり得ない数（非整数、負の数）のシュートがある確率で生じることになる 応答変数がしたがう 　確率分布：ポアソン分布

8 例1 ポアソン回帰 リンク関数：log-link （娘シュート数の平均）=exp（a（母シュート長）+b）
例1　ポアソン回帰 リンク関数：log-link （娘シュート数の平均）=exp（a（母シュート長）+b） こうすることで、娘シュート数の平均が負の数字になることをふせぐ。

9 canonical link function 正準連結関数
canonical link function: 数学的・計算的に便利で、多くの場合自然な選択（canonicalでないlinkも可能） 分布族（family） リンク（link） 正規分布（Normal）） η=μ ポアソン分布（Poisson） η=ln(μ) 2項分布（Binomial） η=ln(μ/(1-μ)) ガンマ分布（Gamma） η=μ-1 逆ガウス分布（Inverse Gaussian） η=μ-2

10 モデルのあてはまりの良さ r2（決定係数）は、応答変数が正規分布に従うときにのみ有効な指標 →もっと一般的な指標が必要
log likelihood（model fitの時に最大化するもの）を基礎にする。ただし、log likelihoodはデータ数によって大きさがかわる Deviance: log likelihoodに-2をかけた数字を基礎に計算した数値 null devianceとresidual devianceの二つの数字をみる… Null Deviance Deviance 対数尤度に-2を かけたもの Residual Deviance 飽和したモデルのDeviance 指定したモデルのDeviance 定数項のみのモデルのDeviance

11 GLMにおける仮説検定 尤度比検定（likelihood ratio test）
説明変数x1, x2,…xnを含むモデルをあてはめたときの最大尤度と説明変数の一部（例えば、x1,…,xp-1, xp+1 …xn）のみを含むモデル（部分モデル）の最大尤度の比を使う検定 変数（例えば、xp）をモデルからぬくことで尤度はどれだけ減少するか？減少が大きければ、抜いた変数の効果は有意である。

12 GLMにおける仮説検定 尤度比検定（likelihood ratio test） 尤度の比＝対数尤度の差なので、対数尤度の差を使います。
帰無仮説（ぬいた変数の効果なし）のもとで、対数尤度の差にｰ2をかけたもの（尤度比検定統計量）がχ2分布（自由度：抜いた変数の個数）に近似的に従うことが分かっているので、尤度比検定統計量とχ2分布を比較し有意性を判定します。

13 GLMにおける仮説検定 個々の変数の有意性 →Wald検定 変数の推定値を変数推定値の標準誤差で割った値をχ2分布（自由度：1）と比較する。
この検定はあまり正確でないことが分かっています。

14 Rを使ったデータ処理 （データ準備+グラフ描き）
L <- c( 113, 90, 57, 65, 62, 75, 80 ) shootNumber <- c( 6, 3, 1, 1, 2, 1, 2 ) plot( L, shootNumber, xlab="Shoot Length (cm)", ylab="Offspring Shoot Number" ) グラフ描き （）の中：横軸に使う変数、縦軸に使う変数、 xlab=横軸のラベルの文字列、 ylab =縦軸のラベルの文字列）

15 Rを使ったデータ処理 （モデルの当てはめ）
poissonRegResult <- glm( formula = shootNumber ~ L, family = poisson( link = "log" ) ) #一般化線型モデルを当てはめるための関数glm()を使って、 #娘シュート数を母シュート長にポアソン回帰する。　 #formula = offspringShootNumber ~ L 　モデルを指定する #family = poisson( link = “log” ) 　応答変数がポアソン分布に従うことと #logリンク関数を使用することを指定する。 #結果は、 poissonRegResultというオブジェクトに代入する。 summary(poissonRegResult ) #poissonRegResultに入っている結果を要約して出力する。 coef( poissonRegResult ) #推定したパラメータを出力する ab

16 Summaryの中身 パラメータの推定値、標準誤差、推定値/標準誤差、Wald testの結果
Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) L * パラメータの推定値、標準誤差、推定値/標準誤差、Wald testの結果 Null deviance: on 6 degrees of freedom Residual deviance: on 5 degrees of freedom AIC: Null deviance：定数項のみのモデルと飽和したモデル（すべての応答変数が残差なく説明されるよう、応答変数の個数だけ説明変数を使ったモデル）の対数尤度の差に-2をかけたもの Residual deviance：指定したモデルと飽和したモデルの対数尤度の差に-2をかけたもの Null Deviance Deviance 対数尤度に-2を かけたもの Residual Deviance 飽和したモデルのDeviance 指定したモデルのDeviance 定数項のみのモデルのDeviance

17 Rを使ったデータ処理 （予測値のグラフへの追加）
shootLength <- seq( 55, 125, length = 15 ) #55, 60, 65 … 115という数列を作る #最小55、最大125で長さ（数字の数）15の数列を作るの意 linearCombination <- coef( poissonRegResult )[ 1 ]　+ coef( poissonRegResult )[ 2 ] * shootLength #推定したパラメータを使って、55, 60, 65 … 115に対応した線型結合を作る pred <- exp( linearCombination ) #線型結合の指数（log-linkの逆数）を計算し、娘シュート数の予測値を計算する lines( shootLength, pred ) #shootLengthを横軸の値、 pred を縦軸の値とした折れ線をグラフに追加する

18 Rを使ったデータ処理 （尤度比検定） reducedModel <- update( poissonRegResult, ~.- L ) #poissonRegResultから変数Lを除いたモデルの当てはめを行い、 #結果をreducedModel に入れる。 summary( reducedModel ) #reducedModel の要約の出力 #null devianceとresidual devianceが等しくなっている anova( reducedModel, poissonRegResult, test = "Chi" ) #poissonRegResult、 reducedModelの対数尤度の差→尤度比検定 Model 1: offspringShootNumber ~ 1 Model 2: offspringShootNumber ~ L Resid. Df Resid. Dev Df Deviance P(>|Chi|) 2つのモデルに含まれている変数の数 2つのモデルのresidual devianceと両者の差（この差が尤度比検定統計量） χ2分布と比較して出したp-value

19 例2 ロジスティック回帰 応答変数がしたがう 確率分布：2項分布 x=10の時→p=0.17 n=1の試行をして、
例2　ロジスティック回帰 x=10の時→p=0.17 n=1の試行をして、 1回成功する確率は0.17 0回成功する（1回失敗する）確率は0.83 応答変数がしたがう 　確率分布：2項分布

20 例2　ロジスティック回帰 リンク関数：logit-link こうすることで、pは[0, 1]に収まる

21 例2 ロジスティック回帰 データ準備・グラフ描き
例2　ロジスティック回帰 データ準備・グラフ描き survival <- c( rep( 0, 8 ), 1, rep( 0, 3 ), 1, 1, 0, 0, 1, 0, rep( 1, 12 ) ) size <- 1:30 plot( size, survival, xlab = "Size", ylab = "Survival", type = "p" )

22 Rを使ったデータ処理 （モデルの当てはめ）
logisticRegResult <- glm( survival ~ size, family = binomial( link = "logit" ) #一般化線型モデルを当てはめるための関数glm()を使って、 #（生存or死亡）を個体サイズにロジスティック回帰する。　 #survival ~ size　モデルを指定する #family = binomial( link = “logit” )　応答変数が二項分布に従うことと #logitリンク関数を使用することを指定する。 #結果は、 logisticRegResultというオブジェクトに代入する。 summary( logisticRegResult ) #logisticRegResultに入っている結果を要約して出力する。 coef( logisticRegResult ) #推定したパラメータを出力する ab

23 Summaryの中身 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ** size ** --- Signif. codes: 0 ‘***’ ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: on 29 degrees of freedom Residual deviance: on 28 degrees of freedom AIC: Null deviance：定数項のみのモデルと飽和したモデル（すべての応答変数が残差なく説明されるよう、応答変数の個数だけ説明変数を使ったモデル）の対数尤度の差に-2をかけたもの Residual deviance：指定したモデルと飽和したモデルの対数尤度の差に-2をかけたもの Null Deviance Deviance 対数尤度に-2を かけたもの Residual Deviance 飽和したモデルのDeviance 指定したモデルのDeviance 定数項のみのモデルのDeviance

24 Rを使ったデータ処理 （予測値のグラフへの追加）
xForSize <- seq( from = 1, to = 30, length = 60 ) #1… 30という数列を作る #最小1、最大30で長さ（数字の数）60の数列を作るの意 linearCombination <- coef( logisticRegResult )[[ 1 ]]+ coef( logisticRegResult )[[ 2 ]] * xForSize #推定したパラメータを使って、1… 30に対応した線型結合を作る pred <- exp( linearCombination ) / ( exp( linearCombination ) + 1 ) #logistic関数（log-linkの逆数）に線型結合を代入し、生存率の予測値を計算する lines( xForSize, pred ) #xForSizeを横軸の値、 #pred を縦軸の値とした折れ線をグラフに追加する

25 Rを使ったデータ処理 （尤度比検定） reducedModel <- update( logisticRegResult , ~.- size ) #logisticRegResultから変数sizeを除いたモデルの当てはめを行い、 #結果をreducedModel に入れる。 summary( reducedModel ) #reducedModel の要約の出力 #null devianceとresidual devianceが等しくなっている anova( reducedModel, logisticRegResult, test = "Chi" ) #logisticRegResult、 reducedModelの対数尤度の差→尤度比検定 Model 1: survival ~ 1 Model 2: survival ~ size Resid. Df Resid. Dev Df Deviance Pr(>Chi) e-06 *** 2つのモデルに含まれている変数の数 2つのモデルのresidual devianceと両者の差（この差が尤度比検定統計量） χ2分布と比較して出したp-value

26 まとめ GLM: 重回帰、ロジスティック回帰、ポアソン回帰などを包括する一般的統計モデル
GLM: link-functionを仮定する GLMのあてはめ、パラメータ推定: 最尤法 GLMにおける仮説検定: 尤度非検定 参考図書 「データ解析のための統計モデリング入門」　久保拓哉　著 岩波書店