ロジスティック回帰分析という解析があるが、聞いたことはあっても実はやったことがなかったのでやってみる。
ロジスティックの式は、昔もやったが、微分方程式
を解いて
となる。
解析するときは、線形回帰として
を解く。これなら直線でプロットできて、これを確率について解けば
というシグモイドカーブが得られる。
試験の得点と合否のデータを取得して、得点と合格率をロジスティック回帰分析する。
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)