ロジスティック回帰

かなり前(最終更新日2012-08-12)にざっと書いて放置していたページを少しずつ更新しています。ロジスティック回帰と変数選択もご覧ください。

合格の予測(1次元)

100点満点の模擬テストを20人が受験しました。点数は次の通りです:

x = c(4,5,17,17,19,19,29,41,45,47,49,54,56,58,70,76,88,88,93,97)

この20人がその後実際の試験を受験して,合否が判明しました。合 = 1,否 = 0 で表すと,次のようになりました:

y = c(0,0,0,0,0,0,1,1,1,0,0,1,1,0,1,1,1,1,1,1)

プロットしてみましょう:

plot(x, y, xlim=c(0,100), ylim=c(-0.2,1.2), yaxt="n")
axis(2, at=c(0,1), labels=c(0,1), las=1)

x から y を予測する式を求めたいのですが,通常の線形回帰ではうまくあてはまりません:

r = lm(y ~ x)
abline(r)

そもそもこれでは何を予測したいのかはっきりしません。予測するものを,模擬テストの点数が x のときの合格確率だとしましょう。確率は負になったり 1 を超えたりしませんので,線形回帰ではうまくありません。うまく 0 から 1 の範囲におさまってくれる関数はないでしょうか。

例えば

\[ y = \frac{1}{1 + \exp(-x)} = \frac{1}{1 + e^{-x}} \]

はどうでしょうか。これをグラフにすると,次のようになります:

curve(1 / (1 + exp(-x)), -5, 5)

これなら縦軸の値は必ず 0 から 1 の範囲に収まりますので,確率と解釈することができます。実際には,このグラフの中心(変曲点の位置)と曲がり具合とを変えられるようにした

\[ y = \frac{1}{1 + \exp(-(w_0 + w_1 x))} \]

という式で表される曲線を使うことにします。

この形の曲線をロジスティック曲線といいます。より一般に,0 から 1 に滑らかに変化する S 字形をつぶしたような曲線をシグモイド曲線といいます。シグモイド曲線は,ニューラルネットの活性化関数として使われてきました(現在は ReLU を使うのが普通です)。

$x$ が与えられたとき $y = 1$ になる確率 $p$ が

\[ p = \frac{1}{1 + \exp(-(w_0 + w_1 x))} \]

で近似できるような係数 $w_0$,$w_1$ を求める方法をロジスティック回帰といいます。この式を書き直すと

\[ \log\left(\frac{p}{1-p}\right) = w_0 + w_1 x \]

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

R ではロジスティック回帰は次のように打ち込むと簡単にできます:

r = glm(y ~ x, family=binomial(link="logit"))
summary(r)

結果は次のようになります:

Call:
glm(formula = y ~ x, family = binomial(link = "logit"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7368  -0.4784   0.1618   0.5249   1.6902  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -3.56668    1.68570  -2.116   0.0344 *
x            0.08319    0.03562   2.335   0.0195 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 27.526  on 19  degrees of freedom
Residual deviance: 15.187  on 18  degrees of freedom
AIC: 19.187

Number of Fisher Scoring iterations: 5

この結果の Coefficients の Estimate のところに係数 $w_0$,$w_1$ が出ます。つまり,

\[ P = \frac{1}{1 + \exp(-(-3.56668 + 0.08319 x))} \]

ということになります。元のデータとこの確率を重ねてグラフにしてみます:

plot(x, y, xlim=c(0,100), ylim=c(-0.2,1.2), yaxt="n")
axis(2, at=c(0,1), labels=c(0,1), las=1)
curve(1 / (1 + exp(-(0.08319 * x - 3.56668))), add=TRUE)

fitted(r) と打ち込めば,各人の合格確率が表示されます:

         1          2          3          4          5          6          7 
0.01575446 0.02773306 0.17147626 0.17147626 0.19543345 0.19543345 0.31521945 
         8          9         10         11         12         13         14 
0.45896264 0.50687704 0.53083424 0.55479144 0.61468444 0.63864164 0.66259884 
        15         16         17         18         19         20 
0.80634203 0.87821363 1.02195682 1.02195682 1.08184982 1.12976422 

合格の判定(2次元)

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

$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

さきほどと同様に,$x_1$,$x_2$ が与えられたとき $y = 1$ になる確率 $p$ が

\[ p = \frac{1}{1 + \exp(-(w_0 + w_1 x_1 + w_2 x_2))} \]

で近似できるようにパラメータ $w_1$,$w_2$ を定めます。

次のように打ち込むとデータが data1 に読み込まれます:

data1 = read.csv("https://oku.edu.mie-u.ac.jp/~okumura/stat/data1.csv")

このデータについてロジスティック回帰を行うには

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

と打ち込みます。続いて

summary(r2)

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

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(r2) と打ち込むと,フィットされた値 $\mathrm{logit}^{-1}(w_0 + w_1 x_1 + w_2 x_2)$ が表示されます。

あやめ(iris)の分類

R で「iris」と打ち込んでみてください。次のような150行のデータが現れます:

    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1            5.1         3.5          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
3            4.7         3.2          1.3         0.2     setosa
...

これは有名な Fisher の iris(あやめ)のデータです。3種類の iris(setosa,versicolor,virginica)それぞれ50個ずつについて,sepal length(萼(がく)の長さ),sepal width(萼(がく)の幅),petal length(花びらの長さ),petal width(花びらの幅)が与えられています(単位は cm)。

まずは何らかの形でプロットしてみましょう。例えば萼と花びらの長さをプロット:

plot(iris$Sepal.Length, iris$Petal.Length, col=iris$Species)
legend("bottomright", legend=levels(iris$Species), pch=1, col=1:3)

主成分分析をするのも手です:

r = prcomp(iris[,1:4])
plot(r$x[,1], r$x[,2], col=iris$Species)
legend("topright", legend=levels(iris$Species), pch=1, col=1:3)

ロジスティック回帰で,これら4つの変数から,まずは setosa かそれ以外かを判別してみましょう。

ヒント:iris というデータに setosa という列を追加して,Species が setosa なら 1,そうでないなら 0 とする命令:

iris$setosa = ifelse(iris$Species == "setosa", 1, 0)

残りは課題です。

カブトムシの問題

次の表は,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) = w_0 + w_1 x$ で,確率 $\pi$ はほぼ $y/n$ になるはずです。

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

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

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

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

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: