【二項分布】ロジスティック回帰について|R
GLMモデルと呼ばれる、一般化線形モデルにはポアソン分布だけではなく、色々の分布が使われます。
今回は二項分布を使った「ロジスティック回帰」について解説いたします。
Rではなく、pythonでコードを見たい方は【分類タスク】ロジスティック回帰の使い方|pythonをご覧ください。
ロジット関数とロジスティック関数
かんたんにGLMの復習を行います。
しっかり復習したい方は、【GLM】一般化線形モデルを解説(統計モデル編)をご覧ください。
ある分布に従う確率変数Yがあるとして、その分布の平均をμとします。
μを変形した形f(μ)にすると、目的変数を線形で表すことができます。
これを一般化線形モデル(GLM)と呼びます。
今回は、f(μ)がロジットリンク関数であるものを扱います。
一般化線形モデルには、3つの特徴がありました。
$$f(θ_{n})=β_{0}+β_{1}X_{1}+…+β_{k}X_{k}$$
①線形予測子
上の式で言うと、右辺にあたります。入力データXの数だけ、パラメーターをかけたものです。
②確率分布
パラメーターβが得られた時のyの分布
③リンク関数
上の式で言うと、左辺のことです。線形予測子が作る関数のことです。Yが形を変えたものと考えてもらって構いません。
確率分布がポアソン分布の時のGLMは、リンク関数はλを対数化したものになります。下のような形でした。
$$logλ_{i}=β_{0}+β_{1}X_{1}+β_{2}X_{2}$$
まずは、今回扱うロジスティック回帰の特徴を説明していきます。
①確率分布:二項分布
②リンク関数:ロジットリンク関数
このようになります。注意として、カウントデータが1個のみなら①の確率分布はベルヌーイ分布になります。
二項分布(binomial distribution)
二項分布とは、独立なベルヌーイ分布をn回繰り返した時の和の分布です。
実現値が有限である点において、ポアソン分布と異なります。
【例】コインを5回投げて表の出た回数をXとすると、実現値は0,1,2,3,4,5の6種類です。表が出る確率を仮にpとすると、実現値0と実現値1は以下のように表せます。
$$(1-p)^5$$
$$5p(1-p)^4$$
このようになっていきます。実現値0から5までの確率の和は以下のようになります。
一般化すると、二項分布の確率密度関数は以下のようになります。nは実現値です。
ロジットリンク関数(logit link function)
先ほどの二項分布と同じコインの例を挙げます。
ごちゃごちゃ説明する前に、コードをみてください。
z<-seq(-6,6,0.1)
#線形予測子/右辺
logit<-function(z) z
#ロジット関数の定義
#リンク関数/左辺
plot(z,logit(z))
#作図
seq(a,b,c)というのは、最小値aで最大値bまでをc刻みで代入するものです。
$$log \frac{q_{i}}{1-q_{i}}=z_{i}$$
上の式がロジット関数です。
右辺は線形予測子、左辺はオッズを対数化したものです。コインの例で言うと、下のようになります。
$$log \frac{コインがi回表の確率}{コインがn-i回裏の確率}=z_{i}$$
対数オッズと、オッズの関係はこのようになります。上側の右辺がよくみる線形予測子の形です。
しかし、出力結果が、「オッズを対数化したもの」だと目的に沿わない場合もあります。
ちゃんと、「コインがi回表になる確率を知りたい!」と言う場合には、ロジット関数をロジスティック関数にする必要があります。
補足:指数型分布族
リンク関数をロジット関数にすると、尤度関数が指数型分布族になるというメリットがあります。
以下のようなeの指数関数の形で確率密度を表せるものは、すべて指数型分布族と呼びます。
例えば、正規分布やガンマ分布、ポアソン分布などの主要な確率分布は全て指数型分布族です。
$$f(x;θ)=ae^{θ^TS(x)-ψ(θ)}$$
例えば、\(N(μ,1)\)の正規分布を指数型分布族の形に直してみます。
$$f(x;μ)=\frac{1}{\sqrt{2π}}e^{-\frac{(x-μ)^2}{2}}=\frac{1}{\sqrt{2π}}e^{-\frac{x^2}{2}+μx-\frac{μ^2}{2}}$$
ここでいうθとは、期待値μのことなので、以下のようになります。
aはxの関数です。aとeの積では、指数同士を足し算しています。
$$a=\frac{1}{\sqrt{2π}}e^{-\frac{x^2}{2}},S(x)=x,ψ(θ)=\frac{θ^2}{2}$$
このような形だと何が嬉しいのかというと、\(ψ(θ)\)が狭義の凸関数になるということです。
つまり、この確率密度関数の対数尤度関数は凹関数(上に凸)になります。
よって尤度方程式の解は(存在するならば)一意的で,それが最尤推定量となる。
このように,指数型分布族 は比較的扱いやすい統計モデルです。(停留点を求めると、それが最尤推定量になるということ)
最尤推定量については、以下をご覧くださいませ。
ロジスティック関数(logistic function)
ロジット関数の逆関数が、ロジスティック関数です。
ロジスティック関数は以下のような形になります。
$$logistic(z_{I})=\frac{1}{1+exp(-z_{i})}$$
線形予測子Ziに−1がかかりました。eの-線形予測子乗に1を足したものが分母となっています。
Rで書いてみましょう。
logistic <- function(z) 1 /(1+exp(-z))
#ロジスティック関数の定義
plot(z,logistic(z))
#作図
ロジット関数と比べて、滑らかになりました。「コインが表になる確率」が上のように推移することがわかりました。
オッズ比に関して詳しく知りたい方は、【ベイズ因子】オッズ比を実例を通して紹介をご覧ください。
【実践】健康寿命と喫煙、肺の大きさ
ロジット関数の意味が大体理解できたと思うので、Rで実例を扱ってみます。
story
2055年に亡くなった方の健康寿命を2040年時点から遡ってみました。 中でも、喫煙経験があった25人と25人の計50人を調査し、肺の大きさも測ります。 N:2040年から2055年までの年数。当然全員15年。 y:2040年から測った健康寿命。15以下。 x:腎臓の長さ s:喫煙因子/喫煙してるかしてないかの二択。 N:15 y:10 x:10.43 s:N 例:2050年までは健康に生きれたが2055年に亡くなった。肺の横の長さは10.43cm。喫煙経験はない人。
health<-read.csv("health.csv")
#データの読み込み
health$y
#健康寿命を表示
健康寿命をとりあえず出力してみました。平均すると2040年から7年くらいですね。
以下は簡単にグラフ化したものです。白の方が喫煙経験なしなので、ある程度健康寿命が伸びそうな傾向が見えます。
では分析を行います。
glm(cbind(y,N-y)~x+s,data=health,family=binomial)
ばらつきの確率分布は、family=binomialで二項分布を使います。
応答変数に関しては、cbind(y,N-y)となっていますが、これはcbind(健康に生きれた年数、残りの健康に生きれなかった年数)を表しています。
推定値としては、以上のようになります。
ロジスティック回帰からわかること
上の例で言うと、腎臓の大きさによる推定パラメータが0.17039でした。
つまり、個体iの腎臓の大きさが「1単位」増大すると、健康寿命のオッズは、
exp(0.17039)倍になることがわかります。
電卓で簡単に計算できますが、1.185767倍でした。腎臓が大きい人は、ちょっと健康寿命が大きい傾向にあることがわかりました。
*二項分布を使ったGLMで使われるリンク関数は、ロジットリンク関数だけでなく、プロビット関数など別の関数もあります。
ロジスティック回帰の場合には、リンク関数は、ロジット関数を使うというだけの話です。