久保みどり本 7章 途中まで。

midori_study_7.r

library("magrittr")
library("dplyr")
library("ggplot2")

# midori_study 7
# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html

rsc.dir <- "~/workspace/kubo_midori_7/"
d <- paste0(rsc.dir,"data.csv") %>% read.csv
head(d)
summary(d)

d4 <- d[d$x ==4,]
dim(d4)
table(d4$y)

qplot(data=d,x,y,geom="violin",group=x)+
  xlab("生存種子数 y_i")+
  ylab("葉数 x_i")+
  geom_jitter(alpha=.5,size=I(10),position = position_jitter(height = 0,width=.5)) # データの分布

# glmによるモデリング
model.glm <- glm(cbind(y,N-y)~x,family = binomial,data=d)
# 予測
d4$pred <- predict(model.glm,d4,type="response")

logistic <- function(z){
  1/(1+exp(-z))
}
logit_q <- function(b1,b2,x){
  b1 + b2*x
}
logit_q_rand <- function(b1,b2,x,r){
  b1 + b2*x + r
}

# rbinom(n = 20, prob = q,size = 8) %>% hist
# pbinom(prob = q,q=0:8,size = 8) %>% plot
# dbinom(0:8,size = 8,prob=q) %>% plot(type="b")

q <- logistic(logit_q(-2.148,0.510,x=4))
expected.df <- data.frame(
  y=0:8,
  prob=dbinom(0:8,size = 8,prob=q),
  N=20
)

# x_i=4
d4 %>%
  group_by(y,x) %>%
  summarize(count=length(x)) -> plot.d4
qplot(data=plot.d4,x=y,y=count)+geom_point(col="white",size=I(10))+
geom_line(data=expected.df,aes(x=y,y=prob*N))

x <- seq(0,8,0.1)
length(x)
z.upr <- logit_q_rand(-2,.7,x, 1)
z.mid <- logit_q_rand(-2,.7,x, 0)
z.lwr <- logit_q_rand(-2,.7,x,-1)
z.ptn.df <- data.frame(
  x=rep(x,3),
  z=c(z.upr,z.mid,z.lwr),
  r=c(rep("r_i > 0",length(x)), rep("r_i = 0 ",length(x)), rep("r_i < 0",length(x)))
)
# pict 7.5 p.152
qplot(data=z.ptn.df,x=x,y=logistic(z),geom="line",col=r)


# p.159

library("glmmML")
model.glmm <- glmmML(cbind(y, N-y)~x, data=d,family = binomial,cluster = id)
model.glmm$sigma
model.glmm$sigma.sd
model.glmm$coefficients

# http://qiita.com/HirofumiYashima/items/173ae2aacbaa339f75b8

calcProb <- function(x, b, r){
  # 生存確率の算出
  1.0 / (1.0 + exp(-1 * (b[1] + b[2] * x + r)))
}

calcL <- function(ylist, xfix, n, b, s){
  # 葉数 x を固定した場合の生存種子数 y の確率分布を算出
  sapply(ylist, function(y) {
   integrate(
      f = function(r) dbinom(y, n, calcProb(xfix, b, r)) * dnorm(r, 0, s),
      lower = s * -10,
      upper = s * 10
    )$value
  })
}

sigma <-model.glmm$sigma
beta <- model.glmm$coefficients
# 生存種子数 y と個体数の分布
yy <- 0:8
plot(yy, table(d[d$x == 4,]$y), xlab="y", ylab="num")
# 予測線
lines(yy, calcL(yy, 4, max(d$N), beta, sigma) * length(d[d$x == 4,]$y), col="red", type="b")