【二項分布】ロジスティック回帰について|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で使われるリンク関数は、ロジットリンク関数だけでなく、プロビット関数など別の関数もあります。

ロジスティック回帰の場合には、リンク関数は、ロジット関数を使うというだけの話です。

FOLLOW ME !