# Exemplo 2.1, pagina 36
# Modelos Lineares Generalizados e Aplicacoes
# Cordeiro, Demetrio e Moral (2024)
# Editora Blucher 
# ABE - Projeto Fisher
#
# Descricao:  Os dados referem-se a um ensaio de toxicidade de rotenona,
# no delineamento completamente casualizado, em que doses (x_i) do 
# inseticida foram aplicadas a n_i insetos (pulgao de crisantemo) e, 
# apos um certo tempo, foi observado o numero (y_i) de insetos mortos.
#
# Objetivo:  O interesse do pesquisador estava na deteminacao das doses
# letais que matam 50% e 90% dos insetos, para recomendacao de aplicacao
# do inseticida no campo.
#

x = c(0,2.6,3.8,5.1,7.7,10.2)
n = c(49,50,48,46,49,50)
y = c(0,6,16,24,42,44)
p = y/n
cbind(x,n,y,p)

plot(x,p,ylim=c(0,1))

ns = sum(n)
ys = c(
rep(1,y[1]),rep(0,n[1]-y[1]),
rep(1,y[2]),rep(0,n[2]-y[2]),
rep(1,y[3]),rep(0,n[3]-y[3]),
rep(1,y[4]),rep(0,n[4]-y[4]),
rep(1,y[5]),rep(0,n[5]-y[5]),
rep(1,y[6]),rep(0,n[6]-y[6]))
xs = c(
rep(x[1],n[1]),
rep(x[2],n[2]),
rep(x[3],n[3]),
rep(x[4],n[4]),
rep(x[5],n[5]),
rep(x[6],n[6]))

# Generalized linear model (GLM)
run1 = glm(ys~xs,family=binomial(link="logit"))
run2 = glm(ys~xs,family=binomial(link="probit"))
run3 = glm(ys~xs,family=binomial(link="cloglog"))

par(mfrow=c(1,1))
xx = seq(min(xs),max(xs),length=1000)
plot(xs+rnorm(ns,0,0.1),ys+rnorm(ns,0,0.025),
     ylab="y=1 se o inseto morreu",xlab="Dose de rotenona",axes=FALSE)
axis(2);box();axis(1,at=x)
lines(xx,1/(1+exp(-run1$coef[1]-run1$coef[2]*xx)),col=2,lwd=2)
lines(xx,pnorm(run2$coef[1]+run2$coef[2]*xx),col=3,lwd=2)
lines(xx,1-exp(-exp(run3$coef[1]+run3$coef[2]*xx)),col=4,lwd=2)
legend("topleft",legend=c("logit","probit","cloglog"),col=2:4,lwd=2,bty="n")
abline(h=0.5,lty=2)
abline(h=0.9,lty=2)

# Bayesian GLM
N = 100
alphas = seq(-4.5,-2,length=N)
betas  = seq(0.3,0.9,length=N)
like   = matrix(0,N,N)
for (i in 1:N)
for (j in 1:N){
  theta = 1/(1+exp(-alphas[i]-betas[j]*xs))
  like[i,j] = sum(ys*log(theta))+sum((1-ys)*log(1-theta))
}
like = exp(like-max(like))

contour(alphas,betas,like,xlab=expression(alpha),ylab=expression(beta),
        drawlabels=FALSE)
points(run$coef[1],run$coef[2],col=2,pch=16)

# SIR-based posterior inference
M     = 10000
alpha = runif(M,-5,-1)
beta  = runif(M,0.2,1.0)
w     = rep(0,M)
for (i in 1:M){
  theta = 1/(1+exp(-alpha[i]-beta[i]*xs))
  w[i]  = sum(ys*log(theta))+sum((1-ys)*log(1-theta))
}
w = exp(w-max(w))
ind = sample(1:M,size=M,replace=TRUE,prob=w)
alpha1 = alpha[ind]
beta1  = beta[ind]

thetas = matrix(0,M,1000)
for (i in 1:1000)
 thetas[,i] = 1/(1+exp(-alpha1-beta1*xx[i]))
qtheta = t(apply(thetas,2,quantile,c(0.05,0.5,0.95)))


par(mfrow=c(2,2))
plot(alpha,beta,xlab=expression(alpha),ylab=expression(beta))
points(alpha1,beta1,col=2)
hist(alpha1,prob=TRUE,xlab="",main=expression(alpha))
hist(beta1,prob=TRUE,xlab="",main=expression(beta))


plot(xs+rnorm(ns,0,0.1),ys+rnorm(ns,0,0.025),
     ylab="y=1 se o inseto morreu",xlab="Dose de rotenona",axes=FALSE)
axis(2);box();axis(1,at=x)
lines(xx,1/(1+exp(-run$coef[1]-run$coef[2]*xx)),col=2,lwd=2)
lines(xx,qtheta[,1],col=3,lwd=2)
lines(xx,qtheta[,2],col=2,lwd=2)
lines(xx,qtheta[,3],col=3,lwd=2)

ind1 = (qtheta[,1]<0.501)&(qtheta[,1]>0.499)
ind2 = (qtheta[,2]<0.501)&(qtheta[,2]>0.499)
ind3 = (qtheta[,3]<0.501)&(qtheta[,3]>0.499)
segments(xx[ind3][1],0.5,xx[ind1][1],0.5,lwd=3)
interval1 = c(xx[ind3][1],xx[ind2][1],xx[ind1][1])

ind1 = (qtheta[,1]<0.901)&(qtheta[,1]>0.899)
ind2 = (qtheta[,2]<0.901)&(qtheta[,2]>0.899)
ind3 = (qtheta[,3]<0.901)&(qtheta[,3]>0.899)
segments(xx[ind3][1],0.9,xx[ind1][1],0.9,lwd=3)
interval2 = c(xx[ind3][1],xx[ind2][1],xx[ind1][1])

interval1
interval2