envel.nb <- function(modelo=fit.model,iden=0,nome=seq(along = model.matrix(modelo)[,1]),sim=100,conf=.90,res="D",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 negativa
# (utilizando a biblioteca MASS).
#
# 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.
#
# 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 e "C" distância de Cook. A opção padrão é a "D";
# 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]
# 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]
# SVETLIZA, C. F. (2002). Modelos Não-Lineares com Resposta Binomial Negativa. Tese de Doutorado. IME-USP, São Paulo.
#
# Exemplos:
# envel.nb(ajuste,sim=1000,conf=.95,maxit=50)
# envel.nb(ajuste,res="C")
#

if(class(modelo)[1] != "negbin") {
	stop(paste("\nA classe do objeto deveria ser binomial negativa e nao ",class(modelo)[1],"!!!\n"))
}
library("MASS")

alfa<-(1-conf)/2
X <- model.matrix(modelo)
n <- nrow(X)
p <- ncol(X)
fi <- modelo$theta
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<-modelo$y

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"
	p<-fi/(m+fi)
	a<-ifelse( y>0 , pbeta(p,fi,pmax(y,1)) , 0)
	b<-pbeta(p,fi,y+1)
	r<-qnorm( runif(n=n,min=a,max=b) )
} else {
	if(res=="D") {
		tipo<-"Resíduo Componente do Desvio"
		r<-resid(modelo,type="deviance")/sqrt(1-h)
	} else {
		if(res=="P") {
			tipo<-"Resíduo de Pearson Padronizado"
			r<-resid(modelo,type="pearson")/sqrt(1-h)
		} 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)
}
fii<-fi
mu<-m
for(i in 1:sim) {
	resp <- rnegbin(n,mu,fii)
	if ( (is.null(version$language) == T && link == "Log: log(mu)") | (is.null(version$language) == F && link == "log") ) {
		fit <- glm.nb(resp ~ X-1,init.theta=fii,maxit=maxit,start=pm)
	} else {
		if ( (is.null(version$language) == T && link == "Square Root: sqrt(mu)") | (is.null(version$language) == F && link == "sqrt") ) {
			fit <- glm.nb(resp ~ X-1,init.theta=fii,link=sqrt,maxit=maxit,start=pm)
		} else {
			if ( (is.null(version$language) == T && link == "Identity: mu") | (is.null(version$language) == F && link == "identity") ) {
				fit <- glm.nb(resp ~ X-1,init.theta=fii,link=identity,maxit=maxit,start=pm)
			} else {
				stop(paste("\nEsta funcao so aceita as ligacoes: canonica (log), raiz quadrada e identidade!!!\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")
	if(res=="Q") {
		y<-fit$y
		fi<-fit$theta
		p<-fi/(m+fi)
		a<-ifelse( y>0 , pbeta(p,fi,pmax(y,1)) , 0)
		b<-pbeta(p,fi,y+1)
		rq<-runif(n=n,min=a,max=b)
	}
	e[,i] <- 
	sort( if(res=="Q") {
		qnorm( rq )
	} 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=="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")
}