読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第50回 プロビット回帰分析

第50回はロジスティックとプロビットの関係について書いていきます。ロジスティック回帰分析のことを調べると、爾汝の交わりのようにつきまとってくるプロビットという言葉。「プロビット」って何者なのでしょうか。

 

・許容値分布

2値変数を扱う場合はロジスティック回帰分析を使うことが多いと思います。まず、プロビットを理解するために許容値分布(Tolerance Distribution)を説明します。

 

成功の確率を以下のように考えましょう。

f:id:teruaki-sugiura:20150822232409p:plain

確率を求めるために、日々このような式を使っていますね。この確率値を複数の変数を用いて分析したいので、一番簡単な線形モデルを作りましょう。

f:id:teruaki-sugiura:20150823122103p:plain

しかし、このままではの値によっては、確率pが0~1に収まらなくなるので、いろいろ式を変形して、0~1に収まるように工夫しなければいけないと、前回の記事でかきました。つまり、確率pを0~1に制限するため、以下の累積分布関数を用います。このとき、確率密度関数f(s)が、許容値分布と呼ばれています。

f:id:teruaki-sugiura:20150822232648p:plain

しれっ、と登場しましたが、累積分布関数って何者なのでしょうか。個人的には、統計学というよりも確率論の段階の講義で説明されるイメージですね。大学や学部によっては、「こんな確率分布ありますよー」で、さらっ、と流されて終わり。なんてこともあるみたいです。詳しく学びたい時は、YouTube筑波大学の「確率論」の講義が無料で観られますで、そちらでご覧ください。では、簡単に説明します。

 

確率密度関数 

離散型を飛ばして連続型の確率密度関数を説明します。確率密度関数をf(x)としましょう。天下り的ですが、「確率」密度関数なので、全域で積分すると「1」になります。言い換えると、「確率密度関数積分すると確率が得られる関数」のことです。

f:id:teruaki-sugiura:20150822232740p:plain

上に書いたような関数f(x)があるならば、f(x)を確率変数Xの「確率密度」と呼び、f(x)の分布に従ってXは発生していると考えることができます。連続型確率密度関数は以下のように定義されます。

f:id:teruaki-sugiura:20150822232835p:plain

・累積分布関数(分布関数)F(x)

次に、確率密度関数と関わりの深い「累積分布関数F(x)」の説明をします。この関数は、確率を求めるための関数であり、「確率変数Xがxという値をとる時、x以下の確率の値を教えてくれる関数」なのです。以下のように定義されます。

f:id:teruaki-sugiura:20150822232910p:plain

以下のような性質を持ちます。

f:id:teruaki-sugiura:20150822232926p:plain

確率密度関数と分布関数の関係

以下の図のような確率密度関数があったとします。そこで、確率密度関数f(x)を区間(−∞, x]で積分すると分布関数F(x)が得られます。つまり、「確率密度関数」を積分すると「分布関数」が得られます。反対に、「分布関数」を微分すると「確率密度関数」を得られという関係にあります。なので、確率密度関数の-∞からxまでの面積が分布関数F(x)ということになります。

f:id:teruaki-sugiura:20150822232959p:plainf:id:teruaki-sugiura:20150822233008j:plain

・プロビット回帰分析

例を使ってプロビット回帰分析を考えます。サッカー選手100人(各カテゴリ20人)を対象に、痛止めの新薬の投与量を変えながら、効果量を探っていったとします。結果を以下の図にまとめました。

 

f:id:teruaki-sugiura:20150822233041p:plain

0.1gのカテゴリで反応した(痛みが消えた選手)2人は痛止めが効くであろう量は、0g~0.1gの範囲の中にあるはずです。 1.0gで反応した3(5-2)人の痛止めが効くであろう量は0.1g~1.0gの範囲に、 10.0gで反応した8人の痛止めが効くであろう量は1.0g~10.0gの範囲にあると考えることができます。10.0gで効果があった8人、と考えていくと、痛止めが効くであろう量は0g~10.0gまでの累積数は13人となります。そう考えると、以下の右側のグラフが得られます(左側のグラフは効果率ではなく、効果人数ですね。申し訳ありません。)。以上より、投与量と反応率の関係は正規分布と累積分布関数のようにも見えますし、シグモイド曲線(緑点線)で近似できそうです。

 

従って、痛み止めが効くであろう量が近似的に正規分布するとき、投与量と効果率は、累積正規分布曲線で近似できます。投与量xと効果率pの関係は「累積分布関数」で表せ、投与量と効果率の関係を線形にするために、リンク関数を累積分布関数の逆関数を設定し、効果率を変換します。プロビット回帰分析では、正規分布を許容値分布として利用します。分析結果はロジスティック回帰分析と似ますが、このようなデータを分析する場合にプロビット回帰分析を行います。こんな感じのグラフになります。

 > curve(pnorm(x,mean=0,sd=1),xlim=c(-10,10),ylim=c(0,1),xlab="linear predictor",ylab="p",main="probit")

f:id:teruaki-sugiura:20150823164623p:plain

・プロビット回帰分析のモデル

 

[0,1]に抑えるため、以下のようにモデルを組みます。ここでは、確率密度関数f(s)で表し、効果率をpで表します。許容値分布は正規分布です。

f:id:teruaki-sugiura:20150823001318p:plain

Φは標準正規分布N(0,1)の分布関数とすると、以下のようになります。 

f:id:teruaki-sugiura:20150822233532p:plain

そして、標準正規分布N(0,1)の分布関数の逆関数を求め、

f:id:teruaki-sugiura:20150822233812p:plain

と置きます。すると、以下のようになります。

f:id:teruaki-sugiura:20150822233742p:plain

これは、リンク関数gを「標準正規分布の『累積分布関数の逆関数』」とする一般化線形モデルといえます。これがプロビット回帰分析の仕組みです。

 

以上でおしまい。

 

広告を非表示にする