diag.norm <- function(modelo=fit.model,iden=c(0,0,0,0,0,0),nome=seq(along = model.matrix(modelo)[,1])) {

#
# Descrição e detalhes:
# A saída terá seis gráficos:
# 1º) Influência na Locação. O gráfico feito é das distâncias de Cook contra os valores ajustados. Utilizou-se o critério
#     de destacar observações maiores do que duas vezes a média de todas as distâncias obtidas;
# 2º) Influência Locação/Escala. A medida C, que é um aperfeiçoamento do DFFIT e também é conhecida como distância de Cook
#     modificada, foi utilizada para medir a influência das observações nos parâmetros de locação e escala. O critério
#     foi o de destacar observações maiores do que duas vezes a média de todas as distâncias obtidas;
# 3º) Influência Local. A influência local consiste em procurar pontos que sob pequenas perturbações causam variações
#     muito grandes nos resultados. O dmax é o autovetor que corresponde ao maior autovalor da matriz do processo de
#     perturbações. Para maiores detalhes veja Paula (2003, págs.50-54 e 65-66). O critério foi o de destacar observações
#     maiores do que duas vezes a média de todos os dmax's;
# 4º) Pontos Alavanca. A matriz chapéu H=X%*%solve(t(X)%*%X)%*%t(X) é a matriz de projeção ortogonal de vetores no
#     subespaço gerado pelas colunas da matrix X. Os pontos remotos nesse subespaço costumam ser considerados alavanca
#     (leverage), por exercer uma forte influência no seu valor ajustado. Ou seja, esses pontos tem um perfil diferente
#     dos demais com relação às variáveis explicativas. Ao fazer predições para um determinado vetor x, pode-se também
#     obter a medida h para esse valor, através de h=t(x)%*%solve(t(X)%*%X)%*%x. Caso esse valor seja grande com relação
#     aos pontos utilizados na estimação do modelo, isso é um indício em que a combinação de valores de x é uma
#     extrapolação, mesmo que o valor separadamente de cada variável esteja dentro dos limites em que o modelo abrange;
# 5º) Pontos Aberrantes. Um ponto é aberrante (discrepante, outlier) se o seu valor estiver mal ajustado pelo modelo.
#     Como os resíduos utilizados tem uma distribuição t-Student de n-p-1 graus de liberdade, em que n é o número de
#     observações e p o número de parâmetros, então adicionamos linhas tracejadas nos quantis 2.5% e 97.5%. Com isso
#     esperamos que cerca de 5% dos pontos possam estar um pouco fora desses limites. Esse gráfico serve como indicação
#     para detectar valores aberrantes marginalmente. Se o objetivo for detectar valores conjuntamente aberrantes
#     deve-se construir o gráfico de envelopes ou utilizar critérios de comparações múltiplas, como o de Bonferroni
#     que consistiria em utilizar os quantis 2.5%/n e 1-2.5%/n, uma vez que estamos fazendo n comparações. Para ter
#     uma idéia, pode-se obter esses valores através dos comandos
#	qt(.025/sum(summary(modelo)$df[1:2]),summary(modelo)$df[2]-1);
#	qt(1-.025/sum(summary(modelo)$df[1:2]),summary(modelo)$df[2]-1);
# 6º) Função de Variância. McCullagh e Nelder (1989, pág.400) sugere o gráfico dos resíduos absolutos contra os valores
#     ajustados para checar se a função de variância adotada é adequada. O padrão esperado é de não encontrarmos nenhuma
#     tendência. Funções de variância erradas irão resultar em tendências dos resíduos com a média. Tendências positivas
#     indicam que a função de variância está crescendo muito devagar com a média, então deve-se aumentar a potência (no
#     caso de uma função de variância da família potência). Uma linha suavizada pelo método lowess robusto é adicionada
#     para ajudar na procura de tendências.
#
# Os dados devem estar disponíveis pelo comando attach( ).
#
# Argumentos obrigatórios:
# modelo: deve-se informar o objeto onde está o ajuste do modelo normal linear, 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 em cada gráfico. O vetor deve conter 6
#	posições de números inteiros. A ordem que deve ser informada é a mesma em que os gráficos são feitos. Os
#	componentes do vetor iguais a 0 indicam que não se quer que identifique pontos, se for um inteiro positivo irá
#	automaticamente nos gráficos respectivos permitir que identifiquemos o número de pontos solicitados e qualquer
#	outro valor (negativo ou decimal) parar nos gráficos e solicitar que especifiquemos o número de pontos a ser
#	destacado. O padrão é c(0,0,0,0,0,0) caso não se entre com nada e c(-1,-1,-1,-1,-1,-1) caso se entre
#	com qualquer coisa que não seja um vetor de 6 posições, como por ex.-1;
# nome: esse argumento só é utilizado caso algum dos componentes do vetor da opção iden não seja 0. Caso não seja
#	informado nada, os pontos identificados serão os números da ordem em que estão no banco de dados (í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.
#
# A função retorna os seguintes valores: ResPearsonStd, Di, Ci, Dmax e h.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Referências:
# MCCULLAGH, P. e NELDER, J. A. (1989). Generalized Linear Models. 2ª ed. Chapman and Hall, London.
# NETER, J., KUTNER, M. H., NACHTSHEIM, C. J. and WASSERMAN, W. (1996). Applied Linear Statistical Models. 4ª ed.
#    Mc Graw Hill, Boston.
# 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:
# diag.norm(ajuste,iden=c(1,5,2,4,3,0),nome=estados)
# diag.norm(ajuste,iden=-1)
#

if( class(modelo)[1]=="lm" || (class(modelo)[1]=="glm" && (modelo$family[[1]]=="Gaussian" | modelo$family[[1]]=="gaussian")) ) {

} else {
	stop(paste("\nA classe do objeto deveria ser lm ou glm (com distribuicao gaussian) !!!"))
}

if(length(iden)<6) {
	iden<-c(-1,-1,-1,-1,-1,-1)
}

X <- model.matrix(modelo)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
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

lms <- summary(modelo)
s <- lms$sigma
r <- resid(modelo)
ts <- r/(s*sqrt(1-h))
di <- (1/p)*(h/(1-h))*(ts^2)
si <- lm.influence(modelo)$sigma
tsi <- r/(si*sqrt(1-h))
#dff <- sqrt(h/(1-h))*abs(tsi) #DFFIT
ci <- sqrt( ((n-p)*h) / (p*(1-h)) )*abs(tsi) #aperfeiçoamento do DFFIT
A <- diag(r)%*%H%*%diag(r)
dmax <- abs(eigen(A)$vec[,1]/sqrt(eigen(A)$val[1]))
m <- fitted(modelo)

par(mfrow=c(2,3))

plot(m,di,xlab="Valor Ajustado", ylab="Distância de Cook",main="Influência na Locação", ylim=c(0,max(di,2*mean(di))), pch=16)
abline(2*mean(di),0,lty=2)
while ( (!is.numeric(iden[1])) || (round(iden[1],0) != iden[1]) || (iden[1] < 0) ) {
	cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
	out <- readline()
	iden[1]<-as.numeric(out)
}
if(iden[1]>0) {identify(m,di,n=iden[1],labels=nome)}

plot(m,ci,xlab="Valor Ajustado", ylab="Distância de Cook Modificada",main="Influência Locação/Escala", ylim=c(0,max(ci,2*mean(ci))), pch=16)
abline(2*mean(ci),0,lty=2)
while ( (!is.numeric(iden[2])) || (round(iden[2],0) != iden[2]) || (iden[2] < 0) ) {
	cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
	out <- readline()
	iden[2]<-as.numeric(out)
}
if(iden[2]>0) {identify(m,ci,n=iden[2],labels=nome)}

plot(m,dmax,xlab="Valor Ajustado", ylab="dmax",main="Influência Local", ylim=c(0,max(dmax,2*mean(dmax))), pch=16)
abline(2*mean(dmax),0,lty=2)
while ( (!is.numeric(iden[3])) || (round(iden[3],0) != iden[3]) || (iden[3] < 0) ) {
	cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
	out <- readline()
	iden[3]<-as.numeric(out)
}
if(iden[3]>0) {identify(m,dmax,n=iden[3],labels=nome)}

plot(m,h,xlab="Valor Ajustado", ylab="Medida h",main="Pontos Alavanca", ylim=c(0,max(h,2*p/n)), pch=16)
abline(2*p/n,0,lty=2)
while ( (!is.numeric(iden[4])) || (round(iden[4],0) != iden[4]) || (iden[4] < 0) ) {
	cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
	out <- readline()
	iden[4]<-as.numeric(out)
}
if(iden[4]>0) {identify(m,h,n=iden[4],labels=nome)}

plot(m,tsi,xlab="Valor Ajustado", ylab="Resíduo Padronizado",main="Pontos Aberrantes", ylim=c(min(tsi)-1,max(tsi)+1), pch=16)
abline(qt(.025,n-p-1),0,lty=2)
abline(qt(1-.025,n-p-1),0,lty=2)
while ( (!is.numeric(iden[5])) || (round(iden[5],0) != iden[5]) || (iden[5] < 0) ) {
	cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
	out <- readline()
	iden[5]<-as.numeric(out)
}
if(iden[5]>0) {identify(m,tsi,n=iden[5],labels=nome)}

plot(m,abs(tsi),xlab="Valor Ajustado", ylab="Resíduo Padronizado Absoluto",main="Função de Variância", pch=16)
lines(lowess(m,abs(tsi)))
while ( (!is.numeric(iden[6])) || (round(iden[6],0) != iden[6]) || (iden[6] < 0) ) {
	cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
	out <- readline()
	iden[6]<-as.numeric(out)
}
if(iden[6]>0) {identify(m,abs(tsi),n=iden[6],labels=nome)}

par(mfrow=c(1,1))
list(ResPearsonStd=tsi,Di=di,Ci=ci,Dmax=dmax,h=h)
}