envel.bino <- function(modelo=fit.model,iden=0,nome=seq(along = model.matrix(modelo)[,1]),sim=100,conf=.90,res="D",quasi=F,quad=T,maxit=20) {

#
# Descrição e detalhes:
# A saída será o gráfico de probabilidade normal com envelopes simulados para um ajuste da distribuição binomial ou
# "quase-binomial" (com ou sem réplicas, com ligação logito, probito, complementar log-log ou Aranda-Ordaz).
#
# A opção res="C" faz o gráfico de probabilidade meio-normal com envelopes simulados utilizando a distância de Cook,
# possibilitando a detecção de pontos simultaneamente aberrantes e/ou influentes.
#
# Vale ressaltar que na "quase-binomial" o método quase-verossimilhança assume apenas a função da variância em relação
# à média, e portanto o ajuste corresponderia a qualquer distribuição que tiver tal função de variância e não apenas a
# binomial. Na geração de envelopes estamos assumindo que a distribuição é a binomial, mas os resíduos do ajuste são
# corrigidos pela super ou subdispersão. Se ajustar pelo R é necessário usar a familia quasibinomial ao invés da
# quasi(variance="mu(1-mu)") pois caso contrário a função não consegue perceber que o ajuste é da distribuição
# 'quase-binomial' e parará a execução.
#
# Atenção: a função não funcionará corretamente se o ajuste possuir offsets! Neste caso é preciso adaptá-la como foi
# feito na função envel.pois
#
# Os dados devem estar disponíveis pelo comando attach( ).
#
# Argumentos obrigatórios:
# modelo: deve-se informar o objeto onde está o ajuste do modelo, caso não seja informado, a função procurará
# 	  o ajuste no objeto fit.model;
# 
# Argumentos opcionais:
# iden: caso deseje, informe o número de observações que irá querer destacar. O padrão é não destacar ninguém (iden=0).
#	Qualquer valor que não seja um inteiro positivo (por ex., negativo ou decimal) fará com que a função pergunte
#	o número de pontos após a execução;
# nome: esse argumento só é utilizado caso seja destacado algum ponto no gráfico. Caso não seja informado nada, os pontos
#	identificados serão os números da ordem em que estão no banco de dados (os índices). Caso se queira, pode-se
#	informar um vetor de nomes ou de identificações alternativas. Obrigatoriamente esse vetor deve ter o mesmo
#	comprimento do banco de dados;
# sim: número de simulações para gerar a banda de confiança. Atkinson sugere um mínimo de 20 simulações.
#      O padrão é de 100;
# conf: nível de confiança do envelope. O padrão é de 90%;
# res: permite-se a escolha dos resíduos. As opções dos resíduos são: "Q" quantil (ver Dunn e Smyth, 1996), "D" componente
#      do desvio, "P" Pearson padronizado, "W" Williams e "C" distância de Cook. A opção padrão é a "D";
# quasi: se estiver usando a ligação de Aranda-Ordaz para a "quase-binomial" é necessário informar a opção quasi=T pois
#	 para esse caso a função não consegue detectar que o ajuste é da "quase-binomial";
# quad: o padrão (quad=T, True) faz um gráfico quadrado, enquanto quad=F (False) faz um gráfico utilizando a área máxima
#       disponível;
# maxit: essa opção é utilizada nos ajustes de cada simulação e indica o máximo de iterações permitidas nos ajustes.
#	 O padrão é maxit=20.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Referências:
# DUNN, K. P., and SMYTH, G. K. (1996). Randomized quantile residuals. J. Comput. Graph. Statist. 5, 1-10
#    [http://www.statsci.org/smyth/pubs/residual.html e http://www.statsci.org/smyth/pubs/residual.ps]
# HOSMER, D. W. e LEMESHOW, S. (2000). Applied Logistic Regression. John Wiley & Sons, New York.
# MCCULLAGH, P. e NELDER, J. A. (1989). Generalized Linear Models. 2ª ed. Chapman and Hall, London.
# PAULA, G. A. (2003). Modelos de Regressão com apoio computacional. IME-USP, São Paulo. [Não publicado,
#    disponível em http://www.ime.usp.br/~giapaula/Book.pdf]
#
# Exemplos:
# envel.bino(ajuste,sim=1000,conf=.95,maxit=50)
# envel.bino(ajuste,res="C")
#

if(class(modelo)[1] != "glm") {
	stop(paste("\nA classe do objeto deveria ser glm e nao ",class(modelo),"!!!\n"))
}
if(modelo$family[[1]] != "Binomial" & modelo$family[[1]] != "binomial" & modelo$family[[1]] != "Aranda") {
	if(modelo$family[[1]] == "quasibinomial" || ( modelo$family[[1]] == "Quasi-likelihood" & modelo$family[[3]] == "Binomial: mu(1-mu)") ) {
		quasi<-T
	} else {
		stop(paste("\nA familia do objeto deveria ser Binomial e nao ",modelo$family[[1]],"!!!\n"))
	}
}

if(length(iden)>1) {
	iden<--1
}

alfa<-(1-conf)/2
X <- model.matrix(modelo)
n <- nrow(X)
p <- ncol(X)
if(is.null(modelo$prior.weights)) {
	ntot<-rep(1,n)
} else {
	ntot<-modelo$prior.weights
}
w <- modelo$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)

#para evitar divisão por 0 ao studentizar os residuos, mas tentando manter o valor exagerado da alavanca
h[round(h,15)==1]<-0.999999999999999

m<-predict(modelo,type="response")
y<-round(modelo$y*ntot,0)

fi <- 1
if(quasi==T) {
	fi <- 1/summary(modelo)$dispersion
}

if(res=="Q") {
	cat("Ao utilizar o residuo Quantil para distribuicoes discretas, sugere-se plotar pelo menos 4 graficos para evitar conclusoes viesadas pela aleatoriedade que esta sendo incluida.\n")
	tipo<-"Resíduo Quantil"
	r<-qnorm( runif(n=n,min=pbinom(y-1,ntot,m),max=pbinom(y,ntot,m)) *sqrt(fi) )
} else {
	if(res=="D") {
		tipo<-"Resíduo Componente do Desvio"
		r<-resid(modelo,type="deviance")*sqrt(fi*(1-h))
	} else {
		if(res=="P") {
			tipo<-"Resíduo de Pearson Padronizado"
			r<-resid(modelo,type="pearson")*sqrt(fi*(1-h))
		} else {
			if(res=="W") {
				tipo<-"Resíduo de Williams"
				r<-sign((y/ntot)-m)*sqrt((1-h)*(( resid(modelo,type="deviance")*sqrt(fi/(1-h)) )^2)+(h*( resid(modelo,type="pearson")*sqrt(fi/(1-h)) )^2))
			} else {
				if(res=="C") {
					tipo<-"Distância de Cook"
					r<-(h/((1-h)*p))*((resid(modelo,type="pearson")/sqrt(1-h))^2)
				} else {
					stop(paste("\nVoce nao escolheu corretamente um dos residuos disponiveis!!!\n"))
				}
			}
		}
	}
}

link<-modelo$family[[2]]

e <- matrix(0,n,sim)
e1 <- numeric(n)
e2 <- numeric(n)

if (is.null(version$language) == T) {
	#No S-Plus, a opção start é para entrar com o preditor linear
	pm<-predict(modelo)
} else {
	#No R, a opção start é para entrar com os coeficientes
	pm<-coef(modelo)
}
mu<-m
for(i in 1:sim) {
	resp <- rbinom(n,ntot,mu)
	if ( (is.null(version$language) == T && link == "Logit: log(mu/(1 - mu))") | (is.null(version$language) == F && link == "logit") ) {
		fit <- glm(cbind(resp,ntot-resp) ~ X-1,family=binomial,maxit=maxit,start=pm)
	} else {
		if ( (is.null(version$language) == T && link == "Probit: qnorm(mu)") | (is.null(version$language) == F && link == "probit") ) {
			fit <- glm(cbind(resp,ntot-resp) ~ X-1,family=binomial(link=probit),maxit=maxit,start=pm)
		} else {
			if ( (is.null(version$language) == T && link == "Complementary Log: log( - log(1 - mu))") | (is.null(version$language) == F && link == "cloglog") ) {
				fit <- glm(cbind(resp,ntot-resp) ~ X-1,family=binomial(link=cloglog),maxit=maxit,start=pm)
			} else {
				if ( modelo$family[[1]] == "Aranda" ) {
					fit <- do.call("glm",list(cbind(resp,ntot-resp) ~ X-1,family=aranda(as.numeric(link)),maxit=maxit,start=pm))
				} else {
					stop(paste("\nEsta funcao so aceita as ligacoes: logito, probito, complementar log-log e Aranda-Ordaz !!!\nLigacao ",link," desconhecida!!!\n"))
				}
			}
		}
	}
	w <- fit$weights
	W <- diag(w)
	H <- solve(t(X)%*%W%*%X)
	H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
	h <- diag(H)
	h[round(h,15)==1]<-0.999999999999999
	m <- predict(fit,type="response")
	y <- round(fit$y*ntot,0)
	e[,i] <- 
	sort( if(res=="Q") {
		qnorm( runif(n=n,min=pbinom(y-1,ntot,m),max=pbinom(y,ntot,m)) )
	} else {
		if(res=="D") {
			resid(fit,type="deviance")/sqrt(1-h)
		} else {
			if(res=="P") {
				resid(fit,type="pearson")/sqrt(1-h)
			} else {
				if(res=="W") {
					sign((y/ntot)-m)*sqrt((1-h)*(( resid(fit,type="deviance")/sqrt(1-h) )^2)+(h*( resid(fit,type="pearson")/sqrt(1-h) )^2))
				} else {
					if(res=="C") {
						(h/((1-h)*p))*((resid(fit,type="pearson")/sqrt(1-h))^2)
					} else {
						stop(paste("\nVoce nao escolheu corretamente um dos residuos disponiveis!!!\n"))
					}
				}
			}
		}
	})
}

for(i in 1:n) {
	eo <- sort(e[i,])
	e1[i] <- quantile(eo,alfa)
	e2[i] <- quantile(eo,1-alfa)
}

med <- apply(e,1,median)

if(quad==T) {
	par(pty="s")
}
if(res=="C") {
	#Segundo McCullagh e Nelder (1989, pág.407) e Paula (2003, pág.57) deve-se usar qnorm((n+1:n+.5)/(2*n+1.125))
	#Segundo Neter et alli (1996, pág.597) deve-se usar qnorm((n+1:n-.125)/(2*n+0.5))
	qq<-qnorm((n+1:n+.5)/(2*n+1.125))
	plot(qq,sort(r),xlab="Quantil Meio-Normal",ylab=tipo, ylim=range(r,e1,e2), pch=16)
} else {
	qq<-qnorm((1:n-.375)/(n+.25))
	plot(qq,sort(r),xlab="Quantil da Normal Padrão",ylab=tipo, ylim=range(r,e1,e2), pch=16)
}
lines(qq,e1,lty=1)
lines(qq,e2,lty=1)
lines(qq,med,lty=2)
nome<-nome[order(r)]
r<-sort(r)
while ( (!is.numeric(iden)) || (round(iden,0) != iden) || (iden < 0) ) {
	cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
	out <- readline()
	iden<-as.numeric(out)
}
if(iden>0) {identify(qq,r,n=iden,labels=nome)}
if(quad==T) {
	par(pty="m")
}
cat("Banda de ",conf*100,"% de confianca, obtida por ",sim," simulacoes.\n")
}