Document

6.GLMの応用範囲を広げる
須崎研究室 村上
目次

6.1
さまざまな種類のデータで応用できるGLM

6.2
例題

6.3
二項分布で表現する「ある・なし」カウントデータ

6.4
ロジスティック回帰とロジットリンク関数
上限のあるカウントデータ
 6.4.1
ロジットリンク関数
 6.4.2
パラメーター推定
 6.4.3
ロジットリンク関数の意味・解釈
 6.4.4
ロジスティック回帰のモデル選択
6.1さまざまな種類のデータで応用できるGLM

第2~5章では、ポアソン分布・対数リンク関数を中心に扱っています。

GLMの特徴は、確率分布・リンク関数・線形予測子の組み合わせを
指定することで様々なタイプのデータを表現できます。

下記の表は、GLMで確率分布ごとによく使われるリンク関数の一覧
です。
確率分布
(離散)
(連続)
乱数生成
よく使う
リンク関数
glm()の
family指定
二項分布
rbinom()
binomial
logit
ポアソン分布
rpois()
poisson
log
負の二項分布
rnbinom()
(glm.nb()関数)
log
ガンマ分布
rganma()
ganma
log
正規分布
rnorm()
gaussian
identity
6.2例題 上限のあるカウントデータ

二項分布は上限のあるカウントデータ(応答変数が𝒴 ∈ 0,1,2 ∙∙∙, 𝑁 の
ような場合)の範囲をとる現象のばらつきを表わすために使います。

この例題では、ある個体の生存確率𝑞𝑖 は体サイズ𝑥𝑖 や施肥処理𝑓𝑖 といっ
た説明変数によってどう変化するかという点について調べます。
6.2例題 上限のあるカウントデータ

例題データ

体サイズ𝑥𝑖

植物個体数100個
 このうち50個体を無処理(𝑓𝑖
= 𝐶)、残り50個体を施肥処理(𝑓𝑖 = T)とする。

観察種子数8個

生存種子数(応答変数)は𝒴𝑖 ∈ 0,1,2 ∙∙∙, 8 の範囲

個体iの種子の生存確率𝑞𝑖
summary(d)
6.2例題 上限のあるカウントデータ
下記の例題データ表より
・体サイズ𝑥𝑖 が大きくなると生存種子数𝒴𝑖 が多くなる。
・施肥処理を行う(𝑓𝑖 = T)と生存種子数𝒴𝑖 が多くなる。
また、例題ではどの個体でも観察種子数𝑁𝑖 が8個である為、生存種子数𝒴𝑖
の増加は生存確率𝑞𝑖 の増加と同じです。
6.3二項分布で表現する「ある・なし」カウントデータ

二項分布は、例題の種子データのように「N個のうち𝒴個が存在してい
た」といった構造のカウントデータを表現するときによく使います。

二項分布の確率分布
𝑁 𝑦
𝑝 𝑦 𝑁, 𝑞 =
𝑞 1−𝑞
𝑦

𝑁−𝑦
この例題では


𝑝 𝑦 𝑁, 𝑞 はある個体での生存確率です。
𝑁
𝑦
はNこの観察種子の中から𝒴個の生存種子を選び出す場合の数です。
6.3二項分布で表現する「ある・なし」カウントデータ
上限N=8で生起確率𝑞が 0.1,0.3,0.8 である二項分布の確率分布
6.4.1ロジットリンク関数

ロジスティック回帰では

確率分布は二項分布

リンク関数はロジットリンク関数を指定します。

二項分布では事象が生起する確率をパラメーターとして指定する必要
があり、この例題では種子の生存確率𝑞𝑖 です。

ロジットリンク関数はパラメーター𝑞𝑖 のような制約と線形予測子をうまく関
連付けるリンク関数です。
6.4.1ロジットリンク関数

ロジスティック関数について
 関数形は𝑞𝑖


