envel.empar1 <- function(x1=stop("Informe a variavel explicativa binaria/dummy de y=1!\n"),x0=stop("Informe a variavel explicativa binaria/dummy de y=0!\n"),iden=0,nome=NULL,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 quando se usa a distribuição
# binomial (com ligação logito) para o caso de emparelhamentos 1 a 1 e apenas uma variavel explicativa binaria/dummy,
# ou seja, na verdade o ajuste da binomial é apenas um artifício e o modelo de interesse é o condicional que é a
# hipergeométrica não-central.
# Para maiores detalhes, ver Paula (2003), pág.167.
#
# 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:
# x1: deve ser um vetor de valores 1 ou 0 da variavel explicativa para y=1 (casos);
# x0: deve ser um vetor de valores 1 ou 0 da variavel explicativa para y=0 (controles);
# 
# 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";
# 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.
#
# O modelo que será ajustado é logito(sucesso)=z (sem intercepto), onde z=x1-x0
#
# 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.empar1(x1=exposto,x0=naoexposto,sim=1000,conf=.95,maxit=50)
# envel.empar1(ajuste,res="C")
#

if(any(x1!=1 & x1!=0)) {
	stop(paste("\nA variavel explicativa x1 precisa ser binaria/dummy, isto e, precisam ter apenas valores 0 ou 1!\n"))
}
if(any(x0!=1 & x0!=0)) {
	stop(paste("\nA variavel explicativa x1 precisa ser binaria/dummy, isto e, precisam ter apenas valores 0 ou 1!\n"))
}

n<-length(x1)

if(is.null(nome)) {
	nome<-1:n
}
if(length(iden)>1) {
	iden<--1
}

if(length(x0)!=n) {
	stop(paste("\nO comprimento das variaveis explicativas nao e o mesmo!\n"))
}

alfa<-(1-conf)/2
X <- x1-x0
y<-rep(1,n)
modelo <- glm(y ~ X-1,family=binomial)
p <- 1
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")

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,1,m),max=pbinom(y,1,m)) )
} 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=="W") {
				tipo<-"Resíduo de Williams"
				r<-sign(y-m)*sqrt((1-h)*(( resid(modelo,type="deviance")/sqrt(1-h) )^2)+(h*( resid(modelo,type="pearson")/sqrt(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"))
				}
			}
		}
	}
}

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) {
	x1g<-x1
	x0g<-x0
	ger <- ifelse(runif(n)-mu>=0,1,0)
	x1g[x1!=x0]<-ger[x1!=x0]
	x0g[x1!=x0]<-1-ger[x1!=x0]
	Xg<-x1g-x0g
	fit <- glm(y ~ Xg-1,family=binomial,maxit=maxit,start=pm)
	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")
	e[,i] <- 
	sort( if(res=="Q") {
		qnorm( runif(n=n,min=pbinom(y-1,1,m),max=pbinom(y,1,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-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")
}