##############
#  R Plotting routine for dugongs Binary regression model
#  See Carlin and Louis (2008), Example 4.4
##############

xgrid <- seq(0.5,32,length=101)
age <- c(1.0, 1.5, 1.5, 1.5, 2.5, 4.0, 5.0, 5.0, 7.0, 8.0, 8.5, 9.0, 9.5, 9.5, 
  10.0, 12.0, 12.0, 13.0, 13.0, 14.5, 15.5, 15.5, 16.5, 17.0, 22.5, 29.0, 31.5)
Y <- c(1.80, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47,
       2.19, 2.26, 2.40, 2.39, 2.41, 2.50, 2.32, 2.32, 2.43,
       2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.70, 2.72, 2.57)
lgage <- log(xgrid)
ctlgage <- lgage - mean(lgage)
z <- as.integer(Y>2.4)

beta0 <- -1.52;  beta1 <- 6.19
p_logit <- exp(beta0 + beta1*ctlgage)/(1 + exp(beta0 + beta1*ctlgage))

beta0 <- -0.79;  beta1 <- 3.39
p_probit <- pnorm(beta0 + beta1*ctlgage)

beta0 <- -1.79;  beta1 <- 4.58
p_cloglog <- 1 - exp(-exp(beta0 + beta1*ctlgage))

plot(xgrid, p_logit, xlab="age", ylab="P(dugong is full grown)",
  pch=20,type="l")
lines(xgrid, p_probit, lty=2)
lines(xgrid, p_cloglog, lty=3)
points(age, z, pch=19)
legend(20, .6, legend=c("logit", "probit", "cloglog"), lty=1:3,
  ncol=1)