#------------------------------------------------------------#
# Para rodar este programa  deixe no objeto fit.model a saída 
# do  ajuste  da regressão com erros binomial e ligação logit. 
# Deixe tambem o vetor de réplicas no objeto ntot e  os dados
# disponíveis através do comando attach(...).  Depois  use  o 
# comando source(...) no R ou S-Plus para executar o programa.  
# A sequência de comandos é a seguinte:
#
#       fit.model <- ajuste
#       ntot <- replicas
#       attach(dados)
#       source("envelr_bino")
#
# A saída será o gráfico de  envelope para o resíduo
# componente do desvio padronizado. Para colocar um título
# no gráficos após a saída use title("..."). Para outras
# ligações substituiro termo family=binomial por
# family=binomial(link=probit) ou por 
# family=binomial(link=cloglog).
#------------------------------------------------------------#
par(mfrow=c(1,1))
X <- model.matrix(fit.model)
k <- nrow(X)
e <- matrix(0,k,100)
tot <- numeric(k)
w <- fit.model$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
td <- sort(resid(fit.model, type="deviance")/sqrt(1-h))
#
for(i in 1:100){
for(j in 1:k){
dif <- runif(ntot[j]) - fitted(fit.model)[j]
dif[dif>=0] <- 0
dif[dif<0]  <- 1
tot[j] <- sum(dif)}
xmat <- cbind(tot,ntot-tot)
fit <- glm(xmat ~ X, family=binomial(link=cloglog))
w <- fit$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
e[,i] <- sort(resid(fit, type="deviance")/sqrt(1-h))}
#
e1 <- numeric(k)
e2 <- numeric(k)
#
for(i in 1:k){
  eo <- sort(e[i,])
e1[i] <- (eo[2]+eo[3])/2
e2[i] <- (eo[97]+eo[98])/2}
#
med <- apply(e,1,mean)
faixa <- range(td,e1,e2)
par(pty="s")
qqnorm(td,xlab="Quantil da N(0,1)",
ylab="Componente do Desvio", ylim=faixa,pch=16, main="", cex=2, cex.axis=1.5, cex.lab=1.5)
#
par(new=TRUE)
qqnorm(e1,axes=FALSE,xlab="",ylab="",type="l",ylim=faixa,lty=1, main="",lwd=2)
par(new=TRUE)
qqnorm(e2,axes=FALSE,xlab="",ylab="", type="l",ylim=faixa,lty=1, main="", lwd=2)
par(new=TRUE)
qqnorm(med,axes=FALSE,xlab="", ylab="", type="l", ylim=faixa, lty=2, main="",lwd=2)
#------------------------------------------------------------#