発達したいブログ

コンサータとかいう人類の叡智の結晶

【自分用メモ】久保拓弥著『データ解析のための統計モデリング入門』第6章

本当に個人メモ

6.1 さまざまな種類のデータで応用できるGLM

表6.1 確率分布とリンク関数の対応一覧

  確率分布 乱数生成

glm()の

family指定

よく使うリンク関数
(離散) 二項分布 rbinom() binomial logit
(離散) ポアソン分布 rpois() poisson log
(離散) 負の二項分布 rbinom() (glm.nb()関数) log
(連続) ガンマ分布 rgamma() gamma logかな?
(連続) 正規分布 rnorm() gaussian identity

 

6.2 例題:上限のあるカウントデータ

・上限があるので、ポアソンではなく二項分布を使う

・例題で使われるデータの説明

  • ある植物の個体iそれぞれにおいて、「Ni個の観察種子のうち生きていて発芽能力があるのはyi個、死んだのはNi-yi個」という観測データがあるとする。
  • 固体ごとに体サイズや生育環境に左右されて生き残る確率は上下するとする
  • 全部で100個体調べたとする。(個体1~個体100)
  • どの個体も観察種子数は8とする
  • 個体iの種子の生存確率はqiとする 
  • ※生存確率 = ある個体iから得られた1個の種子が生きている確率

・説明変数の説明

  • 個体の大きさ:変数xi
  • 肥料あげた(T) / あげない(C):変数f
  • つまり生存確率qiが大きさxiとか施肥fiでどう変化するか調べる
  • ・プロットすると下図f:id:takuro52006:20170502011738p:plain

データのサマリは以下

  f:id:takuro52006:20170505134532p:plain

・なんとなくわかること

  • 体サイズxiが大きくなると生存種指数yiが多くなるっぽい
  •  肥料をやると(fi=T)生存種子数yiが多くなるっぽい

 6.3 二項分布で表現する「あり・なし」カウントデータ

・「N個のうちy個が生存していた」みたいなデータでは、二項分布がよく使われる

・パラメタである生存確率qiの値によって、分布の形がいろいろかわる

 

6.4 ロジスティック回帰とロジットリンク関数

・ロジスティック回帰は、二項分布を使ったGLMのひとつ

6.4.1 ロジットリンク関数

・パラメータとして生存確率qiを指定する必要がある

・確率なので、0qi≦1 

・こういうときロジットリンク関数が使える

・ちなみにロジスティック関すの関数形は、以下になる

  qi = logistic(zi) = 1 / {1 + exp( -zi)}

・プロットすると、下図

f:id:takuro52006:20170502014642p:plain

・これの、logistic(z)のとこがqiで、z=β1i+β2ixiになるイメージ

・また、qiを変形すると、以下になる

 zi = log{ qi / (1-q) }

・この右辺のlog { qi / (1 - q) }をロジット関数とよぶ

 logit( qi ) = log{ qi / (1-q) }

・要するに、ロジット関数とロジスティック関数は、逆関数の関係にある

6.4.2 パラメーター推定

・yiがでる確率を全部かけると尤度関数ができるので、これを対数尤度関数にすると、logL( {Bi} )となる。(数式はめんどうなので割愛)

※qi は、二項分布の公式でだせるよ。 nCr x**r (1-x)**rみたいなの。

・そこまでわかれば、以下のglm()関数がかける。

 glm(cbind(y, N - y)  ~ x + f, data = d, family = binomial )

・そうすると、 {β1, β2, β3} = {-19.536, 1.952, 2.022}となる。

 

6.4.3  ロジットリンク関数の意味・解釈

・ロジスティック関数の逆関数であるロジット関数はzi = logit(qi) = log { qi / (1-qi) } となりこれが線形予測子に等しい

※ロジスティック関数はqi = logistic(-zi) = 1 / { 1 - exp(zi) }

・qi / (1-qi) = exp(線形予測子)= exp(β1 + β2xi + β3fi) = exp(β1) + exp(β2xi) + exp(β3fi)となる。

・このqi / 1-qiは、オッズ(odds)とよばれる

・オッズはこの場合は(生存する確率) / (生存しない確率) という比になる

・個体iの大きさが1単位増加すると、exp{1.95 (xi +1) }= exp(1.95xi) * exp(1.95) とかけるので、生存確率のオッズはexp(1.95)倍、つまり約7倍くらい増加する

・同様に施肥の場合は、exp(2.02)倍なので、7.5倍くらい増える

・このようにロジットリンク関数で生存確率を定義すると、さまざまな要因と応答事象のオッズの解釈が簡単にある

 

6.4.4  ロジスティック回帰のモデル選択

・あてはめはしたが、これが種子の生存数をもっとも良く予測するモデルかは不明

・説明変数はどちらか一方とか、あるいはどちらも使わないとか、のが予測はよいかもしれない

AICによるモデル選択

モデル  k

logL

deviance

-2logL

residual

deviance

AIC

一定 1 -321.2 642.4 499.2 644.4
f 2 -316.9 633.8 490.6 637.8
x 2 -180.2 360.3 217.2 364.3
x+f 3 -133.1 266.2 123.0 272.2
フル 100 -71.6 143.2 0.0 343.2

 

 6.5  交互作用項の入った線形予測子

交互作用項を追加してみる。交互作用項とは、体サイズxiと施肥処理の効果fiの「積」の効果

・交互作用項をいれた線形予測子は、以下のようになる

 logit(qi) = β1 + β2xi + β3fi + β4xifi

・glmで計算すると、{β1,β2,β3,β4} = {-18.5233,1.8525,-0.0638,0.2163}

・施肥処理しなかった場合:logit(qi) = 18.5 + 1.85xiくらい

・施肥処理した場合:logit(qi) = -18.5 -0.0638 + (1.85 + 0.216)xi  = -18.6 + 2.07xi くらい

AICは274くらいなので、x + fモデルよりも悪化している

・交互作用項をいれるときに注意することは、むやみに交互作用項をいれないこと

 

6.6  割算値の統計モデリングはやめよう

・ロジットリンク関数を使ったロジスティック回帰を使う利点のひとつは、生起確率を推定するときに(観測データ) / (観測データ)という割り算値を作り出す必要がなくなること

・観測値どうしの割り算とか、観測値の変数変換とか、こねくり回して応答変数をつくるのと、間違った結果を導きかねない

・例えば、観測データどうしの割り算だと以下のようなことが発生

  • 情報が失われる:例えば1000打数300安打の打者と、10打数3安打の打者をどちらも同じ「3割打者」として扱うと、「確からしさ」の情報が失われる
  • 変換された値の分布:分子、分母にそれぞれ誤差が入った数どうしを割算して作られた割算の値はどんな確率分布に従う?カウントデータに1を加えて対数変換すれば正規分布になる???
6.6.1 割算値いらずのオフセット項わざ

・ロジスティック回帰のように「N個のうちy個で事象が生じる確率」を明示的にあつかう二項分布を使うことで、割算値の使用は回避できる

・人口密度のような「単位面積あたりの個体」みたいなのは個体数/面積とか使われがちだが、このときオフセット項わざ使える

・人口密度を、xiの関数として定義し、人口密度 = exp(β1 + β2xi)とする

・以下のようにモデル化できる

 平均個体数λi = 面積Ai × 人口密度 = 面積Ai × exp(β1 + β2xi)

・変形すると、以下

 λi = exp(β1 + β2xi + log面積Ai)

・これで、zi = β1 + β2xi + logAiを線形予測子とする対数リンク関数になる

・logAiのように、パラメーターがつかないlogAiのような項をオフセッット項とよぶ。線形予測子にlogAiといゲタをはかせているいめーじ

・Rでは、オフセット項以下のように指定可能

  glm(y ~ x, offset = log(A), family = poisson, data =d)

 

6.7 正規分布とその尤度

正規分布は連続値のデータを統計モデルで扱うための確率分布

f:id:takuro52006:20170503230018p:plain

ポアソン分布と同じく平均値のパラメーターμをもつ

ポアソン分布とは違って標準偏差のパラメーターσでデータのばらつきも指定できる

・Rで図示する場合は以下

 > y <- seq(-5, 5, 0.1)

 >plot (y, dnorm(y, mean = 0, sd = 1) type = "l" )

・Rで確率密度を計算した場合は、以下

 pnorm(1.8, 0, 1)  - pnorm(1.2, 0, 1)

正規分布最尤推定について

 ・N人あらなる集団の身長データをY = {yi}とする。iさんの身長がyi。

 ・あるyiがyi -0.5Δy ≦y ≦yi + 0.5Δyである確率は、確率密度関数p(yi|μ,σ)と区間幅Δyの積であると近似できる。

 ・従って、正規分布つかった尤度関数は、以下となる   

    f:id:takuro52006:20170505124036p:plain

 ・また、対数尤度関数は log  L (μ,σ) とかける

 

6.8 ガンマ分布のGLM

・ガンマ分布は確率変数のとりうる範囲が0以上の連続確率分布

f:id:takuro52006:20170503234450p:plain

・s= 1のとき、指数分布になる

・平均はs / r、分散はs / r**2となるので、分散 = 平均 / r

・期間1 / rごとに1回おこるランダムな事象がs回起こるまでの時間の分布を表す

・架空植物50個体の葉の重量xiと花の重量yiの関係を調べる。

・花の重量yiが、平均uiのガンマ分布にしたがっている

・何らかの理由があって以下を仮定 

 ui = Axi**b

・右辺でA = exp(a)とおいて、全体を指数関数でまとめる

 ui = exp( a ) xi**b = exp(a + blogxi)

・両辺の対数をとると、以下となる

 logμi = a + b logχi

・これで対数リンク関数と線形予測子a + b logxiで平均μiが与えられた。なお、線形予測子の説明変数はxiではなくlogxiとなり、推定すべきパラメタは切片aと傾きb

・Rのglm()でパラメーター推定するときは、以下

 glm(y  ~ log(x), family = Gamma(link = "log"), data = d)