[[radvance]]

[[ロジスティック回帰(ロジット分析)>http://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%82%B9%E3%83%86%E3%82%A3%E3%83%83%E3%82%AF%E5%9B%9E%E5%B8%B0]]とプロビット分析は、被説明変数が二者択一になる状況を表すのによく使われるモデルです。

　例えば世帯の年収と、自家用車を1台以上持っているかどうかを示すデータが何百世帯分かあるとします。自家用車を持っている=1、持っていない=0であらわすものとします。

　年収が高い世帯ほど、自家用車を持っていそうですよね。でも100%持っているわけではありません。年収が低くても、仕事に要るからと自動車をもっている人もいるでしょう。

+ x1 + x2 + x3 + …が大きくなるとyはだんだん大きくなって1に近づく。
+ x1 + x2 + x3 + …が小さくなるとyはだんだん小さくなって0に近づく。
+ 0 < y < 1を保ち、0以下や1以上の値をとらない。

　となるような関数y=f(z)があれば、z = a + b1x1 + b2x2 …という普通の回帰分析に似た式を経由して、y = a + b1x1 + b2x2 …と書くよりリアルな関係を分析できそうです。もっともらしく滑らかで平均点に対し対称形な関数として、ロジスティック関数とプロビット関数(正規分布累積関数の逆関数)がよく使われます。このふたつの関数形それぞれを仮定して、zとx1、x2…からb1、b2…を求める分析をロジット分析、プロビット分析といいます。

　以下ではプロビット分析の例を示しますが、ロジット分析の場合はprobitをlogitに置き換えればほぼ同じ結果になります。

　Rでふたつの分析を行う関数はglmで、R-2.8.0ではすでに基本パッケージに含まれていて、指名して読み込む必要もありません。たぶんそれ以降のバージョンもそうでしょう。ただ式の当てはまりを判断するための追加機能がパッケージpsclに含まれています。2008年以降の並河の講義でR-2.8.0以降を受け取った人は、その中に含まれています。自分でダウンロードした人は、Rの公式サイトからpsclパッケージをダウンロードする必要があります。

　以下のサンプルプログラムで説明します。

 library("pscl")
 
 x1 <- c( 1.5, 2.3, 3.8,   4.2,  5.6,   6.3)
 x2 <- c( 1.0, 3.9, 8.6, 15.25, 28.6, 32.68)
 x3 <- c( 8.5, 4.0, 9.2,   6.3,  6.8,   6.0)
 y  <- c(   1,   0,   1,     1,    0,     1)
 
 plot(z,x3)
 eq1 <- glm(y ~ x1 + x2, family=binomial(probit))
 summary (eq1)
 pR2(eq1)
 
 eq2 <- glm(y ~ x2 + x3, family=binomial(probit))
 summary (eq2)
 pR2(eq2)

　先頭1行は本当は実行のたび繰り返す必要はありません。読者の皆さんはcsvファイルなどを使って、別のファイルからデータを読み込むのが普通でしょう。その方法については初級的なコンテンツを見てください。

　実行結果の先頭部分はこうなっています。

 Call:
 glm(formula = y ~ x1 + x2, family = binomial(probit))
 
 Deviance Residuals: 
       1        2        3        4        5        6  
  1.2157  -1.5076   0.3398   0.6813  -1.2359   1.0224  
 
 Coefficients:
             Estimate Std. Error z value Pr(>|z|)
 (Intercept)  -1.9142     3.1309  -0.611    0.541
 x1            1.3711     1.7316   0.792    0.428
 x2           -0.1985     0.2483  -0.800    0.424
 
 (Dispersion parameter for binomial family taken to be 1)
 
     Null deviance: 7.6382  on 5  degrees of freedom
 Residual deviance: 6.9028  on 3  degrees of freedom
 AIC: 12.903
 
 Number of Fisher Scoring iterations: 6
　
　z統計量は普通の回帰分析でt値に相当するもので、Pr(>|z|)は各係数がそれだけ0から離れている確率を示します。x1もx2もあまり説明力はないようです。

　係数の解釈には注意してください。x1が1増えても、1.3711増えるのは上の説明中でのzであってyではありません。確率が何%上がるかはプロビット関数にzの値を代入して求める必要があります。

　さて、

     Null deviance: 7.6382  on 5  degrees of freedom
 Residual deviance: 6.9028  on 3  degrees of freedom
 AIC: 12.903

　ですが、式全体の当てはまりを言っているのだろうと見当はつきます。ヘルプファイルによると、glmの推定にはiteratively reweighted least squares (IWLS)と呼ばれるアルゴリズムを使っています。これはIWLSでなくIRLSと略されることのほうが多いようで、すぐに[[説明が見つかります>http://en.wikipedia.org/wiki/Iteratively_re-weighted_least_squares]]。[[日本語がよければこっち>http://ibisforest.org/index.php?IRLS%E6%B3%95]]。この繰り返しアルゴリズムで尤度関数を最大化し、b1、b2…の最尤推定量を求めています。

　AICは[[赤池情報量基準>http://ja.wikipedia.org/wiki/%E8%B5%A4%E6%B1%A0%E6%83%85%E5%A0%B1%E9%87%8F%E8%A6%8F%E6%BA%96]]で、同じような式をいろいろ推定したときは、AICがいちばん高いものを選びます。ただR2のように、1にどれくらい近いかで当てはまりのよさを判断する使い方はできません。

　ふたつのdevianceを尤度比検定して、b1、b2…がすべて0であるという仮説をカイ2乗検定することができます。[[具体的な方法はこちらを参照してください>http://www.ats.ucla.edu/stat/R/dae/probit.htm]]。しかし多くの統計パッケージでは、決定係数に似た当てはまりのよさを判断する数著出力しているようです。

　Rの場合、そうした機能はpsclのpR2関数に入っています。eq1への出力は次の通りです。

         llh     llhNull          G2    McFadden        r2ML        r2CU 
 -3.45142012 -3.81908501  0.73532978  0.09627041  0.11534272  0.16019432 

　NcFaddenの擬似R2(McFadden's pseudo r-squared)がよく判断に使われるようです。これは0から1の間の値を取ります。

　x3はyとほどほどに相関を持つよう作ってあります。サンプルを走らせるとグラフが出てくるので確認してください。x2、x3でyを説明するプロビットモデルのeq2は、McFadden's pseudo r-squaredが0.34とぐっと高くなり、x3のz統計量は1.232となります。
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS