m11.3

d<-Book%>%dplyr::select(Included2017,ADG,BW_start,Country)
table(d$Country, d$Included2017)
d$Included2017<-as.factor(d$Included2017)
levels(d$Included2017)<-list("1"="Yes", "0"="No")
fit<-glm(Included2017~ ADG + BW_start + Country, 
         data=d, 
         family="binomial",x=TRUE, y=TRUE)
summary(fit)
d$Included2017<-as.numeric(d$Included2017)-1
d<-d[complete.cases(d), ]
d$Country<-as.factor(d$Country)

m11.3 <- quap(
  alist(
    Included2017 ~ dbinom( 1 , p ) ,
    logit(p) <- a + b*ADG +c*BW_start + d[Country] ,
    a            ~ dnorm( 2  , 1 ),
    b            ~ dnorm( 0.5  , 1 ),
    c            ~ dnorm( 0.5  , 1 ),
    d[Country] ~ dnorm( 3  , 1 )
  ) , data=d)
set.seed(1999)
prior <- extract.prior( m11.3 , n=1e4 )
p <- sapply(1:5, function(k) inv_logit( prior$a + 
                                            prior$b +
                                            prior$c +
                                            prior$d[,k] ) )
dens( p , adj=0.1 )
mean( abs( p[,1] - p[,2] ) ) # Country 1 vs 2 
dens( abs( p[,1] - p[,2] ) , adj=0.1 )
dens( abs( p[,2] - p[,3] ) , adj=0.1 )
precis(m11.3, depth=2)
plot(precis(m11.3, depth=2))
post <- extract.samples(m11.3)
par(mfrow = c(4, 2))
dens(post$a)
dens(post$b)
dens(post$c)
dens(post$d[1,])
dens(post$d[2,])
dens(post$d[3,])
dens(post$d[4,])
dens(post$d[5,])
diffs <- list(
  db12 = post$d[,1] - post$d[,2],
  db13 = post$d[,1] - post$d[,3],
  db14 = post$d[,1] - post$d[,4],
  db15 = post$d[,1] - post$d[,5])
plot(rethinking::precis(diffs) )
mean( exp(post$d[,1]-post$d[,2]) )
diff_d <- post$d[,1] - post$d[,2]
diff_p <- inv_logit(post$d[,1]) - inv_logit(post$d[,2])
precis( list( diff_d=diff_d , diff_p=diff_p ) )
post <- extract.samples(m11.2)
p_incl <- inv_logit( post$a )
plot( precis( as.data.frame(p_incl) ) , xlim=c(0,1) )