envel.mult <- function(modelo=fit.model,iden=NULL,nome=seq(along = model.matrix(modelo)[,1]),sim=100,conf=.90,res="D",tit=T,quad=T,maxit=100) {

#
# Descrição e detalhes:
# A saída será o gráfico de probabilidade normal com envelopes simulados para um ajuste da distribuição multinomial
# (através da biblioteca nnet). Para o resíduo componente do desvio, a saída é apenas 1 gráfico de probabilidade
# meio-normal, enquanto que para o resíduo de Pearson a saída será J-1 gráficos de probabilidade normal,
# sendo J o número de possíveis respostas da multinomial.
#
# 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.
#
# Note que a função multinom da biblioteca nnet faz o ajuste tanto considerando apenas um vetor de respostas com o nível
# de cada resposta, quanto utilizando uma matriz de respostas com J colunas. A diferença é que utilizando apenas
# um vetor, consideramos as respostas como um modelo Bernoulli (multivariado), enquanto que ao utilizar uma matriz
# consideramos um modelo Multinomial. O modelo Bernoulli (uni/multivariado) é chamado pro alguns autores de
# "Binomial/Multinomial sem réplicas", enquanto o Multinomial/Binomial de "Binomial/Multinomial com réplicas". Outros
# autores ao invés de usar os termos com e sem réplicas, usam dados desagrupados e agrupados.
#
# Nota: essa função foi programada para tentar recuperar a(o) matriz(vetor) resposta. No entanto, dependendo da forma
#	como o modelo foi ajustado pode ser que ocorra algum problema. Nesse caso, talvez seja necessário reajustar
#	o modelo de uma forma mais adequada. Normalmente, o problema costuma ocorrer apenas quando o ajuste é feito
#	utilizando alguma operação na variável resposta, como por ex. lrm(cbind(resp1,resp2,resp3) ~ ...). Se primeiro
#	for criado uma variável resp<-cbind(resp1,resp2,resp3), então o ajuste lrm(resp ~ ...) não deve causar problemas.
#
# 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 feito pela biblioteca nnet, 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. Para res="P" e res="C" deve-se informar um vetor com J-1 posições;
# 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: "D" componente do desvio,
#      "P" Pearson padronizado e "C" distância de Cook. A opção padrão é a "D";
# tit: como padrão (tit=T, True) o programa coloca o titulo de cada submodelo caso seja escolhido res="P" ou res="C". Se
#      não desejar que seja impresso nada como título, deve-se utilizar tit=F (False). Essa opção não é avaliada quando
#      res="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=100.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Referências:
# 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]
# SILVA, J. A. P. (2003). Métodos de Diagnóstico em Modelos Trinomiais. Dissertação de Mestrado. IME-USP, São Paulo.
# THOMPSON, L. A. (2003). S-PLUS (and R) Manual to Accompany Agresti’s Categorical Data Analysis (2002). 2ª ed.
#    [Não publicado, disponível em http://math.cl.uh.edu/~thompsonla/Splusdiscrete2.pdf]
#
# Exemplos:
# envel.mult(ajuste,sim=1000,conf=.95,maxit=50)
# envel.mult(ajuste,res="C")
#

if(class(modelo)[1] != "multinom") {
	stop(paste("\nA classe do objeto deveria ser multinom (biblioteca nnet) e nao ",class(modelo),"!!!\n"))
}

alfa<-(1-conf)/2
X <- model.matrix(modelo)
n <- nrow(X)
p <- ncol(X)
ntot<-modelo$weights

mu<-modelo$fitted.values
y<-do.call("as.array",list(as.name(dimnames(attr(modelo$terms,"factors"))[[1]][1])))
niveis<-modelo$lab
J<-length(niveis)

if(res=="D") {
	if(is.null(iden)) {
		iden<-0
	}
} else {
	if(is.null(iden)) {
		iden<-rep(0,J-1)
	} else {
		if(length(iden)<(J-1)) {
			iden<-rep(-1,J-1)
		}
	}
}

if(is.na(dim(y)[2])) {
	y2<-matrix(0,n,J)
	for(j in 1:J) {
		y2[y==niveis[j],j]<-1
	}
	y<-y2
}

#IO<-modelo$Hess
#if(is.null(IO)) {
#	modelo<-multinom(y~X-1,Hess=T,trace=F)
#	IO<-modelo$Hess
#}
#VarO<-solve(IO)

#S<-kronecker(diag(J-1),X)
V<-matrix(0,(J-1)*n,(J-1)*n)
for(j in 1:(J-1)) {
	for(k in 1:(J-1)) {
		if(j == k) {
			V[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]<-diag(mu[,j]*(1-mu[,j]))
		} else {
			V[(n*(j-1)+1):(n*j),(n*(k-1)+1):(n*k)]<-diag(-mu[,j]*mu[,k])
		}
	}
}
#IF<-t(S)%*%V%*%S
#VarF<-solve(IF)
#M<-S%*%VarF%*%t(S)%*%V

h<-matrix(0,n,J-1)
#m<-matrix(0,n,J-1)

for(j in 1:(J-1)) {
	VT<-V[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]
	h[,j]<-diag( sqrt(VT)%*%X%*%solve(t(X)%*%VT%*%X)%*%t(X)%*%sqrt(VT) )
#	m[,j]<-diag( M[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)] )
}
hm<-apply(h,1,mean)

if(res=="D") {
	tipo<-"Resíduo Componente do Desvio"
	r<-numeric(n)
	for(j in 1:J) {
		r<-r+ifelse( y[,j]==0,0,y[,j]*log(y[,j]/(ntot*mu[,j])) )
	}
	r<-sqrt(2*r/(1-hm))
	e <- matrix(0,n,sim)
} else {
	if(res=="P" | res=="C") {
		if(res=="C") {
			tipo<-"Distância de Cook"
		} else {
			tipo<-"Resíduo de Pearson Padronizado"
		}
		r<-matrix(0,n,J-1)
		for(j in 1:(J-1)) {
			r[,j]<-(y[,j]-ntot*mu[,j])/sqrt(ntot*mu[,j]*(1-mu[,j])*(1-h[,j]))
			if(res=="C") {
				r[,j]<-(h[,j]/((1-h[,j])*p))*(r[,j]^2)
			}
		}
		e <- array(0,dim=c(n,J-1,sim))
	} else {
		stop(paste("\nVoce nao escolheu corretamente um dos residuos disponiveis!!!\n"))
	}
}

#Funções para gerar vetores de multinomiais para o S-Plus e o R. Ambas foram retiradas do manual de Thompson (2003, pág.7)
if (is.null(version$language) == T) {
	#S-Plus
	rmultinomial<-function (n, p, rows = max(c(length(n), nrow(p)))) {
	n <- rep(n, length = rows)
	p <- p[rep(1:nrow(p), length = rows), , drop = FALSE]
	sapply(1:rows,function(i,n,p) {
		k <- length(p[i,])
		tabulate(sample(k, n[i], replace = TRUE, prob = p[i,]), nbins = k)
		},n=n,p=p)
	}
} else {
	#R
	rmultinomial<-function (n, p, rows = max(c(length(n), nrow(p)))) {
		rmultinomial.1 <- function(n, p) {
			k <- length(p)
			tabulate(sample(k, n, replace = TRUE, prob = p), nbins = k)
		}
		n <- rep(n, length = rows)
		p <- p[rep(1:nrow(p), length = rows), , drop = FALSE]
		sapply(1:rows, function(i) rmultinomial.1(n[i], p[i, ]))
		}
}
mu0<-mu
library(Matrix)
i<-1
tot<-0
while(i <= sim) {
	tot<-tot+1
	if(tot>(10*sim)) {
		stop(paste("\nA funcao descarta ajustes que produzem X'VX singular. Ja foram quase 10 vezes o numero de simulacoes solicitadas!\n"))
	}
	y <- t(rmultinomial(ntot,mu0))
	fiti <- multinom(y~X-1,maxit=maxit,trace=F) #Hess=T,
	mu<-fiti$fitted.values

	#S<-kronecker(diag(J-1),X)
	V<-matrix(0,(J-1)*n,(J-1)*n)
	for(j in 1:(J-1)) {
		for(k in 1:(J-1)) {
			if(j == k) {
				V[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]<-diag(mu[,j]*(1-mu[,j]))
			} else {
				V[(n*(j-1)+1):(n*j),(n*(k-1)+1):(n*k)]<-diag(-mu[,j]*mu[,k])
			}
		}
	}
	#IF<-t(S)%*%V%*%S
	#VarF<-solve(IF)
	#M<-S%*%VarF%*%t(S)%*%V

	h<-matrix(0,n,J-1)
	#m<-matrix(0,n,J-1)

	sing<-F	
	for(j in 1:(J-1)) {
		VT<-V[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]
		XVX<-t(X)%*%VT%*%X
		if (rcond.Matrix(XVX) < 10^(-15)) {
			sing<-T
		} else {
			h[,j]<-diag( sqrt(VT)%*%X%*%solve.Matrix(XVX)%*%t(X)%*%sqrt(VT) )
		#	m[,j]<-diag( M[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)] )
		}
	}

	if(sing==F) {
		hm<-apply(h,1,mean)

		if(res=="D") {
			for(j in 1:J) {
				e[,i]<-e[,i]+ifelse( y[,j]==0,0,y[,j]*log(y[,j]/(ntot*mu[,j])) )
			}
			e[,i]<-sqrt(2*e[,i]/(1-hm))
			e[,i]<-sort(e[,i])
		} else {
			if(res=="P" | res=="C") {
				for(j in 1:(J-1)) {
					e[,j,i]<-(y[,j]-ntot*mu[,j])/sqrt(ntot*mu[,j]*(1-mu[,j])*(1-h[,j]))
					if(res=="C") {
						e[,j,i]<-(h[,j]/((1-h[,j])*p))*(e[,j,i]^2)
					}
					e[,j,i]<-sort(e[,j,i])
				}
			} else {
				stop(paste("\nVoce nao escolheu corretamente um dos residuos disponiveis!!!\n"))
			}
		}
		i<-i+1
	}
}

if(quad==T) {
	par(pty="s")
}
if(res=="D") {
	e1 <- numeric(n)
	e2 <- numeric(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)
	#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=c(min(e1[1],r),max(e2[n],r)), 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)}
} else {
	if(res=="P" | res=="C") {
		if ((J-1)>2) {
			if ((J-1)>8) {	
				par(mfrow=c(3,ceiling((J-1)/3)))
			} else {
				par(mfrow=c(2,ceiling((J-1)/2)))
			}
		} else {
			par(mfrow=c(1,ceiling((J-1))))
		}
		for(j in 1:(J-1)) {
			e1 <- numeric(n)
			e2 <- numeric(n)
			for(i in 1:n) {
				eo <- sort(e[i,j,])
				e1[i] <- quantile(eo,alfa)
				e2[i] <- quantile(eo,1-alfa)
			}
			med <- apply(e[,j,],1,median)
			if(res=="P") {
				qq<-qnorm((1:n-.375)/(n+.25))
				plot(qq,sort(r[,j]),xlab="Quantil da Normal Padrão",ylab=tipo, ylim=c(min(e1[1],r[,j]),max(e2[n],r[,j])), pch=16)
			} else {
				#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[,j]),xlab="Quantil Meio-Normal",ylab=tipo, ylim=c(min(e1[1],r[,j]),max(e2[n],r[,j])), pch=16)
			}
			lines(qq,e1,lty=1)
			lines(qq,e2,lty=1)
			lines(qq,med,lty=2)

			nomej<-nome[order(r[,j])]
			r[,j]<-sort(r[,j])
			while ( (!is.numeric(iden[j])) || (round(iden[j],0) != iden[j]) || (iden[j] < 0) ) {
				cat("Digite o num.de pontos a ser identificado (0=nenhum) e <enter> para continuar\n")
				out <- readline()
				iden[j]<-as.numeric(out)
			}
			if(iden[j]>0) {identify(qq,r[,j],n=iden[j],labels=nomej)}

			if(tit==T) {
				title(niveis[j+1])
			}
		}
		par(mfrow=c(1,1))
	} else {
		stop(paste("\nVoce nao escolheu corretamente um dos residuos disponiveis!!!\n"))
	}
}
if(quad==T) {
	par(pty="m")
}

cat("Banda de ",conf*100,"% de confianca, obtida por ",sim," simulacoes.\n")
}