ロジスティック回帰分析

MikuHatsune2013-05-28

ロジスティック回帰分析という解析があるが、聞いたことはあっても実はやったことがなかったのでやってみる。
ロジスティックの式は、昔もやったが、微分方程式
\frac{dN}{dt}=r(\frac{K-N}{N})N
を解いて
N(t)=\frac{1}{1+e^{rK(t_0-t)}}
となる。
解析するときは、線形回帰として
log(\frac{p}{1-p})=\alpha + \beta x
を解く。これなら直線でプロットできて、これを確率pについて解けば
p(t)=\frac{1}{1+e^{-(\alpha + \beta x)}}
というシグモイドカーブが得られる。
 
試験の得点と合否のデータを取得して、得点と合格率をロジスティック回帰分析する。

a <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
g0 <- glm(admit ~ gre, data=a)
g0$coefficients
  (Intercept)           gre 
-0.1198407147  0.0007441564 
# 係数は、線形回帰
# シグモイドカーブ関数に変換する
p <- function(x){
	1/(1+exp(-sum(g0$coefficients * c(1, x))))
}

s0 <- seq(-5000, 5000, length=1000)
plot(a$gre, a$admit, xlim=range(s0))
points(s0, mapply(p, s0), type="l", col=2, lwd=3)

x軸の幅の都合で、ほぼ直線のように見えるが、幅を広げるとシグモイドカーブになっている。


 
と思っていたら、

glmでロジスティック回帰は family=binomial 指定が必要です

という指摘を受けたので、glmのヘルプを見ると

Usage:

     glm(formula, family = gaussian, data, weights, subset,
         na.action, start = NULL, etastart, mustart, offset,
         control = list(...), model = TRUE, method = "glm.fit",
         x = FALSE, y = TRUE, contrasts = NULL, ...)
     
Arguments:

  family: a description of the error distribution and link function to
          be used in the model. This can be a character string naming a
          family function, a family function or the result of a call to
          a family function.  (See ‘family’ for details of family
          functions.)

デフォルトで family = gaussian, family引数は何かというと、誤差の分布のようである。
合否の二値なので binomial を指定しなければならない。これはLASSOでもはまった。

g0 <- glm(admit ~ gre, data=a, family=binomial)
s0 <- seq(0, 1500, length=1000)
plot(a$gre, a$admit, xlim=range(s0))
points(s0, mapply(p, s0), type="l", col=2, lwd=3)