= 𝑙𝑜𝑔𝑖𝑠𝑡𝑖𝑐 𝑧𝑖 =
1
1+𝑒𝑥𝑝 −𝑧𝑖
変数𝑧𝑖 は線形予測子𝑧𝑖 =𝛽1 + 𝛽2 𝑥𝑖 + ⋯
ロジスティック曲線について

logistic(z)は確率qを表わす

𝑧が小さくなると0に近づく

𝑧が大きくなると1に近づく

𝑧 = 0の時、𝑞 = 0.5

生存確率𝑞𝑖 が𝑧𝑖 のロジスティック関数と仮定すると
線形予測子 𝑧𝑖 がどんな値でも 0 ≤ 𝑞𝑖 ≤ 1 になります。
logistic <- function(z) 1 /(1+exp(-z))
z <- seq(-6,6,0.1)
plot(z,logistic(z),type = "l")
6.4.1ロジットリンク関数

ロジスティック関数と線形予測子の関係を調べる

生存確率𝑞𝑖 が体サイズ𝑥𝑖 だけに依存していると仮定すると
線形予測子は𝑧𝑖 =𝛽1 + 𝛽2 𝑥𝑖 となります。

このときに𝑞𝑖 と𝑥𝑖 の関係がパラメーター𝛽1 と𝛽2 に依存している様子を
下に示しています。
6.4.1ロジットリンク関数

ロジスティック関数と線形予測子の関係を調べる

前述のロジスティック関数を変形すると
𝑞𝑖
𝑙𝑜𝑔
= 𝑧𝑖
1 − 𝑞𝑖

この左辺のことをロジット関数といいます。
𝑞𝑖
𝑙𝑜𝑔𝑖𝑡 𝑞𝑖 = 𝑙𝑜𝑔
1 − 𝑞𝑖

ロジット関数はロジスティック関数の逆関数です。

また、ロジスティック関数の逆関数がロジット関数です。
6.4.2パラメーター推定

統計モデルをデータに当てはめてパラメーターを推定する
 尤度関数から
𝐿 𝛽𝑖
=
𝑖

𝑁𝑖 −𝑦𝑖
対数尤度関数が得られます。
𝑙𝑜𝑔𝐿 𝛽𝑖 =

𝑁𝑖 𝑦
𝑞 𝑖 1 − 𝑞𝑖
𝑦𝑖
𝑖 𝑙𝑜𝑔
𝑁𝑖
𝑦𝑖
+ 𝑦𝑖 𝑙𝑜𝑔 𝑞𝑖 + 𝑁𝑖 − 𝑦𝑖 𝑙𝑜𝑔 1 − 𝑞𝑖
このlogLを最大にする推定値のセット 𝛽𝑗 を探し出すのが最尤推定です。
6.4.2パラメーター推定

Rによる推定計算

応答変数cbind(y,N-y)はcbind(生存数,死亡数)を表わします。
glm(cbind(y,N-y)~x+f,data=d,family=binomial)

この結果からおおよそ 𝛽1 , 𝛽2 , 𝛽3 = −19.5,1.95,2.02 となっています。
これらの推定値を使うと下に示しているような予測の曲線が描けます。
6.4.3ロジットリンク関数の意味・解釈

ロジスティック関数の逆関数であるロジット関数は
𝑞𝑖
𝑙𝑜𝑔𝑖𝑡 𝑞𝑖 = 𝑙𝑜𝑔
1 − 𝑞𝑖

これは線形予測子に等しいので
𝑞𝑖
1−𝑞𝑖
= exp (線形予測子)
= exp(𝛽1 + 𝛽2 𝑥𝑖 + 𝛽3 𝑓𝑖 )
= exp(𝛽1 ) exp(𝛽2 𝑥𝑖 ) exp(𝛽3 𝑓𝑖 )

左辺の𝑙𝑜𝑔
𝑞𝑖
はオッズと呼ばれる量で
1−𝑞𝑖
この場合だと(生存する確率)/(生存しない確率)の比です。
6.4.3ロジットリンク関数の意味・解釈

また、ロジスティック回帰のモデルでは
このオッズはexp(パラメーター×要因数)に比例しています。

つまり、推定値 𝛽2 , 𝛽3 を代入して定数である𝑒𝑥𝑝 𝛽1 を省略すると
𝑞𝑖
∝ exp(1.95𝑥𝑖 )exp 2.02𝑓𝑖
1 − 𝑞𝑖
6.4.3ロジットリンク関数の意味・解釈

生存確率のオッズに対する植物個体の体サイズ𝑥𝑖 の影響を調べます

個体𝑖が「1単位」増加した場合
𝑞𝑖
∝ exp(1.95 𝑥𝑖 + 1 )exp 2.02𝑓𝑖
1 − 𝑞𝑖
∝ exp(1.95𝑥𝑖 )exp(1.95)exp 2.02𝑓𝑖

上記のように変化し、生存のオッズはexp(1.95) = 約7倍増加します。

同様に「肥料なし」(𝑓𝑖 = 0)に比べ、
「肥料あり」 (𝑓𝑖 = 1)のオッズはexp(2.02) = 約7.5倍増加します。

このように、ロジットリンク関数で生存確率を定義することによって、
様々な要因と応答事象のオッズの解釈が簡単になります。
6.4.3ロジットリンク関数の意味・解釈

また、前述のオッズの対数をとると
𝑞𝑖
= 𝛽1 + 𝛽2 𝑥𝑖 + 𝛽3 𝑓𝑖
1 − 𝑞𝑖


となります。右辺の𝛽1 + 𝛽2 𝑥𝑖 + 𝛽3 𝑓𝑖 は線形予測子𝑧𝑖 そのものです。
オッズと「リスク」について考える

ここでの「リスク」とは(近似的には)オッズ比のことです。
6.4.3ロジットリンク関数の意味・解釈

ある人間集団におけるナントカ病の発病と生活習慣に基づいた
データ解析をしたとします。

個人𝑖の生活習慣𝑋の効果を表わす係数が𝛽𝑠 であるとする。

発病確率をロジスティック回帰で調べ、その最尤推定値𝛽𝑠 = 1.95が
得られたとする。

この場合
Xの𝑜𝑑𝑑𝑠
𝑒𝑥𝑝 X・非Xの共通部分 × 𝑒𝑥𝑝 1.95 × 1
=
非Xの𝑜𝑑𝑠𝑠 𝑒𝑥𝑝 X・非Xの共通部分 × 𝑒𝑥𝑝 1.95 × 0
= 𝑒𝑥𝑝1.95

病気になるオッズ比(リスク)は exp(1.95) = 約 7倍と見積もられます。
6.4.4ロジスティック回帰のモデル選択

第4章と同じくAICによるモデル選択を行う

6.4.2項の体サイズ𝑥𝑖 と施肥処理𝑓𝑖 を説明変数とするモデルが最も
良い予測をするモデルなのか不明です。

ネストしてるモデルたち、すなわち、説明変数をどちらか一方だ
け使ったモデル、あるいはどちらも使わないモデルのほうが良い
予測を得られるかもしれません。
6.4.4ロジスティック回帰のモデル選択

第4章と同じくAICによるモデル選択を行う

RのstepAIC()関数を使うとネストしているモデルのAICを比較し、
最少AICモデルを選択できます。
fit.xf <- glm(cbind(y, N - y) ~ x + f, data = d,family = binomial)
library(MASS)
stepAIC(fit.xf)
6.4.4ロジスティック回帰のモデル選択

前述の比較結果のまとめ

f,x,x+fモデルはそれぞれ、施肥処理𝑓𝑖 のみ、体サイズ𝑥𝑖 のみ、
両者に依存するモデルを表わしています。

この表から、体サイズと施肥処理を同時に組み込んだx+fモデルが
AICの観点から最良になりました。