#source("E:\\windows\\Unicamp\\Disciplinas\\1_semestre_2014\\ME 714\\Programas\\envel_binoM.txt")

#---------------------------------------------------------------#
# 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  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
#     attach(dados)
#     source("envel_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áfico após
# a saída use title("..."). Para outras ligações substituir
# o termo family=binomial por family=binomial(link=probit) ou por
# family=binomial(link=cloglog).
#---------------------------------------------------------------#
#
par(mfrow=c(1,1))
X <- model.matrix(fit.model)
n <- nrow(X)
p <- ncol(X)
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 <- resid(fit.model,type="deviance")/sqrt(1-h)
e <- matrix(0,n,100)
#
for(i in 1:100){
nresp <- rbinom(n,m,fitted(fit.model))
fit <- glm(cbind(nresp,m-nresp) ~ X, family=binomial(link="logit"))
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(n)
e2 <- numeric(n)
#
for(i in 1:n){
  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="Percentil da N(0,1)",
ylab="Resíduo Componente do Desvio", ylim=faixa, pch=16,main="Gráfico de quantil-quantil normal")
#
par(new=T)
#
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1,main="")
par(new=T)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1,main="")
par(new=T)
qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2,main="")
#---------------------------------------------------------------#