0_データ解析のための統計データモデリング ロジスティック回帰をggplotで.md

データ解析のための統計データモデリング　ロジスティック回帰をggplotで

• 書籍
  • http://www.amazon.co.jp/dp/400006973X
    • 6.4 ロジスティック回帰とロジットリンク関数
• ウェブ
  • http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html

logistic <- function(z) 1/(1 + exp(-z))
z <- seq(-6, 6, 0.1)
# plot(z,logistic(z),type='l')
library("ggplot2")
qplot(x = z, y = logistic(z)) + geom_line()

# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
d <- read.csv("./data4a.csv")
qplot(data = d, x = x, y = y, col = f)

head(d)

##   N y     x f
## 1 8 1  9.76 C
## 2 8 6 10.48 C
## 3 8 5 10.83 C
## 4 8 6 10.94 C
## 5 8 1  9.37 C
## 6 8 1  8.81 C

head(cbind(d$y, d$N - d$y))

##      [,1] [,2]
## [1,]    1    7
## [2,]    6    2
## [3,]    5    3
## [4,]    6    2
## [5,]    1    7
## [6,]    1    7

glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)

## 
## Call:  glm(formula = cbind(y, N - y) ~ x + f, family = binomial, data = d)
## 
## Coefficients:
## (Intercept)            x           fT  
##      -19.54         1.95         2.02  
## 
## Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
## Null Deviance:	    499 
## Residual Deviance: 123 	AIC: 272

dim(d)

## [1] 100   4

train.row <- sample(100, 50)
d.train <- d[train.row, ]
d.test <- d[-train.row, ]

d.glm <- glm(cbind(y, N - y) ~ x + f, data = d.train, family = binomial)
print(d.glm)

## 
## Call:  glm(formula = cbind(y, N - y) ~ x + f, family = binomial, data = d.train)
## 
## Coefficients:
## (Intercept)            x           fT  
##      -19.71         1.97         1.87  
## 
## Degrees of Freedom: 49 Total (i.e. Null);  47 Residual
## Null Deviance:	    266 
## Residual Deviance: 80.6 	AIC: 157

# 残り半分のデータの予測。type='response'
# 付けないとデフォルトで線形予測値ななる模様 ?predict.glm
d.pred <- predict(d.glm, d.test[, c(3, 4)], type = "response", se = T)
d.res <- data.frame(d.test, pred = d.pred$fit * d$N, ucl = (d.pred$fit + d.pred$se.fit) * 
    d$N, lcl = (d.pred$fit - d.pred$se.fit) * d$N)

## Warning: row names were found from a short variable and have been
## discarded

library("ggplot2")
p <- qplot(y = y, x = x, col = f, data = d.res)
# p<-p+facet_grid(~f)
p <- p + geom_smooth(aes(ymax = ucl, ymin = lcl), data = d.res, stat = "identity")
# 勝手に線引かれてしまうのはなぜ？
p

data4a.csv

rpubsにアップロードできなかったのでひとまずここに。