ロジスティック回帰

かなり前(最終更新日2012-08-12)にざっと書いて放置していたページです。とりあえずMathJaxで書き直しました。より新しいロジスティック回帰と変数選択のほうが詳しいと思います。

合格の判定

問題

次のような表があるとします。

$y$$x_1$$x_2$
03.660.1
14.152.0
03.762.5
04.960.6
14.454.1
03.663.6
14.968.0
04.838.2
14.159.5
04.347.3

ここで $y$ は入試の合否(合格なら1,不合格なら0),$x_1$ は内申点,$x_2$ は模試偏差値を表すとしましょう。

通常の重回帰分析では,

\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 \]

が(最小2乗法の意味で)できるだけ正確に成り立つように $\beta_0$,$\beta_1$,$\beta_2$ を求めます。しかし,合否は0か1しか値がありません。そこで,1である確率を $\pi$ と書けば,ロジスティック回帰では次の式が(最尤法の意味で)できるだけ正確に成り立つような $\beta_0$,$\beta_1$,$\beta_2$ を求めます:

\[ \log\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 \]

上のような回帰式の左辺に現れる関数をリンク関数といいます。この場合のリンク関数はロジット関数と呼ばれ,$\mathrm{logit}(\pi)$ と書くことがあります。ちなみに $\pi/(1-\pi)$ はオッズ(odds)と呼ばれる量ですので,ロジット関数はオッズの対数ということになります。

上の式を確率 $\pi$ について解けば次のようになります:

\[ \pi = \frac{\exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2)}{1 + \exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2)} = \frac{1}{1 + \exp[-(\beta_0 + \beta_1 x_1 + \beta_2 x_2)]} \]

この右辺は $\mathrm{logit}^{-1}(\beta_0 + \beta_1 x_1 + \beta_2 x_2)$ とも書きます。この $\mathrm{logit}^{-1}$ という関数(logitの逆関数)はロジスティック関数と呼ばれるものの一種です。

データの読み込み方(RStudio使用)

data1.csv というファイルをダウンロードして,RStudioの右上の「Workspace」ペインで「Import Dataset」→「From Text File...」とたどって,読み込ませます。

データの読み込み方(クリップボード経由)

  1. WindowsならRに次のように打ち込んでください(リターンキーはまだ押しません):
    data1 = read.table("clipboard", header=TRUE)
    MacならRに次のように打ち込んでください(リターンキーはまだ押しません):
    data1 = read.table(pipe("pbpaste"), header=TRUE)
  2. このページの上の表または同じデータの入ったExcelの表を(ヘッダ部分「y,x1,x2」も含めて)マウスで選択し,右クリックで「コピー」を選びます。
  3. Rで,さきほど押さなかったエンターキーを押します。
  4. 念のため,Rで data1 と打ち込んで,正しいデータになっていることを確かめます。

データの読み込み方(Excel経由)

  1. Excelにデータを打ち込み,CSV形式で保存してください。
  2. Rで次のように打ち込んでください:
    data1 = read.csv(file.choose())
    ファイル選択の窓が現れますので,さきほど保存したCSVファイルを選んでください。Rで「不完全な最終行が見つかりました」という警告が出ることがありますが,無視してください。
    もし日本語を含む表が文字化けするなら次のようにして文字コードを指定します:
    data1 = read.csv(file.choose(), fileEncoding="CP932")
  3. 念のため,Rで data1 と打ち込んで,正しいデータになっていることを確かめます。

通常の重回帰分析

Rで重回帰分析をするには,

lm(y ~ x1 + x2, data=data1)

と打ち込みます(lm は linear model の意味です)。もっと良い方法は,

result1 = lm(y ~ x1 + x2, data=data1)

のように結果を変数(オブジェクト)に入れておき,

summary(result1)

として結果を表示させることです。このほうが複数のモデルを比較するときにも便利です。この場合は次のように表示されます:

Call:
lm(formula = y ~ x1 + x2, data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6548 -0.3092 -0.2713  0.5010  0.7095 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.69441    2.24318  -0.755    0.475
x1           0.29598    0.37222   0.795    0.453
x2           0.01483    0.02159   0.687    0.514

Residual standard error: 0.552 on 7 degrees of freedom
Multiple R-squared: 0.1113,	Adjusted R-squared: -0.1426 
F-statistic: 0.4385 on 2 and 7 DF,  p-value: 0.6616 

つまり $y = -1.69441 + 0.29598x_1 + 0.01483x_2$ で予測されることがわかります(ただしどの係数も有意ではありません)。

ロジスティック回帰

上で説明したロジスティック回帰を行って結果を result2 というオブジェクトに代入するには

result2 = glm(y ~ x1 + x2, data=data1, family=binomial(link="logit"))

と打ち込みます。続いて

summary(result2)

と打ち込めば次のように表示されます:

Call:
glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), 
    data = data1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4754  -0.8584  -0.8007   1.1905   1.5719  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.44589    9.12237  -1.035    0.300
x1           1.27158    1.49423   0.851    0.395
x2           0.06424    0.08739   0.735    0.462

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 13.460  on 9  degrees of freedom
Residual deviance: 12.345  on 7  degrees of freedom
AIC: 18.345

Number of Fisher Scoring iterations: 4

つまり $\mathrm{logit}(\pi) = -9.44589 + 1.27158x_1 + 0.06424x_2$ で予測されることがわかります(ただしどの係数も有意ではありません)。

また,fitted(result2) と打ち込むと,フィットされた値 $\mathrm{logit}^{-1}(\beta_0 + \beta_1 x_1 + \beta_2 x_2)$ が表示されます。

カブトムシの問題

次の表は,Annette J. Dobson and Adrian G. Barnett, An Introduction to Generalized Linear Models, 3rd ed. (CRC Press, 2008) の p.127 にあるデータです。$x$ はある薬品の量(対数),$n$ はその薬品を与えたカブトムシの数,$y$ はそのうち死んだ数です。薬品の量が増えると死ぬ確率がどのように高くなるかを調べます。

$x$$n$$y$
1.690759 6
1.72426013
1.75526218
1.78425628
1.81136352
1.83695953
1.86106261
1.88396060

モデルは $\mathrm{logit}(\pi) = \beta_0 + \beta_1 x$ で,確率 $\pi$ はほぼ $y/n$ になるはずです。

さきほどと同様にしてデータ(data2.csv)を data2 に読み込んで,次のように打ち込みます。

result3 = glm(cbind(y,n-y) ~ x, data=data2, family=binomial(link="logit"))

関数 glm() に与えるモデル式の左辺は,成功した個数のベクトル $y$ と失敗した個数のベクトル $n - y$ をコラム結合 cbind() した8行2列の行列です。

summary(result3) と打つと次のような結果が表示されます:

Call:
glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"), 
    data = data2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5941  -0.3944   0.8329   1.2592   1.5940  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -60.717      5.181  -11.72   <2e-16 ***
x             34.270      2.912   11.77   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 284.202  on 7  degrees of freedom
Residual deviance:  11.232  on 6  degrees of freedom
AIC: 41.43

Number of Fisher Scoring iterations: 4

その他のモデル

関数 glm() の名前は一般化線形モデル(generalized linear model,GLM)から来ています。family= で指定できる確率分布およびリンク関数のおもなものを挙げておきます:

詳しくは glmfamily のヘルプをご覧ください。

glmが収束しないとき

より収束しやすい glm2 というパッケージが開発されています。詳しくは The R Journal, Volume 3/2 の Ian C. Marschner, “glm2: Fitting Generalized Linear Models with Convergence Problems” というペーパーをご覧ください(オープンアクセスです)。

多項ロジスティック回帰

Rで多項ロジスティック回帰をするパッケージはいくつかあるようですが,ここではVGAMパッケージをご紹介しておきます。詳しくは Thomas W. Yee, The VGAM Package for Categorical Data Analysis, Journal of Statistical Software, Vol.32, No.10 (2010) をご覧ください。


Last modified: